evaluate the variation in the power of a known periodicity

Consider the following example:

t = 1/24:1/24:365;
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');

This demonstrates the dominant periodicity in the data set. However, if I have a time series for one year worth of data, similar to that shown above, is it possible to look at how a specific frequency changes through time. For example, the time series above shows the hourly variation of a given variable throughout the year. Therefore if we knew that the dominant period was driven by a diurnal cycle (i.e. once per day) then we would know that the dominant periodicity would be 1/24. So, if we know the dominant periodicity to be 1/24, is it possible to see how the power of this specific periodicity changes in time (i.e. throughout the year)?

[Walter Roberson]

Related to http://www.mathworks.co.uk/matlabcentral/answers/45637-analyse-specific-frequencies-in-a-time-series (I'm not sure what the difference is?)

[Star Strider]

The power spectral density should give you the power of the circadian variation as well as that of the annual variation at the appropriate frequencies. The periodogram code you provided, when I used the data I calculated earlier, shows a peak at what in your plot is about 2.4 cycles/hour. One cycle/24 hours is 0.0417 cycles/hour. The predominant annual frequency would be 1.14E-4 cycles/hour.

And thank you, Walter.

[Richard]

I must be wording this incorrectly. What I was hoping of achieving here was something along the lines of using a moving window over the data (over each day, in this case 24 rows), finding the power of the frequency relating to a periodicity of 1/24, and then seeing how much this varies from each window i.e. how much it varies each day. Sorry if this isn't clear.

[Star Strider]

I'm sorry, but I still don't understand what you mean by ‘power of the frequency’. ‘Power’ (or energy) is generally considered to be the square of the amplitude, and this is conserved by the Fourier transform (periodogram or spectrogram) of the squared amplitude ( i.e. the power in a given period remains the same).

You are already (and appropriately with your data) assuming a circadian periodicity. If you want to know the maxima and minima on each day and when they occur, I thought I already provided that ability. You can get the amplitude and time variations from the DayMax and DayMin matrices, then do the requisite maths on them to get the information on the differences between the maxima and minima and the differences in the times they occur, if you want those data.

You can certainly do Fourier transforms on the individual columns of the DayT matrix (or their individual squared amplitudes if you prefer), but it's not obvious to me what information you would get from that. (Then again I'm not a hydrologist or a climatologist, so I'm quite likely missing something.) Those data are already rectangularly-windowed with a 24-hour window, so you shouldn't need to do anything else with them. (The rectangular window effect can produce relatively high-amplitude sidelobes due to the abrupt transitions, so bear that in mind when you interpret your data.)

Answer 1

[Wayne King]

You'll need to use spectrogram() for that. Typically, to get a nice (not blocky looking) output from spectrogram, you overlap your moving windows, but you don't have to do that. You can set the NOVERLAP to 0.

I'm not sure what you intend with your example. It seems like your example has a time vector sampled in increments of 1/24, but you really want your time vector in hours

   t = 1:(365*24);

The problem with just using a 24-hour window when your sampling is once an hour is that you only have 24 data points.

Are you sure you can't use a bit longer window, say 72 hours?

 t = 1:(365*24);
 x = cos((2*pi)/12*t).*(t<72)+1/2*sin((2*pi)/12*t).*(t>72 & t<=96)+1/4*cos((2*pi)/12*t).*(t>96 & t<=144)+randn(size(t));
% compute spectrogram -- moving window of 72 samples, no overlap
[y,f,t,p] = spectrogram(x,72,0,128,1);

The power estimates for the moving windows are contained in the columns of the matrix p

   plot(f,p(:,1))