FFT code on time series?

15 views (last 30 days)
Sam
Sam on 19 Feb 2013
I have run through the example used here: http://www.mathworks.co.uk/products/matlab/examples.html;jsessionid=20a658a2dcd514f152ada895e6f8?file=/products/demos/shipping/matlab/fftdemo.html and am now trying to do the same with my data. I have a time series that consists of 5127 data points which were collected every 12 hours for just over 7 years. It is of ocean current water movement and I am trying to look at cycles within it.
I have tried to do this and the graph comes up blank. I suspect it is something to do with this part f = 1000/251*(0:127);
I don't know why they have used 1000 or 127.
Would someone be able to explain this? Thanks
Ok sorry, This is what I used. I have a time variable and a measurement of water volume. By following the example above I used this, but have clearly done something wrong:
t=0:.5:2578;
>> plot (Thermocline(1:5157))
>> Y=fft(Thermocline,5157);
>> Pyy=Y.*conj(Y)/5127;
>> f=1000/5127*(0:2563);
plot(f,Pyy(1:2564))
Pyy - 1x5157
Max ans =
1.6713e+06
ans =
4.4727e-07
I am unsure of the range as I don't fully understand what Pyy is. If its time then its half days. If its volume then the range is 18.7 Sv and if its periodicity then I expect to see a strong annual signal (720).
  2 Comments
Oleg Komarov
Oleg Komarov on 19 Feb 2013
Without a code example is hard to understand what you're plotting against what.
Image Analyst
Image Analyst on 19 Feb 2013
Edited: Image Analyst on 19 Feb 2013
Tell us this:
whos Pyy
max(Pyy(:))
min(Pyy(:))
Also what range is Pyy in and is that the range you expected it to be in?

Sign in to comment.

Sign in to answer this question.

Accepted Answer

Youssef Khmou
Youssef Khmou on 20 Feb 2013
hi Sam,
in case you have verified the properties of the signal in the above comment,we can proceed as the following :
We have a signal of length L=3950 taken every 12 hours for 2563.5 days. So formally the sampling rate is Fs= 3950/2563 =~1.54 Hz .
For these types of times series , the range is dynamically changing therefore we plot the frequency in dB to expand the low values and compress the high ones :
1)Solution 1:
Fs=1.55; %
L=length(Thermocline);
N=ceil(log2(L));
FFTherm=fft(Thermocline,2^N);
f=(Fs/2^N)*(0:2^(N-1)-1);
Power=FFTherm.*conj(FFTherm);
figure,
semilogy(f,Power(1:2^(N-1)))
xlabel(' Frequency (Hz)'), ylabel(' Magnitude (w)'),
title(' Power Spectral Density'), grid on;
2) Solution : 2
% Directly you type in the C.P :
psd(Thermocline)
So the frequency is in the range [0,...,0.1] Hz. To see if this result, physically, makes sens, try to check :
"Stochastic climate models, Part I1, Application to sea-surface temperature anomalies and thermocline variability" by CLAUDE FRA NKIG NOUL :
And another paper: "MEAN AND EDDY DYNAMICS OF THE MAIN THERMOCLINE" by GEOFFREY K. VALLIS
http://www.princeton.edu/~gkv/papers/thermocline.pdf
These two papers give a general idea on how much an estimate of the frequency must be .
I hope this helps .
KHMOU Youssef.
  1 Comment
Sam
Sam on 20 Feb 2013
Yes it does, thank you very much for all of your help, I understand it now, thanks for thaking the time for such detailed answers!

Sign in to comment.

More Answers (1)

Youssef Khmou
Youssef Khmou on 19 Feb 2013
hi Sam,
In the demo, they used 1000 and 127 because :
t and y are of length 251, and they computed FFT resulting in same number of points which is 251, so the spectrum is "TWO" sided ( Symmetric) and they created a frequency axis based on a number N satisfying the relation 2^N>=length(y) which is N=8 and this is for adjusting the frequency axis to match the exact frequency .
for a value 1000: its the sampling frequency Fs=1000 Hz , you can anderstand that by looking at the vector t :
t = 0:.001:.25; % t=0:1/Fs......
So to adjust the frequency axis you have to do :
F=Fs/Length(vector)*(0:2^(N-1)-1) ;
Second Part : concerning your data , i think you have to increase the sampling rate to 10 or even 20, to make sure that the Nyquist condition is met.
  9 Comments
Show 7 older comments
Sam
Sam on 19 Feb 2013
I imagine there was a better way to do that but couldn't see one.
Youssef Khmou
Youssef Khmou on 19 Feb 2013
Edited: Youssef Khmou on 20 Feb 2013
hi Sam it is ok, i tried to copy the vector, so to make sure the process of copying was correct verify these properties :
1.Length=3950.
2. Expectation=-17.7831.
3. Variance = 10.6303 .

Sign in to comment.

Categories

Find more on Spectral Measurements in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!