Filter löschen
Filter löschen

FFT code on time series?

33 Ansichten (letzte 30 Tage)
Sam
Sam am 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 Kommentare
Oleg Komarov
Oleg Komarov am 19 Feb. 2013
Without a code example is hard to understand what you're plotting against what.
Image Analyst
Image Analyst am 19 Feb. 2013
Bearbeitet: Image Analyst am 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?

Melden Sie sich an, um zu kommentieren.

Melden Sie sich an, um diese Frage zu beantworten.

Akzeptierte Antwort

Youssef Khmou
Youssef Khmou am 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 Kommentar
Sam
Sam am 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!

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Youssef Khmou
Youssef Khmou am 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 Kommentare
7 ältere Kommentare anzeigen
Sam
Sam am 19 Feb. 2013
I imagine there was a better way to do that but couldn't see one.
Youssef Khmou
Youssef Khmou am 19 Feb. 2013
Bearbeitet: Youssef Khmou am 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 .

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Spectral Measurements finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by