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.