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.

I have a dataset of amplitude versus time $(t,A(t))$ and I need to extract the dominant frequency and amplitude, and also get the amplitude at one other specific frequency. My data looks like this:

data = {{-9.75, 13.76}, {-9.5, 14.352}, {-9.25, 15.66}, {-9.,  16.506}, {-8.75, 17.768}, {-8.5, 17.218}, {-8.25, 15.794}, {-8.,  13.18}, {-7.75, 11.58}, {-7.5, 10.524}, {-7.25, 8.428}, {-7., 6.544}, {-6.75, 4.408}, {-6.5, 2.586}, {-6.25,   0.274}, {-6., -2.194}, {-5.75, -4.982}, {-5.5, -6.224}, {-5.25, -9.698}, {-5., -12.22}, {-4.75, -13.986}, {-4.5, -15.372}, {-4.25, -15.1}, {-4., -15.47}, {-3.75, -14.088}, {-3.5, -13.388}, {-3.25, -12.424}, {-3., -12.506}, {-2.75, -11.83}, {-2.5, -9.886}, {-2.25, -8.066}, {-2., -7.434}, {-1.75, -5.23}, {-1.5, -2.418}, {-1.25,   0.252}, {-1., 2.726}, {-0.75, 5.184}, {-0.5, 7.668}, {-0.25,  8.684}, {0., 9.7}, {0.25, 11.866}, {0.5, 13.534}, {0.75, 15.05}, {1., 17.512}, {1.25, 17.99}, {1.5, 16.84}, {1.75,  15.154}, {2., 12.682}, {2.25, 10.358}, {2.5, 9.314}, {2.75, 9.07}, {3., 7.866}, {3.25, 5.244}, {3.5, 2.27}, {3.75, -0.564}, {4., -2.012}, {4.25, -3.078}, {4.5, -5.484},  {4.75, -8.834}, {5., -11.234}, {5.25, -13.162}, {5.5, -15.11}, {5.75,  -16.684}, {6., -16.588}, {6.25, -14.99}, {6.5, -14.336}, {6.75, -12.956}, {7., -12.44}, {7.25, -11.72}, {7.5, -10.854}, {7.75, -8.384}, {8., -6.082}, {8.25, -3.848}, {8.5, -1.772}, {8.75, 0.552}, {9., 3.186}, {9.25, 5.154}, {9.5, 6.976}, {9.75, 8.602}, {10., 9.042}}

And I would like to define a function that takes in account the $x$ scaling for proper scaling of the spectrum. The scaling has been discussed in other questions. But there are some issued that are stopping me.

What I have is (mainly copied from here) :

xyF[d_] := Block[{n, t, y, dt, ft, fy},
  n = Length[d];
  t = d[[All, 1]];
  y = d[[All, 2]];
  dt = Tally[Differences[t]][[1, 1]];
  fy = Abs@Fourier[y];
  ft = RotateRight[Range[-n/2, n/2 - 1]/(n dt), n/2];
  Sort[Transpose[{ft, fy}], (#1[[1]] < #2[[1]]) &]
  ]

My questions are:

  • What is the scaling if the number of points n is odd instead of even?

  • How do I obtain a high-resolution estimation of the dominant frequency?

  • How do I get the intensity of a frequency not explicitly given? (Other than interpolating, given that close to the peak that interpolation may be poor.)

share|improve this question
    
Are these position and heights? The common application of frequency is to sound but you can work with wavenumber and wavelength if your are looking at spatial data. –  Hugh Mar 18 at 15:44
    
@Hugh, I changed the question to time for easier understanding. It is a series of Amplitudes taken in a delay line, time and distance are here equivalent at the fixed speed of light. –  rhermans Mar 18 at 16:14
1  
Possibly related –  bobthechemist Mar 19 at 2:52
    
Probably I should clarify that in the third question I mean at a frequency that is not necessarily the peak of the Fourier-transformed data. –  rhermans Mar 19 at 10:46

2 Answers 2

up vote 11 down vote accepted

I will take the data as a time history. First I assume that the x-values are equally spaced and work out the time increment and frequency increment and then plot the data.

tinc = data[[2, 1]] - data[[1, 1]];
finc = 1/(tinc Length[data]);
ListPlot[data]

Mathematica graphics

This looks like almost two cycles of a sine wave with a frequency of about 0.1 and an amplitude of 17. There is a phase of about an eighth of a cycle. I now take the Fourier transform. As the data is close to a sine wave I think it is best to select FourierParameters to reflect this. We want the amplitude of the sine wave to be readable from the spectra. I choose the parameters to be {-1,-1}. This means that should we be lucky and there is a whole number of cycles in the time interval then the peak in the spectra will be half the amplitude of the sine wave. The half comes from the fact that the spectrum has two peaks (at positive and negative frequencies which are wrapped) and that we want the peak value not the rms. (We get a Sqrt[2] from each of these factors.) Calculate the Fourier transform and the corresponding frequency values.

ft = Fourier[data[[All, 2]], FourierParameters -> {-1, -1}];
ff = Table[f, {f, 0, 1/tinc - finc, finc}];
ListLinePlot[Transpose[{ff, Abs[ft]}], PlotRange -> All]

Mathematica graphics

If we expand the range we can see that we are lucky and the peak appears to coincide with a spectral line. This is not the usual case. Normally the peak lies between two spectral values. Expanding the range...

Mathematica graphics

As we are looking for a sine wave and the number of data points is small I suggest we get most accuracy by fitting in the time domain. The Fourier transform gives us a good starting estimate of the frequency. It is essential to have this otherwise NonlinearModelFit will not work. I also include an unknown for the mean value.

nlmf = NonlinearModelFit[data, 
   a + b Cos[2 \[Pi] f t] + c Sin[2 \[Pi] f t], {a, b, c, { f, 0.1}}, 
   t];

The parameter table shows that the frequency is slightly smaller than 0.1. The standard error on the frequency shows we have about 3 figures of accuracy.

nlmf["ParameterTable"]

Mathematica graphics

We can take the best fit and regenerate the data and over-plot with the original data.

fd = Table[{t, nlmf["BestFit"]}, {t, data[[All, 1]]}];
ListLinePlot[{data, fd}]

Mathematica graphics

If we take away the fitted data we can see if there is anything else to find in the data.

diff = Transpose[{data[[All, 1]], data[[All, 2]] - fd[[All, 2]]}];
ListLinePlot[diff]

Mathematica graphics

It is tempting to see a frequency in this time history but it might be random. The best strategy here is to look again at the Fourier transform.

ftdiff = Fourier[diff[[All, 2]], FourierParameters -> {-1, -1}];
ListLinePlot[Transpose[{ff, Log[10, Abs[ftdiff]]}], 
 PlotRange -> {{0, 1/(2 tinc)}, {-4, 1}}]

Mathematica graphics

I would say this is a flattish spectrum so we have reduced the time history to a white noise process which is the best we can do.

Edit for additional calculations.

The OP wishes to find out what would be the spectral value at frequencies other than those calculated by the original Fourier transform. Let us suppose that the set of frequencies that are required have frequency increment finc1. As an example we will take finc1 = 0.012. Let the number of points in the new spectrum be 1001. Then the sample rate is equal to the frequency increment times the number of points or 12.012 and the time increment is the reciprocal of this i.e. 0.0832501.

finc1 = 0.012;
nn1 = 1001;
tinc1 = 1/(finc1 nn1)

We now calculate an new sampled time history based on the best fit sine wave we have estimated.

data1 = Table[
   a Cos[2 \[Pi] f (n - 1) tinc1] + b Sin[2 \[Pi] f (n - 1) tinc1] /. 
    nlmf["BestFitParameters"], {n, 1, nn1}]; 

If we take the Fourier transform of this we will get the spectrum at the frequencies we require. We also calculate the frequency values.

ft1 = Fourier[data1, FourierParameters -> {-1, -1}];
ff1 = Table[(n - 1) finc1, {n, 1, nn1}];   

We can now plot the spectra with the frequency increment required. I have chosen a value for the frequency increment which deliberately does not get the lucky case we had before where the peak was close to one frequency value. Now we have the typical case where the peak does not correspond to a frequency so there are spectral values for many frequencies. I have expanded the frequency range and chosen the vertical range to be similar to the first spectrum. The actual best fit amplitude is 15.4995 with our spectral scaling being half of this value. We can see how misleading a spectrum can be when the frequency increment does not coincide with the frequency of the time history.

ListPlot[Transpose[{ff1, Abs[ft1]}][[1 ;; 30]], PlotRange -> {All, {0, 8}}]

Mathematica graphics

I hope I have understood what the OP wanted. There may be some confusion over the use of the word amplitude. The best fit sine wave has an amplitude and we then have ordinates of the time history and ordinates of the spectrum.

share|improve this answer
1  
Pleasure to read :) +1 –  Sektor Mar 18 at 22:40
    
This is a fantastic walk through on the general solution, a pleasure to read as @Sektor writes, and I thank you for that. Unfortunately, it does not answer my specific questions: "What is the scaling if the number of points n is odd instead of even?" and "How do I get the intensity of a frequency not explicitly given?" –  rhermans Mar 19 at 10:39
    
@rhermans as to the point regarding frequencies not explicitly given, isn't Hugh arguing that the FT analysis suggests the signal consists of the dominant frequency and white noise? I guess I'm confused on what you mean by 'not explicitly given'. –  bobthechemist Mar 19 at 11:59
    
@bobthechemist, Fitting the time-domain curve will give me a frequency and amplitude, which is part of what I need. It will not give me the amplitude at a nominal frequency very close to that which is not explicitly one of the discrete points for which Fourier gives an explicit value, that is my question. Say the fit tells me the measured frequency is 9.9 Hz, but nominally the frequency was supposed to be 10Hz, and the Fourier discrete transformation, gives me amplitudes at 9.5 Hz and 10.5 Hz. I have half I want, the 9.9 Hz frequency an its amplitude, but I don't have the amplitude at 10 Hz. –  rhermans Mar 19 at 12:13
    
@rhermans I see, so perhaps you are looking to do something like zerofilling the original data to increase the resolution of your frequency domain dataset? –  bobthechemist Mar 19 at 12:19

For a set of data $\{t_i,A(t)\}$ its possible to define a continuous function, properly scaled, that is the Fourier transform of the data set.

1/Sqrt[n] Sum[y[[r]] Exp[(2 \[Pi] I)/n (r - 1) (s n dt)], {r, n}]

Then an implementation for the absolute value of the Fourier transform could be:

cxyFtM[d_, s_] := Block[{n, t, y, dt},
  n = Length[d];
  t = d[[All, 1]];
  y = d[[All, 2]];
  dt = Mean[Differences[t]];
  Abs[1/Sqrt[n] Sum[y[[r]] Exp[(2 \[Pi] I)/n (r - 1) (s n dt)], {r, n}]]
  ]

Let's define some data

ListPlot[testdat = Rest@Table[{t, Sin[1.25 2 \[Pi] t] + Sin[0.25 2 \[Pi] t]}, {t, 0, 10., 1/20}]]

Mathematica graphics

Now we can compare the continuous and discrete results

Show[
 Plot[cxyFtM[testdat, s], {s, 0, 2}, PlotPoints -> 100, 
  PlotRange -> {{0, 2}, {0, 10}}],
 ListPlot[xyF[testdat], PlotStyle -> Red]
 , Epilog -> {Line[{{0.25, 0}, {0.25, 10}}], 
   Line[{{1.25, 0}, {1.25, 10}}]}
 , Frame -> True
 ]

Mathematica graphics

Amplitude is available at any arbitrary frequency, peak position now can be calculated analytically, independently of number of data points.

The estimation of the first peak, a frequency close to zero, is imperfect, but the second peak is bang on. We also see that the discrete estimation of the spectrum between the two main peaks is artificially high.

No intention to diminish the excellent answer by Hugh, just wanted to leave this other approach documented.

share|improve this answer
    
+1 - keep 'em comin :D –  Sektor Nov 28 at 7:31
    
@rhermans This is a nice approach. It is disappointing that the maximum is off for the frequency at 0.25. I suspect that this is due to there being several frequencies. Is the frequency at the maximum only when there is a single frequency in the spectrum? –  Hugh Nov 29 at 22:16
    
@Hugh My guess is that is related to the sampling limit, for frequencies close to half the sampling rate or close to the inverse of the acquisition time. –  rhermans Nov 29 at 22:19

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.