spectral analysis for time series data
显示 更早的评论
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?
采纳的回答
Wayne King
2012-6-28
编辑：Wayne King
2012-6-28
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 个）
Wayne King
2012-6-28
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');
类别
在 帮助中心 和 File Exchange 中查找有关 Spectral Estimation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
