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.

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

Another approach is to work with a continuous periodogram and try to isolate frequency peaks from that. I'll demonstrate with the running example involving two frequencies.

ff2 = 0.3 + 3.5 Cos[2.6 t + .8] + 1.2 Cos[5.3 t - .4];
obstime = 80;
data2 = Table[ff2, {t, 0, obstime, .1}];
normdata2 = data2 - Mean[data2];

pgram[w_?NumberQ, vals_, gap_] := 
 Abs[Total[vals*Exp[I*w*gap*Range[0, Length[vals] - 1]]]^2]/
  Length[vals]

gap = obstime/Length[normdata2];

Plot[pgram[w, normdata2, gap], {w, 0, 10}, PlotRange -> All]

It is clear that there are a couple of peaks and we can get refinements from NMaximize as below. I remark that it is not hard to find good bounds/starting points just by keeping values using the EvaluationMonitor option to Plot.

{m1, f1} = 
 NMaximize[{pgram[w, normdata2, gap], 2 <= w <= 4}, w]
{m2, f2} = NMaximize[{pgram[w, normdata2, gap], 4 <= w <= 6}, w]

(* Out[1434]= {2453.603870204632, {w -> 2.603084800951954}}

Out[1435]= {301.3987842825805, {w -> 5.30742917643148}} *)

Amplitudes can be estimated using a non-squared variant of the above.

amp[w_?NumberQ, vals_, gap_] := 
 2*Abs[Total[vals*Exp[I*w*gap*Range[0, Length[vals] - 1]]]]/Length[vals]

Offhand I don't know why the factor of 2 is needed in the normalization. And maybe this is not a good method anyway. But ehre is what it gives.

NMaximize[{amp[w, normdata2, gap], 2 <= w <= 4}, w]
NMaximize[{amp[w, normdata2, gap], 4 <= w <= 6}, w]

(* Out[1443]= {3.500386188934522, {w -> 2.60308480095994}}

Out[1444]= {1.226830277890204, {w -> 5.307429138313871}} *)

We even get the frequencies again.

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

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:

Hi, Hans,

From my experience the main problem is to determine frequencies as accurate as possible and to do this usually window functions and peak interpolation is required. Take at a look at the following papers of Hunter and Laskar, where algorithms for accurate frequency identifications are described (as well as for parameters, i.e. amplitudes and phases). I personaly prefer to use Laskar's method for frequencies and the NonlinearModelFit[] for parameters, but LinearModelFit[] or LeastSquares[] can be used as well, if you search for $a_1 \cos(f_1 t) + b_1 \sin(f_1 t) + \dots$ instead of $a_1 \cos(f_1 t + \phi_1) + \dots$.

See also my attached package for signal reconstruction. It's probably not very friendly to read, but you can give it a try. You'll also need to replace CompilationTarget->"C" to CompilationTarget->"WVM" in FREQUENCY.m and remove compiler initialization.

I.M.

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

Thanks a lot. I will have to unhurriedly have a close look at it. But what about the example in my notebook?

Would you mind to give me your email? I could send you the original problem-

Regards Hans

sorry, attached the wrong notebook. take this one

Attachments: