FFT code on time series?
Show older comments
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
on 19 Feb 2013
Without a code example is hard to understand what you're plotting against what.
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?
Accepted Answer
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
These two papers give a general idea on how much an estimate of the frequency must be .
I hope this helps .
KHMOU Youssef.
More Answers (1)
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
Sam
on 19 Feb 2013
Youssef Khmou
on 19 Feb 2013
hi, depends on your data : the vector "Therocline", how that vector was measured ? you do not need to declare t at all if Fs in known in "Thermocline"
Sam
on 19 Feb 2013
Youssef Khmou
on 19 Feb 2013
Edited: Youssef Khmou
on 19 Feb 2013
ok, try this :
% Given a vector Thermocline
Fs=10; % an essay
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,
plot(f,Power(1:2^(N-1)))
xlabel(' Frequency (Hz)'), ylabel(' Magnitude (w)'),
title(' Power Spectral Density'), grid on;
Do you visualize 2 or 3 peaks at least ?
Sam
on 19 Feb 2013
Youssef Khmou
on 19 Feb 2013
ok can you upload the data , i mean the vector ?
Sam
on 19 Feb 2013
Sam
on 19 Feb 2013
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 .
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
