5
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • $\begingroup$ What was the sampling frequency? $\endgroup$
    – Nasser
    Commented Feb 3, 2014 at 10:54
  • $\begingroup$ @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 :/ $\endgroup$
    – holistic
    Commented Feb 3, 2014 at 11:07
  • $\begingroup$ Ok, so fs=0.357 hz? (1/2.8) $\endgroup$
    – Nasser
    Commented Feb 3, 2014 at 11:14
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    Commented Feb 3, 2014 at 11:25
  • $\begingroup$ 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? $\endgroup$
    – holistic
    Commented Feb 3, 2014 at 11:38

1 Answer 1

Reset to default
10
$\begingroup$

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.

Improve this answer
$\endgroup$
16
  • $\begingroup$ Quick question: Is all this method with py and nUniquePts a commonly used thing? $\endgroup$
    – Öskå
    Commented Feb 3, 2014 at 11:40
  • 1
    $\begingroup$ @Ö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/… $\endgroup$
    – Nasser
    Commented Feb 3, 2014 at 11:43
  • $\begingroup$ 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? $\endgroup$
    – holistic
    Commented Feb 3, 2014 at 11:44
  • $\begingroup$ I guess "localhost" has to be replaced by "12000.org" :) $\endgroup$
    – Öskå
    Commented Feb 3, 2014 at 11:45
  • 1
    $\begingroup$ @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). $\endgroup$
    – Nasser
    Commented Feb 3, 2014 at 12:55

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.