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

If have some sample data I know not much about. You can download it here

enter image description here

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:

enter image description here

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?

share|improve this question
    
What was the sampling frequency? –  Nasser Feb 3 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 :/ –  holistic Feb 3 at 11:07
    
Ok, so fs=0.357 hz? (1/2.8) –  Nasser Feb 3 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. –  Nasser Feb 3 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? –  holistic Feb 3 at 11:38
show 1 more comment

1 Answer 1

up vote 10 down vote accepted

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

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

Mathematica graphics

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]

Mathematica graphics

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

share|improve this answer
    
Quick question: Is all this method with py and nUniquePts a commonly used thing? –  Öskå Feb 3 at 11:40
1  
@Ö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/… –  Nasser Feb 3 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? –  holistic Feb 3 at 11:44
    
I guess "localhost" has to be replaced by "12000.org" :) –  Öskå Feb 3 at 11:45
1  
@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). –  Nasser Feb 3 at 12:55
show 11 more comments

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.