How do you find the frequency and amplitude from Fourier?

Question

Edit

This is a popular question with several answers. This edit organises the answers and gives links.

1. Fourier just gives y values (ordinates) if you wish to read off frequencies and amplitudes you need to add an x axis (abscissa) and scale the ordinates. Details for doing this and other basic information on Fourier are given here.

There are other methods for finding the frequency and amplitude besides Fourier these are:

2. A method due to Daniel Lichtblau that looks for periodicities in the time history answer below.

3. A method using Prony series. Code from Daniel Lichtblau answer below.

4. Straightforward fitting of sine wave to data. There are pitfalls in formulating this problem which are discussed here.

5. A method using auto-correlation due to Szabolics. Answer below.

The original question was how to use the output from Fourier when the computed values do not lie at a single point. This is answered in two ways below.

Original Question

The problem is that a typical harmonic time history will not have an integer number of cycles in the data. For example, with the following time history you get a numerical Fourier transform.

sr = 400; (* sample rate *)
nn = 200;  (* number of points *)
f = 17.2;   (* frequency in Hz *)
A = 1; B = 0.1; (* Coefficients of Cos and Sin *)
th = N @ Table[A Cos[2 Pi f t] + B Sin[2 Pi f t], {t, 0, (nn - 1)/sr, 1/sr}];
ft = Fourier[th, FourierParameters -> {-1, -1}];
n = Position[#, mx = Max[#]] & [Abs[ft]][[1, 1]];
ListPlot[Abs[ft[[1 ;; n + 10]]], PlotRange -> All]

This produces a plot with several points making up the peak. The frequency resolution is 2 Hz and there is no point at the harmonic frequency of 17.2 Hz. Clearly the frequency is bracketed between the 9 th and 10 th point (16 Hz and 18 Hz). How do we find the exact frequency and amplitude? I have considered extending by zeros but although this increases the resolution it still results in points on either side of the frequency. Also, I have considered calculating points with non-integer values. However, the resulting spectrum does not have a maximum at the frequency of the time history. My actual data includes noise and other harmonics hence my use of Fourier which is good at concentrating the data I am interested in around a few points in the spectrum. I have worked out one very poor way of finding the frequency and amplitude and will post this when I have worked out how to use StackExchange more fully.

Answer 1

I also will give a different method, based heavily on one from this prior period estimation post. Before getting to the somewhat thorny code let me list its advantages and disadvantages.

On the plus side:

(1) It does not require knowing the form of the sinusoid. It need not be a sinusoid.

(2) It is fairly fast.

(3) It does not require equal spacing of x values. It sometimes does better if they are random.

But alas:

(4) It is not always perfectly accurate. But the example below seems promising.

So here we go.

nn = 1000;
sr = 2*nn;
f = 17.743;
A = 1; B = 0.1;
SeedRandom[1111];
xvals = Prepend[Accumulate[RandomReal[2/sr, nn - 1]], 0.];
(*xvals=Range[0.,(nn-1)/sr,1./sr]; for even spacing, e.g. if an FT is desired*)
th2 = N@Table[A Cos[2 Pi f t] + B Sin[2 Pi f t], {t, xvals}];
th3 = Transpose[{xvals, th2}];
ListPlot[th3]

We now sort differences of consecutive y values, get corresponding differences of x values, and discard those that are too close to be a period apart.

sortd = th3[[Ordering[th3[[All, 2]]]]];
xdiffs = Differences[sortd[[All, 1]]];
ydiffs = Abs[Differences[sortd[[All, 2]]]];
ynorms = Map[Norm, Partition[sortd[[All, 2]], 2, 1]];
weighteddiffs = ydiffs/ynorms;
sortdxdiffs = Abs[xdiffs[[Ordering[weighteddiffs]]]];
gaps = Select[sortdxdiffs[[1 ;; 10]], # >= .03 &]

(* Out[1725]= {0.394451, 0.224346, 0.11375, 0.281804, 0.169081, \
0.225444, 0.223596, 0.0563591, 0.0563603, 0.11272} *)

This gives good contenders for separation by an integer multiple of periods lengths. And the frequency should be the reciprocal of the period.

candidates = 1/gaps

(* Out[1747]= {2.53517, 4.45741, 8.79121, 3.54857, 5.91434, 4.43569, \
4.47236, 17.7434, 17.743, 8.87154} *)

Also we know, from e.g. a Fourier transform based on evenly spaced x values, that the period is around 17. So take those values that are in that ballpark.

best = Select[candidates, 16 <= # <= 18 &];
Mean[best]

(* Out[1765]= 17.7432 *)

Here is a perhaps safer way, more in line with what is at the referenced link. We start just prior to where we computed gaps above. This time we bin the x differences and use the largest bin. The idea is that those should correspond to a full number of periods (or at worst half periods), whereas other close differences could have arisen by accident of a y value on an "upward" arm could have been very close to one on a "downward" arm, and so the corresponding x values separation would not be meaningful.

bins = BinLists[sortdxdiffs, {.01, .4, .002}];
cts = Map[Length, bins];
binpos = Position[cts, Max[cts]][[1, 1]];
bin = bins[[binpos]];
Length[bin]

Now remove the smallest and largest as they could have incidental "outliers".

sbin = Sort[bin][[20 ;; -20]];
bvals = {Length[sbin], Min[sbin], Max[sbin], Mean[sbin]}

(* Out[1769]= 104

Out[1771]= {66, 0.0563391, 0.0563822, 0.0563595} *)

Finally estimate the frequency, or a "nice" fraction thereof, as before. Since we expect it to be around 17 we can always multiply by a suitable integer to put it in that range, at least if we went about this correctly. But the result below indicates that we did, and we do not need to make such an adjustment.

1/bvals[[-1]]

Out[1772]= 17.7432

This is to be expected, insofar as we expect more "close" y values from separations by one period than by multiples of one period.

Answer 2

A similar question was recently raised on Wolfram Community.

http://community.wolfram.com/groups/-/m/t/528474?p_p_auth=33hR9Jgo

One of my responses there showed how one might code up a (reasonably) numerically stable version of Prony's method to recover frequencies. I'll do similar with this example.

Repeat the setup:

sr = 400;(*sample rate*)
nn = 200;(*number of points*)
f = 17.2;(*frequency in Hz*)
A = 1; B = \
0.1;(*Coefficients of Cos and Sin*)
th = N@Table[
    A Cos[2 Pi f t] + B Sin[2 Pi f t], {t, 0, (nn - 1)/sr, 1/sr}];

We do not need the FT in this case. We proceed to find out how many frequencies are needed, then on to recover them.

mat = Most[Partition[th, Floor[nn/2], 1]];
keep = Length[SingularValueList[mat]]
mat2 = Most[Partition[th, keep, 1]];
rhs = Drop[th, keep];
soln = PseudoInverse[mat2].rhs
roots = xx /. NSolve[xx^keep - soln.xx^Range[0, keep - 1] == 0, xx]
Map[Norm, roots]
freqs = Log[roots]

(* Out[466]= 2

Out[469]= {-1., 1.92744735658}

Out[470]= {0.96372367829 - 0.26690198932 I, 
 0.96372367829 + 0.26690198932 I}

Out[471]= {1., 1.}

Out[472]= {-9.08995101412*10^-16 - 
  0.270176968209 I, -9.08995101412*10^-16 + 0.270176968209 I} *)

Now recover the (complex) exponentials with amplitudes. We simply absorb the phase shifts into the amplitudes.

newmat = Map[roots^# &, Range[0, nn - 1]];
coeffs = PseudoInverse[newmat].th
newf = Chop[TrigExpand[ExpToTrig[coeffs.Exp[freqs*t*sr]]]]

(* Out[493]= {0.5 + 0.0500000000001 I, 0.5 - 0.0500000000001 I}

Out[494]= 0. + 1. Cos[108.070787283 t] + 0.1 Sin[108.070787283 t] *)

Amplitudes obviously match. Also notice that we did in fact recover the frequency.

Chop[freqs*sr/(2*Pi)]

(* Out[489]= {0. - 17.2 I, 0. + 17.2 I} *)

Anyway, I thought it would be useful to have some form of Prony implementation code on MSE.

Answer 3

There are probably better ways to go about this but one possibility is to use memoization to avoid repeated Fourier transforms. Also you only need a few elements from the FT. The code below will compute and cache the values needed. It can be memory intensive though so if you have many such computations to do it would be a good idea to clear the definitions every so often.

Clear[ff, expons, tabl];
tabl[f_, a_, b_, sr_, nn_] := 
 tabl[f, a, b, sr, nn] = 
  Table[a Cos[2. Pi f t] + b Sin[2. Pi f t], {t, 0., (nn - 1.)/sr, 
    1./sr}]
expons[nn_, n_] := 
 expons[nn, n] = Exp[-2.*Pi*I*(n - 1)/nn*Range[0., nn - 1.]]
ff[f_, a_, b_, sr_, nn_, n_] := 
 ff[f, a, b, sr, nn, n] = 1/nn*expons[nn, n].tabl[f, a, b, sr, nn]

modelH2[sr_, nn_][{f_?NumericQ, A_?NumericQ, B_?NumericQ}][n_, ri_] :=
  Module[{ft, op},
  ft = ff[f, A, B, sr, nn, n];
  Which[ri == 0, op = Re[ft], ri == 1, op = Im[ft]]; op]

Now you can run your fit but as modified below, and everything else as in your posts.

FindFit[data, modelH2[sr, nn][{fr, a, b}][x, y], 
  {{fr, 2*(n - 1)}, {a, 1.2}, {b, 1.2}}, {x, y}]

It uses half the number of function evaluations and most are recalled from prior ones, so it's around 20x faster.

Caveat: Most frequencies above around 7 seem fine (smaller than that and data cannot grab so large a range around the peak). But 23.1 causes both methods to give a negative value for b and a frequency of 376.9. It seems that some rare cases will need to be subtracted from sr.

--- edit ---

Right...of course...there should be two peaks, and if the one on the right is a hair larger than the one on the left, then we get a seemingly "large" frequency. So it might be good to place this line after the one that calculates the position n of the max peak.

If[n > nn/2, n = nn - n + 1]

--- end edit ---

Answer 4

I am answering my own question (is this good form?). It is a very poor, brute force, method for finding the frequency and amplitude. It is based on using FindFit to generate a match between the time history and the spectrum. One problem is that FindFit cannot handle complex data so it is necessary to split the spectral data into two vectors one for real data and one for imaginary data. I do this using the function GetData. The model is a disaster for speed because I have to do a full Fourier transform and then extract the one value I need. This is done by modelH. I have tried to calculate an individual Fourier coefficient by just multiplying by the appropriate cos and sine terms but this turned out to take longer than taking the full Fourier transform.

    GetData[d_, n1_, n2_] := 
    Flatten[Table[{{n, 0, Re[d[[n]]]}, {n, 1, Im[d[[n]]]}}, {n, n1, n2, 1}], 1];

    modelH[sr_, nn_][{f_?NumericQ, A_?NumericQ, B_?NumericQ}][n_, ri_] := 
    Module[{ft, op},
    ft = Fourier[
     Table[A Cos[2 Pi f t] + B Sin[2 Pi f t], {t, 0, (nn - 1)/sr, 1/sr}], 
     FourierParameters -> {-1, -1}];
     Which[
     ri == 0, op = Re[ft[[n]]], ri == 1, op = Im[ft[[n]]]]; op]

Making the data in the question again and then using FindFit gets the correct answer. A further problem is that my model does not work when supplied to NonLinearModelFit so I can't get any statistics on the quality of fit.

    sr = 400;
    nn = 200;
    f = 17.2;
    A = 1; B = 0.1;
    th = N@Table[ 
       A Cos[2 Pi f t] + B Sin[2 Pi f t], {t, 0, (nn - 1)/sr, 1/sr}];
    ft = Fourier[th, FourierParameters -> {-1, -1}];
     n = Position[#, mx = Max[#]] & [Abs[ft]][[1, 1]];
    data = GetData[ft, n - 4, n + 4];
    FindFit[data,  modelH[sr, nn][{fr, a, b}][x, y], {{fr, 17.54}, {a, 1.2}, {b, 1.2}},  {x, y}]

Answer 5

You can try a naive method to estimate the period: compute the auto-correlation (which will be much smoother than the original data) and measure the distances between its peaks. We can't easily find peaks in noisy data, but we can in smooth data.

---

Here's an example. Let's make some noisy weird-shaped sample data:

trafo = Function[x, 1/2 x^2 + x^3 + x];
data = Table[trafo@Sin[x] + RandomReal[.1 {-1, 1}], {x, 0, 100, 0.05}];
ListPlot[data]

Notice that the period is 2Pi/0.05 = 125.664.

Let's compute the correlation of the full dataset with half of it. (Yes, this is a rather arbitrary choice)

corr = ListCorrelate[Take[data, Quotient[Length[data], 2]], data];
ListPlot[corr]

Notice that while the data is noisy, the auto-correlation is rather smooth, so it's easier to find maxima. Let's then measure the distance of consecutive maxima (throwing away the first an last result, as those are distances from the edge of the dataset).

gaps = Length /@ Rest@Most@Split[UnitStep@Differences[corr], #1 <= #2 &]

(* ==> {125, 126, 126, 125, 126, 126} *)

We get very similar distances, which reassures us that nothing went wrong so far.

Now let's compute the average of these integers for a more accurate estimate:

Mean[N[gaps]]

(* ==> 125.667 *)

This is very accurate (the true period is 125.664).

This method will work if: 1. you have at least 6 periods in the dataset 2. the period is very steady and not changing within the dataset 3. there isn't too much noise (though it's fairly tolerant).

Answer 6

The approach is to get a symbolic expression for the numerical Fourier transform. Help gives the series that Fourier sums. Mathematica has the ability to sum certain series to give an algebraic solution. It can sum this one. I was amazed! The important point is that this is an exact result and includes the minor effects that occur when numerical Fourier transforms are taken such as slight aliasing when the number of points per cycle is close to 2. The sum works out of the box as follows:

ClearAll[A, B, nn, r, s, \[Beta], n1, \[Alpha]];
s1 = 1/nn Sum[(A Cos[2 \[Pi] \[Beta]  (r - 1)/nn] + 
   B Sin[2 \[Pi] \[Beta]  (r - 1)/nn] ) E^(-2 \[Pi] I (r - 1) (s - 1)/nn) , {r, 1, nn}];
s2 = FullSimplify[s1, {{nn, s} \[Element] Integers, 0 < s < nn}];
rs = FullSimplify[ComplexExpand[Re[s2]]];
is = FullSimplify[ComplexExpand[Im[s2]]];

In the above A and B are the amplitudes of the in-phase and out-of-phase components and \[Beta] is the frequency at a point number which will not be an integer. The point number in time is r and the point number in frequency is s. The number of points is nn. If you look at the results they are long and not particularly informative. The real and imaginary parts may be formed into functions.

ClearAll[rFT, iFT];
rFT[A_, B_, \[Beta]_, nn_, s_] := Evaluate[rs];
iFT[A_, B_, \[Beta]_, nn_, s_] := Evaluate[is];

Here is a test example.

sr = 400;
nn = 200;
f = 17.2;
Ac = 2; Bs = 1.; 
th = N@Table[ 
Ac Cos[2 \[Pi] f t] + Bs Sin[2 \[Pi] f t], {t, 0, (nn - 1)/sr, 1/sr}];
ft = Fourier[th, FourierParameters -> {-1, -1}];
n = Position[#, mx = Max[#]] & [Abs[ft]][[1, 1]];
ListPlot[Abs[ft[[1 ;; n + 10]]], PlotRange -> All]

I take data from around the peak to perform the fitting. I will take 5 points but it is unclear how many to take. It is possible to fit the absolute value of the data or to fit the complex data. Fitting the complex data is better because it uses more of the available information. One problem with Mathematica's FindFit is that it does not handle complex data so it is necessary to construct a two-dimensional data set from the real and imaginary parts. Give the two-dimensional data {x,y} values where x is the point number and y = 0 for real and y = 1 for imaginary. I form the data using formData.

formData[d_, n1_, n2_] := 
  Flatten[Table[{{n, 0, Re[d[[n]]]}, {n, 1, Im[d[[n]]]}}, {n, n1, n2, 1}], 1];
data = formData[ft, n - 2, n + 2];
FindFit[data, (1 - y) rFT[A, B, fr nn/sr, nn, x] + 
  y iFT[A, B, fr nn/sr, nn, x], {{A, 3}, {B, 2}, {fr, 17}}, {x, y}]

{A -> 2., B -> 1., fr -> 17.2}

The method works when used with NonlinearModelFit which enables statistics to be identified.