spectral analysis for time series data

27 visualizaciones (últimos 30 días)
Richard
Richard el 28 de Jun. de 2012
I am attempting to analyse a time series with spectral analysis and am working through the example shown in:
Here it states:
Compute the power spectral density, a measurement of the energy at various frequencies, using the complex conjugate (CONJ). Form a frequency axis for the first 127 points and use it to plot the result. (The remainder of the points are symmetric.)
Is there any relevance as to why you form a frequency axis for the first 127 points, why 127? and also in the line
f = 1000/251*(0:127);
where does the 1000 come from?
Sorry for what is probably very basic to most matlab users, but my knowledge of data analysis in the frequency domain is minimal.
From this example I am trying to detect any periodicities in my data, which is composed of hourly measurements recorded for one week (24 * 7 = 168 measurements), I aim to show the diurnal component of the temperature variation. So far I have:
clear all
StartDate = '2011-07-01 00:00';
EndDate = '2011-07-07 23:00';
DateTime=datestr(datenum(StartDate,'yyyy-mm-dd HH:MM'):60/(60*24):...
datenum(EndDate,'yyyy-mm-dd HH:MM'),...
'yyyy-mm-dd HH:MM');
DateTime=cellstr(DateTime);
DecDay = datenum(DateTime)-datenum(2011,0,0);
t = 0:25/length(DecDay):(25-0.1488);
x = sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));
Y = fft(y,length(y));
With regards to the example shown in the link where would I go from here?

Respuesta aceptada

Wayne King
Wayne King el 28 de Jun. de 2012
Editada: Wayne King el 28 de Jun. de 2012
Because you have a real-valued signal, the power spectral density is an even function of frequency. Therefore, there is no need to keep all 251 values in the PSD estimate. For an odd length input (251) if you keep the first round(251/2)+1 you have PSD estimates from 0 frequency (the first value) up to almost the Nyquist frequency.
Personally, I think a better decision would be to have kept floor(251/2)+1
If the length of the input is even (imagine here it was 250), then if you kept from 1 to 250/2+1, you would have PSD estimates from 0 frequency to the Nyquist.
The 1000 is the sampling frequency. The spacing between the DFT bins is Fs/N where Fs is the sampling frequency and N is the length. In this case 1000/251
Again, I would have done:
df = 1000/251; freq = 0:df:500;
So I would have done:
t = 0:.001:.25;
Fs = 1000;
x = sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));
Y = fft(y,251);
Y = Y(1:floor(251/2)+1);
Y = 1/(251*Fs)*(Y.*conj(Y));
df = 1000/251;
freq = 0:df:500;
plot(freq,Y);
Here is an example for your situation:
t = 1:168;
x = cos((2*pi)/12*t)+randn(size(t));
% if you have the signal processing toolbox
[Pxx,F] = periodogram(x,rectwin(length(x)),length(x),1);
plot(F,10*log10(Pxx)); xlabel('Cycles/hour');
ylabel('dB/(Cycles/hour');
% if not
xdft = 1/168*fft(x);
xper = abs(xdft(1:length(x)/2+1)).^2;
df = 1/168;
freq = 0:df:1/2;
plot(freq,10*log10(xper))
xlabel('Cycles/hour');
ylabel('dB/Cycles/hour');
Note that MATLAB uses a convention of additionally scaling one-sided PSD estimates by 2 at frequencies beside 0 and the Fs/2, which I did not do with my fft() example.
  1 comentario
Richard
Richard el 28 de Jun. de 2012
Thank you very much for your detailed answer, great. One final question, how do I convert the xaxis to be hours not cycles/hour. I have tried altering the plot command to plot(1/F... but this messes up the periodogram.

Iniciar sesión para comentar.

Más respuestas (1)

Wayne King
Wayne King el 28 de Jun. de 2012
One way
plot(F,10*log10(xper));
set(gca,'xtick',[1/20 1/12 1/8 1/6 1/4 1/2]);
H = 1./(get(gca,'xtick'));
set(gca,'xticklabel',{H});
xlabel('Period in Hours');

Categorías

Más información sobre Parametric Spectral Estimation en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by