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:

http://www.mathworks.co.uk/products/matlab/examples.html;jsessionid=20a658a2dcd514f152ada895e6f8?file=/products/demos/shipping/matlab/fftdemo.html

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?

Answer 1

[Wayne King]

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.

[Richard]

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.

Answer 2

[Wayne King]

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');

[Richard]

Great