MATLAB: Evaluate the variation in the power of a known periodicity

signal processing

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)?

Best Answer

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))
Related Question