I asked Wen-Shin Lee about this and she suggested that it is a natural for Prony's method. This is an area of expertise for her so I decided to attempt an implementation. I will remark that it is a close cousin to sparse polynomial interpolation methods of Ben-Or and Tiwari, and Zippel, and appears in numerous guises in the literature.

What follows is not necessarily best possible but it should be at least reasonably stable numerically (uses singular values under the hood, in two places). For signals with noise it can be modified to behave, e.g. by using a loose tolerance when computing the initial matrix rank.

I will illustrate with the example that has three frequencies (one of which is zero, that is, there's a constant term). The code and expositions I have looked at use time units of 1 so I have a timescale value to compensate for that. Also I show some intermediate results to give some idea of what is being calculated. But really an understanding of this will require some delving into the literature on this method.

func3freq = 0.3 + 3.5 Cos[2.6 t + .8] + 1.2 Cos[5.3 t - .4];
obstime = 80;
timescale = 1/10;
data3 = Table[func3freq, {t, 0, obstime, timescale}];
len = Length[data3];

Set up initial Prony matrix. We won't actually use it, we just need the rank. That tells us the numer of exponentials we are going to fit.

mat = Most[Partition[data3, Floor[len/2], 1]];
Dimensions[mat]
keep = Length[SingularValueList[mat]]

(* Out[321]= {401, 400} *)

(* Out[322]= 5 *)

Now that we know the rank we'll set up the actual matrix used to determine frequencies.

mat2 = Most[Partition[data3, keep, 1]];
Dimensions[mat2]
rhs = Drop[data3, keep];
Length[rhs]

(* Out[334]= {796, 5} *)

(* Out[336]= 796 *)

Find the Prony polynomial coefficients. Then solve for the complex exponentials as roots of the Prony polynomial.

soln = PseudoInverse[mat2].rhs
roots = xx /. NSolve[xx^keep - soln.xx^Range[0, keep - 1] == 0, xx]
Map[Norm, roots]
freqs = Log[roots]

(* Out[341]= {1., -4.6583940973, 8.99362652134, -8.99362652134, \
4.6583940973} *)

(* Out[342]= {0.862807070515 - 0.505533341205 I, 
 0.862807070515 + 0.505533341205 I, 0.966389978135 - 0.257080551892 I,
  0.966389978135 + 0.257080551892 I, 1.} *)

(* Out[343]= {1., 1., 1., 1., 1.} *)

(* Out[344]= {-1.02917674383*10^-13 - 0.53 I, -1.02917674383*10^-13 + 
  0.53 I, 3.16698056668*10^-13 - 0.26 I, 
 3.16698056668*10^-13 + 0.26 I, -2.40141240226*10^-13} *)

Now do a least-squares fit to get the amplitudes. Recover the trig polynomials. Massage to put into the expected form.

newmat = Map[roots^# &, Range[0, obstime, timescale]/timescale];
Dimensions[newmat]
coeffs = PseudoInverse[newmat].data3
newf = Chop[TrigExpand[ExpToTrig[coeffs.Exp[freqs*t/timescale]]]]
Expand[TrigFactor[newf]]

(* Out[351]= {801, 5} *)

(* Out[352]= {0.552636596429 + 0.233651005386 I, 
 0.552636596429 - 0.233651005386 I, 1.21923674127 - 1.25537315885 I, 
 1.21923674127 + 1.25537315885 I, 
 0.300000000025 + 1.38460601482*10^-18 I} *)

(* Out[353]= 0.300000000025 + 2.43847348254 Cos[2.6 t] + 
 1.10527319286 Cos[5.3 t] - 2.5107463177 Sin[2.6 t] + 
 0.467302010771 Sin[5.3 t] *)

(* Out[354]= 0.300000000025 + 3.49999999956 Sin[2.37079632674 + 2.6 t] + 
 1.20000000005 Sin[1.17079632681 + 5.3 t] *)

This is an essentially perfect recovery of the original signal.

Hi, Hans,

Yes, it may be problematic to run that package on v7, could you provide more details on the error messages? Meanwhile, are you sure you need Fourier analysis for your data? What is the process behind your data?

Perhaps, you might be interested in exploring other methods for data analysis like:

• Wavelet analysis (ref1, ref2), but wavelet staff is available in Mathematica for v8+
• Independent component analysis and it's variations (ref1), there is an implementation of FastICA for Mathematica somewhere on library.wolfram.com and some functions are available in Mathematica, e.g. PrincipalComponents[] and KarhunenLoeveDecomposition[] for v8+.
• Hilbert-Huang transform (ref1), I remember seeing an implementation of the sifting process used in this method somewhere on mathematica.stackexchange.com as well as discrete Hilbert transform.
• Dynamic Mode Decomposition (ref1), I'm not aware of any implementation in Mathematica and will be grateful if someone can provide one

This list is by no means complete, but it can be a good starting point.

P.S. I've modified the code to run (hopefully) on v7, see attached.

I.M.

Attachments:

I would like to make my final remark to this discussion:

But the "process" yielding the data contains exactly two vibrations, so 60-some seems a bit unsatisfactory.

It is important to notice that you do in fact not have just one ore two frequencies in your examples - as their construction suggests. There are discontinuities which lead to a whole variety of other frequencies. This is not some technical limitation of the Fourier method but a deep principle of nature.

And I think Ivan is pretty right when he says: "From my experience the main problem is to determine frequencies as accurate as possible "

Well, this is equivalent to the advice: "Measure as long as you possibly can!".

My general suggestion is: Choose one single good book about Fourier transform and its applications. There a several definitions of Fourier transform, and throughout a single book this definition is consistent - in contrast to sources in the internet! When working with Mathematica make carefully sure that you are using that same definition (by adjusting FourierParameters).

Regards -- Henrik

Hallo Hans,

sorry for my late response! Fourier transformation is a too important topic, so your question should not be left unanswered!

I have tried to construct a hopefully meaningful example of Fourier decomposition, maybe it is helpful. Generally when dealing with Fourier it is always important to know exactly its specific meaning, and this is defined by the option FourierParameters. Have a look under "Details and Options" in the documentation.

Regards -- Henrik

Attachments:

Hallo Hans,

I have tried to do the same using your example.

Regards -- Henrik

Attachments:

Hallo Hans,

Would you think adding zeros at the beginning and end of the data (to force it being the same at the ends of the interval) would help?

I do not think so; by adding zeros you create jumps between your data and the zeros. Maybe you add data which act as interpolation, but this is just a first (and stupid?) guess. An acquisition of a finite amount of data (on -> off) is equivalent to the application of a HeavisidePi "cut out" function. One could use a Gaussian instead. But in fact: This topic is endless! An excellent source is e.g.:

The Fourier transform and its applications /? Ronald N. Bracewell. 3rd ed. Boston : McGraw Hill, c2000.

And what about my (2nd) problem with two functions (see Fourier_MC.nb )

Well, it works exactly the same way! Just plug in the other function...

Regards -- Henrik

Hello Henrik and Ivan.

@Henrik: I fed my 2nd function in your method, and it works. But: using your really nice Manipulate-object it turns out that I need 60-some or so cosine-functions to obtain a good reproduction of the original data (assumed to be results of (error-free) measurements). But the "process" yielding the data contains exactly two vibrations, so 60-some seems a bit unsatisfactory.

Look, in my opinion the process is as follows: make your measurements or take the experimental data, apply "Fourier" and plot the spectrum. Look at the spectrum and get the number of peaks. This should be the number of periodic processes involved. In my example there are two peaks suggesting that there are two vibrations or two periodicities (as it is the case in fact). How to find them?

After all our discussion I got the impression that there is no way but to estimate frequencies (which seems to work quite well), amplitudes and phases (that works by far worse) and do a fit, and I assume yes, following Ivan, at best a linear fit with sines and cosines instead of a nonlinear fit with phases. (but I am not quite sure about this).

And I think Ivan is pretty right when he says: "From my experience the main problem is to determine frequencies as accurate as possible "

@Ivan The revised FREQUENCIES still produces lots of error messages. Of course I will send you all of them. Could you give me your email-address? And if you want to I send you the original problem as well.

Kind regards Hans

Hello Henrik,

here are my final remarks to your final remarks.

First of all thank you again for your remarks and the work you spent on the problem.

But : of course there are two and only two frequencies in the problem. It is constructed exactly in this way. We could think of some physical process known to have two periodicities and described by that "Law". The task is to find the parameters of the "Law", namely frequencies , amplitudes and phases. And indeed, in the "spectrum" we see two peaks. Great.

So, when you say there are lots of frequencies these are artefacts of both, observation and / or computation (finite observation interval, finite observation time and so on). If you will that is in fact a technical limitation and not a deep principle of nature,

So, what I learnt here to get what I want is to plot a spectrum, get the number of peaks, extract the frequencies as exactly as possible and do a Fit for amplitudes and phases, and perhaps the frequencies as well.

Concerning the book about Fourier.... eons ago I read Ziessow "On-line Rechner in der Chemie" with some chapters about FFT, and Lighthill, "Fourier-Analysis", the latter not being really practical. But I feel these books are not of much use for my very practical problem.

Thanks a lot again and Kind regards

Hans