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

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: