# Extracting frequencies with Fourier transformation

I know this question came up so many times, but after browsing through many of the threads, I'm still not sure, if I'm doing it correctly with my dataset, so I'll give it a try.

Now I want to look for possible dominant frequencies and it tried it like this:

ListLinePlot[Log[10, Abs[Fourier[data]]], PlotRange -> Automatic]


and I get this:

Correct me if I'm wrong, but I don't see any dominant frequencies in here. Do you guys come to the same conclusion? Honestly, I think I'm doing it all wrong because I'm really not sure which of the many functions of mathematica to use.

EDIT: Now I'm totally confused. I just looked again in some of my books and what they did there was drawing the frequencies on the x-Axes and the Power on the y-Axes. Then they normalised it to the total power. How to do this?

• What was the sampling frequency? Feb 3, 2014 at 10:54
• @Nasser: It's data from an fMRI scan. So a scan was done every 2800 ms. Sorry, I'm just trying to get more into this stuff :/ Feb 3, 2014 at 11:07
• Ok, so fs=0.357 hz? (1/2.8) Feb 3, 2014 at 11:14
• btw, this looks like a very low sample rate. You can adjust the code below if you think it is different. 0.357 samples per second is very low but you can check it on your end. If you meant it was 2800 samples per second, then change fs below to 2800. Feb 3, 2014 at 11:25
• I'm not really sure, but every 2,8 sec a scan was done, i.e. data was collected..so the sampling rate must be 2,8 sec? Feb 3, 2014 at 11:38

Assuming sample rate is 2.8 second per sample, then may be this:

SetDirectory[NotebookDirectory[]];
data = Import["testdata.txt", "List"];
ListPlot[data, Joined -> True]


py = Fourier[data, FourierParameters -> {1, -1}];
nSamples = Length[data];
nUniquePts = Ceiling[(nSamples + 1)/2];
py = py[[1 ;; nUniquePts]];
py = Abs[py];
py[[1]] = 0;  (*zero out the DC frequency, for now, outliner, messes up the scale *)
py = py/nSamples;
py = py^2;
If[OddQ[nSamples], py[[2 ;; -1]] = 2*py[[2 ;; -1]], py[[2 ;; -2]] = 2*py[[2 ;; -2]]];
fs = 0.357;  (*2.8 second per sample*)
f = (Range[0, nUniquePts - 1] fs)/nSamples;
ListPlot[Transpose[{f, py}], Joined -> True,
FrameLabel -> {{"|H(f)|", None}, {"hz", "Magnitude spectrum"}},
ImageSize -> 400, Frame -> True, RotateLabel -> False,
GridLines -> Automatic, GridLinesStyle -> Dashed, PlotRange -> All]


The above is the magnitude spectrum. So you can see now where the frequencies with most energy are.

• Quick question: Is all this method with py and nUniquePts a commonly used thing?
– Öskå
Feb 3, 2014 at 11:40
• @Öskå I do not know if it very commonly used t, but I use this method all the time. In Matlab some use this method. localhost/my_notes/mma_matlab_control/KERNEL/… Feb 3, 2014 at 11:43
• Thank you @Nasser! Why is it that you did not use Fourier or some similar function? I mean..where is the fouriertransform in all this? Feb 3, 2014 at 11:44
• I guess "localhost" has to be replaced by "12000.org" :)
– Öskå
Feb 3, 2014 at 11:45
• @holistic if you mean the py/nSamples then this is something commonly done, to average the power. Its been a while since I did dsp, so I do not remember the details now. But the above is a common way to do it (scale the psd by number of samples). Feb 3, 2014 at 12:55