evaluate the variation in the power of a known periodicity

4 views (last 30 days)
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)?
  4 Comments
Star Strider
Star Strider on 13 Aug 2012
Edited: Star Strider on 13 Aug 2012
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.)

Sign in to comment.

Accepted Answer

Wayne King
Wayne King on 13 Aug 2012
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))

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by