2
$\begingroup$

Here is an explanation of fft.fftfreq

I have data with random noise:

T = 1.5; n = 1000; ω = 20 π;
x = Array[Cos[ω # + RandomReal[]] + .3 RandomReal[] &, n, {-T/2, T/2}]

and DirichletWindow:

w = Array[DirichletWindow[4/3 #] + (RandomReal[] - 0.5)/10000 &, n, {-T/2, T/2}]

By multiplying x and w I get:

ListPlot[{xw}, Filling -> None, Joined -> True, DataRange -> {-T/2, T/2}]

enter image description here

From here, I want to LogPlot FFT of xw and w but I don't get the right solution:

XW = Fourier[xw, FourierParameters -> {1, -1}];
W = Fourier[w, FourierParameters -> {1, -1}];
ListLogPlot[{Abs[XW], Abs[W]}, Filling -> None, Joined -> True, PlotRange -> {{-30, 30}, {1, 1000}}]

enter image description here

or:

ListLogPlot[{Abs[XW], Abs[W]}, Filling -> None, Joined -> True, DataRange -> {-30, 30}]

enter image description here

Closer to right solution is first plot, but I can't get symetric plot. In python exist np.fft.fftfreq, which returned float array f contains the frequency bin centers in cycles per unit of the sample spacing. Is there eny simmilar function in Mathematica?

$\endgroup$
  • $\begingroup$ It is not quite clear to me what you are trying to do. Are you looking for the influence of a window or are you just trying to get a correct frequency axis? $\endgroup$ – Hugh Dec 19 '14 at 17:32
  • $\begingroup$ @Hugh I'm trying to get a correct frequency axis. It's simmilar to my first listlogplot of XW but with peaks at x=0, x=10 and x=-10. I also want plotrange on negative x-axis. $\endgroup$ – oarvi Dec 19 '14 at 18:54
  • $\begingroup$ Here is a question and answer that shows how to build the correct frequency vector: mathematica.stackexchange.com/a/33209/1783 $\endgroup$ – bill s Dec 19 '14 at 21:48
  • $\begingroup$ Also think about Periodogram $\endgroup$ – bobthechemist Dec 20 '14 at 15:04
8
$\begingroup$

According to your comment you want peaks at -10, 0 and +10 Hz in the frequency domain. As you want a peak at 0 Hz this means that you need a constant value in the time domain. I will do this without random noise as it will be simpler to check. I will use symbols with lower case first letters as this avoids clashes with Mathematica functions. This is how I generate the time history.

With[{tt = 1.5, n = 1000, f = 10, mean = 0.5},
 x = Table[{t, mean + Cos[2 \[Pi] f t]}, {t, -(tt/2 - tt/n), tt/2, tt/
 n}];
 ]

If you check there are exactly 1000 points. The time increment between points is 1.5/1000. In order to get exactly a whole number of cycles I have to start one increment beyond -tt/2. The time history is given by

ListLinePlot[x, PlotRange -> All, Frame -> True, 
FrameLabel -> {"Time / s", "x"}, PlotTheme -> "Scientific"]

Mathematica graphics

Now we take the Fourier transform. We just need the ordinates from the time domain so we do

y = Fourier[x[[All, 2]], FourierParameters -> {-1, -1}];

You must choose which FourierParameters you like. This is my preference. We also have to make a frequency axis which I think is the main issue from your post. The frequency increment is equal to the sample rate divided by the number of points. The sample rate is the reciprocal of the time increment. Thus

sr = 1000/1.5; (* sample rate *)
finc = sr/1000; (* frequency increment *)
f = Table[f, {f, 0, sr - finc, finc}];  

A standard plot of the spectrum going from zero to the sample rate less one frequency increment would be

s1 = Transpose[{f, Abs[y]}];
ListLinePlot[s1, PlotRange -> All, Frame -> True, 
FrameLabel -> {"Frequency / Hz", "Spectral level"}, 
PlotTheme -> "Scientific"]

Mathematica graphics

However, I note that you would like the spectrum to go from negative to positive values. So we need to wrap the frequency domain appropriately and make a new frequency axis

f2 = Table[f, {f, -(sr/2 - finc), sr/2, finc}];
n1 = Length[s1]/2;
s2 = Transpose[{f2, Abs[Join[y[[n1 + 2 ;; -1]], y[[1 ;; n1 + 1]]]]}];

ListLinePlot[s2, PlotRange -> All, Frame -> True, 
FrameLabel -> {"Frequency / Hz", "Spectral level"}, 
PlotTheme -> "Scientific"]

Mathematica graphics

If we home in around the origin we can see we have the peaks where we wish.

ListLinePlot[s2, PlotRange -> {{-20, 20}, All}, Frame -> True, 
FrameLabel -> {"Frequency / Hz", "Spectral level"}, 
PlotTheme -> "Scientific"]

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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