9
$\begingroup$

I'm trying to identify the frequencies in my time history samples, and I can see a frequency in the time history, but can't see it in its Fourier transform. Here it is :

the sample data:

dt = 0.01;(*0.01 second per sample*)
ls={7.18384,9.08503,7.13301,9.03243,7.23692,8.82911,7.48153,8.50053,7.8291,8.09453,8.22514,7.67123,8.60473,7.29656,8.90489,7.02926,9.07228,6.91356,9.07469,6.96968,8.90404,7.19156,8.58122,7.54573,8.15062,7.97689,7.67665,8.41549,7.23096,8.78898,6.88426,9.03305,6.69172,9.10231,6.68698,8.97821,6.87272,8.6724,7.22283,8.22643,7.68225,7.70488,8.17857,7.18711,8.62957,6.75303,8.9597,6.47266,9.10889,6.39174,9.04594,6.52704,8.77139,6.85901,8.32088,7.33824,7.7584,7.88776,7.16829,8.41967,6.64141,8.84366,6.2609,9.0868,6.08919,9.10045,6.15615,8.87392,6.45525,8.43236,6.94108,7.83828,7.53853,7.17843,8.15053,6.55427,8.67686,6.06288,9.02603,5.78403,9.13458,5.76455,8.97352,6.01177,8.55858,6.48966,7.94428,7.12458,7.22101,7.8152,6.49678,8.4485,5.88528,8.91739,5.48272,9.13786,5.35709,9.06375,5.53107,8.69409,5.98259,8.07629,6.64157,7.29757,7.40537,6.47512,8.14882,5.73434,8.74894,5.19349,9.09985,4.93916,9.13459,5.01757,8.83338,5.41929,8.23096,6.08644,7.41062,6.91285,6.49317,7.76769,5.61863,8.50777,4.92325,9.00806,4.51966,9.1747,4.47503,8.96814,4.80297,8.40434,5.455,7.55984,6.33198,6.55581,7.29271,5.54458,8.18165,4.6821,8.847,4.10633,9.17184,3.91187,9.08704,4.13559,8.59066,4.74753,7.74309,5.65489,6.66645,6.71462,5.51928,7.75486,4.47976,8.6024,3.70977,9.10911,3.33653,9.17844,3.42289,8.77952,3.96388,7.95795,4.87716,6.82561,6.02223,5.55073,7.21349,4.32513,8.25592,3.34263,8.96977,2.75835,9.22521,2.67336,8.96055,3.10556,8.19668,3.99626,7.03496,5.20525,5.6427,6.54315,4.23034,7.78851,3.01522,8.73391,2.19205,9.20894,1.89461,9.11841,2.17992,8.45091,3.00859,7.29079,4.25733,5.80068,5.727,4.20335,7.18224,2.74288,8.37806,1.64996,9.10926,1.10137,9.23383,1.19232,8.70843,1.91736,7.5871,3.1692,6.02681,4.75307,4.25335,6.41448,2.53923,7.8804,1.14864,8.89922,0.307694,9.28655,0.154088,8.95041,0.725381,7.91675,1.93765,6.31868,3.60713,4.38985,5.46625,2.41628,7.21338,0.707136,8.55287,-0.470752,9.2482,-0.919297,9.15831,-0.561298,8.26386,0.562301,6.67532,2.27715,4.61492,4.31759,2.39053,6.34975,0.340327,8.03907,-1.21045,9.08947,-2.01384,9.30485,-1.9283,8.61208,-0.958302,7.08603,0.756819,4.93354,2.94639,2.47094,5.26353,0.0706556,7.32281,-1.89254,8.77714,-3.10504,9.358,-3.36416,8.93794,-2.61487,7.53742,-0.963082,5.34406,1.33854,2.66827,3.9233,-0.0836731,6.37173,-2.49189,8.2693,-4.17005,9.28401,-4.84777,9.20861,-4.3986,8.01348,-2.88195,5.83756,-0.523299,2.99222,2.304,-0.10749,5.14496,-2.98081,7.52648,-5.18357,9.04878,-6.32959,9.4304,-6.18952,8.56459,-4.73628,6.56344,-2.13332,3.70414,1.31944,0.401559,5.19872,-2.85418,9.0291,-5.5357,12.3241,-7.17454,14.6869,-7.4465,15.8834,-6.29832,15.9227,-3.97558,15.044,-0.986409,13.6463,2.07958,12.1451,4.72194,10.8576,6.6742,9.93404,7.90659,9.37837,8.56294,9.10124,8.84598,8.99267,8.93815,8.96116,8.95636,8.95544,8.95677,8.95391,8.95455,8.95365,8.95341,8.95371,8.95173,8.95356,8.95088,8.95381,8.9499,8.95321,8.94933,8.95282,8.94887,8.9513,8.94838,8.95018,8.94857,8.94853};

Fourier transform function:

DFT[A_, ht_] := RotateRight[ht/Sqrt[2 \[Pi]]*Fourier[RotateLeft[A, Length[A]/2 - 1], 
              FourierParameters -> {1, 1}], Length[A]/2 - 1];
(*shift the zero frequency to the center*)

plot the time history:

ListPlot[ls, PlotRange -> All, Joined -> True, DataRange -> {0, dt*Length[ls]}, Axes -> False, Frame -> True]

enter image description here

We can see there are a fast frequency with period about 0.02s and a slow frequency with period about 1.5s, but in the Fourier transform we only see the fast frequency(except the zero frequency)

ListPlot[Abs[DFT[ls, dt]]^2, PlotRange -> All, Joined -> True, DataRange -> {- 1/dt/2, 1/dt/2}, Axes -> False, Frame -> True]

enter image description here

So where is the low frequency?

Update

As Simon and bill suggest, the slow oscillation is the beating of two close high frequencies. Since the Fourier transform resolution is 2Pi/(N*dt), where N is the number of sample points, so if I increase the resolution by increasing the number of sample points I should see two separated peaks. So I tried to increase the number of sample points, but I can only see one peak all the time. Here is how I did it:

w = 5.0; dt = 0.66125;
f[x_] := Sin[w x]

ls = Table[f[x], {x, dt, 200 dt, dt}];

we see beating in the plot

ListPlot[ls, PlotRange -> All, Joined -> True, DataRange -> {dt, 200 dt}]

enter image description here

but one on peak in the Fourier transform:

ListPlot[Abs[DFT[ls, dt]]^2, PlotRange -> All, Joined -> True, DataRange -> {-1/dt/2, 1/dt/2}, Axes -> False, Frame -> True]

enter image description here

If we increase the number of sample points, we still only see one peak:

ls2 = Table[f[x], {x, dt, 800 dt, dt}];

ListPlot[Abs[DFT[ls2, dt]]^2, PlotRange -> All, Joined -> True, DataRange -> {-1/dt/2, 1/dt/2}, Axes -> False, Frame -> True]

enter image description here

So where is the problem?

Improve this question
$\endgroup$
7
  • $\begingroup$ If you have v9, try Periodogram[ls, SampleRate -> 100]. The 0 Hz spike has distinct shoulders, but the 0.66 Hz signal cannot be disambiguated. This suggests some windowing is needed, but playing with it, I can't seem to pull it out, either. $\endgroup$
    – rcollyer
    Commented Aug 20, 2013 at 16:34
  • 4
    $\begingroup$ I think this is more of a signal processing question than a Mathematica one. What you see in your data is beating between two high frequencies. You can completely zero the low frequencies and you'll still see it, e.g. ft = Fourier[ls]; ft[[;; 150]] *= 0; ft[[-150 ;;]] *= 0; ListLinePlot[Re[InverseFourier[ft]]] $\endgroup$
    – Simon Woods
    Commented Aug 20, 2013 at 17:04
  • $\begingroup$ @SimonWoods if it is because of beating, I would expect to see two separated peaks if I increase the Fourier transform resolution, but I did't see it. See my updates. $\endgroup$
    – xslittlegrass
    Commented Aug 20, 2013 at 19:30
  • 1
    $\begingroup$ Sorry, I was oversimplifying. The "other" frequency peak is an alias reflected on the other side of the Nyquist frequency 1/(2dt). The beat frequency is twice the difference between your signal frequency and Nyquist fbeat = 2(w/(2Pi)-1/(2dt)) = 0.079 $\endgroup$
    – Simon Woods
    Commented Aug 20, 2013 at 20:12
  • $\begingroup$ @SimonWoods does the "other" frequency really exist in the (improperly sampled) signal or not? Which is correct to say, (1) even we can see a beating in the time history, the "other" frequency may not really exist in the signal, or (2) there are two frequencies exists in the signal and we can see their beating, but we can see only one of them. And could you point me to some references of the Nyquist frequency, especially where you get "an alias reflected on the other side of the Nyquist frequency"? Thanks a lot! $\endgroup$
    – xslittlegrass
    Commented Aug 20, 2013 at 20:39

2 Answers 2

Reset to default
9
$\begingroup$

What is happening in your second example (with the single sine wave giving the "beating") is that you have exceeded the Nyquist frequency: what you are seeing is called aliasing. Here's a simple way to explore this (using your DFT function):

dt = 0.66125; 
f[w_, x_] := Sin[w x];
Manipulate[ls = Table[f[w, x], {x, dt, 200 dt, dt}];
   Column[{ListPlot[ls, PlotRange -> All, Joined -> True, DataRange -> {dt, 200 dt}], 
           ListPlot[Abs[DFT[ls, dt]]^2, PlotRange -> All, Joined -> True, 
           DataRange -> {-1/dt/2, 1/dt/2}, Axes -> False,  Frame -> True]}], {w, 0, 10}]

enter image description here

If you play with the slider, you'll see a single sine wave up to the Nyquist frequency (in both the time and frequency plots). When you get higher than the Nyquist frequency, the signal can do many things, among them show the beating that you noticed. Observe that the frequency line begins to descend as the frequency increases (just after you pass Nyquist). When this has certain relationships to the frequency of the sine wave, you get the beating effect.

So -- in your first question, where we could not see the source of the data (you just gave us a list of values) we had assumed you sampled correctly (i.e., below Nyquist). Simon's guess about beating was a good one. In the second question, you have created a sine wave and are sampling it improperly... there are many things that can go wrong once you do this.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ hi; I think there is a typo: below the Nyquist rate should be above the Nyquist rate $\endgroup$
    – Nasser
    Commented Aug 20, 2013 at 20:10
  • $\begingroup$ @Nasser -- I think I corrected it... I was improperly using "rate" and "frequency" as synonyms. $\endgroup$
    – bill s
    Commented Aug 20, 2013 at 20:16
  • $\begingroup$ @bills I've read the reference you linked, but I'm still not quiet understand. Isn't that the Fourier transform should contains all the frequency information in the signal? How is that mathematically possible to have a beat but only one frequency? Thanks very much for your help. $\endgroup$
    – xslittlegrass
    Commented Aug 20, 2013 at 21:49
  • $\begingroup$ What you are thinking of is the (continuous-time) Fourier transform, which indeed would have only one frequency. Mathematica's Fourier command is the DFT (discrete Foureir transform), and if you break the Nyquist frequency, all sorts of things can happen (as you can see above). Below the Nyquist frequency, the analogy between the (continuous-time) Fourier transform and the DFT is close. Said another way, if you were sampling fast enough, this would not happen. $\endgroup$
    – bill s
    Commented Aug 20, 2013 at 22:47
6
$\begingroup$

Simon has hit the problem on the head (those are not real frequency components you see). But there is a way to discover what these frequencies are, by taking the envelope of the signal and transforming (DFTing) the envelope. 

data = ls - MeanFilter[ls, 5];
decay = 0.06; rise = 0.2;
filt[z_, u_] := Max[decay z + (1 - decay) u, rise z + (1 - rise) u];
env = Drop[FoldList[filt, 0, Abs[data]], 1];

The first line removes the "DC" components from the data and the following lines are a method of extracting the envelope, taken from my answer to this question. (Detail: the Drop is needed because your DFT function doesn't work for odd-length sequences).

enter image description here

Taking the DFT (using your function) of this signal env, and plotting

ListPlot[Abs[DFT[env, dt]]^2, PlotRange -> All, Joined -> True, 
      DataRange -> {-1/dt/2, 1/dt/2}, Axes -> False, Frame -> True]

enter image description here

shows the low frequency part of the signal that you are looking for.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks! That's very helpful. But as you showed the beating frequency is about 6Hz, that means the two frequency should separated apart about 12Hz, but why it is not shown in the Fourier transform? Could you have a look at my updates? $\endgroup$
    – xslittlegrass
    Commented Aug 20, 2013 at 19:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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