Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]

Mathematica graphics

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.

share|improve this question
1  
This is not a Mathematica question, it is an algorithmic question. There are many ways of approximating the frequencies of a signal more accurately: interpolation in the frequency domain (basically fitting a polynomial to the values surrounding the peaks), taking multiple FFTS (sometimes called a "phase-vocoder" in audio), zero-padding the input signal. This is an important and well-studied problem. –  bill s Feb 23 at 15:47
    
You'll find the answer here: stats.stackexchange.com/questions/16230/… –  Szabolcs Feb 23 at 15:51
    
Look up sinc interpolation or Goertzel's algorithm... also see this answer: Obtain a signal's peak value if it's frequency lies between two bin centers I agree with others that this is not a Mathematica question, but a DSP one. –  rm -rf Feb 23 at 16:03
    
Another possible duplicate: mathematica.stackexchange.com/q/38217/7167 –  bobthechemist Feb 23 at 16:13
    
Maybe a weighted average from those positions? aft = Abs[ft]/Max[Abs[ft]]; lower = {9, 10}; v1 = (lower - 1).aft[[lower]]/Total[aft[[lower]]] 2*v1 Out[139]= 8.61412 Out[140]= 17.2282 –  Daniel Lichtblau Feb 23 at 19:58

5 Answers 5

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]

enter image description here

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.

share|improve this answer
    
Thanks for this. I tried it out on my real data. I had 4001 points which had some noise and a small amplitude noise frequency away from my frequency of interest. Your method only gave four candidates of which the best was 31.6456. My method with the symbolic formula for the Fourier transform gave 32.78996 with a standard error of 0.002. This was the same data that I use to test Szabolcs method (in comments after original question). I could not run the second part of your method because I did not have binpos. Is it missing from your code? –  Hugh Feb 27 at 23:35
    
Sorry about that. I must have messed up a bit in cut/paste mode. I put in the errant line. My guess is that this second part will be more robust than the first. –  Daniel Lichtblau Feb 27 at 23:52
    
If you could make your test data available in some way (that is easily imported into Mathematica, which presumably is already the case since you use it), that would be useful. –  Daniel Lichtblau Feb 27 at 23:59

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 ---

share|improve this answer
    
@ DanielLichtblau Many thanks for this. The speed up will be very very helpful. Minutes instead of hours. I think I may have found another method which I will add. –  Hugh Feb 25 at 21:49
    
decent answer. as usual. thank you for that. –  Stefan Feb 25 at 23:32

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}]
share|improve this answer
1  
It's fine to answer your own question. Most especially when the problem is nontrivial and others have not done so. –  Daniel Lichtblau Feb 24 at 16:07
    
@Daniel Lichtblau Thank you. What I need is help to make my method go faster. –  Hugh Feb 24 at 17:16
    
By the way, I think your code above is missing a line ft = Fourier[th, FourierParameters -> {-1, -1}];. It is in the original post but really the response should have full stand-alone code unless indicated otherwise. –  Daniel Lichtblau Feb 24 at 20:08
    
@DanielLichtblau Woops! Thanks Daniel I have fixed it now. –  Hugh Feb 25 at 21:42

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.

share|improve this answer
    
@bill s who kindly responded to the original post suggested that this was not a Mathematica question but an algorithm question. It appears that Mathematica's ability to find symbolic sums has made this problem straightforward. Thanks to everyone who helped. –  Hugh Feb 25 at 23:11

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]

enter image description here

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]

enter image description here

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).

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.