Evaluating the variation in the power of a known periodicity

Question

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

Answer 1

Use a moving window.

For example, consider these synthetic data (created in R):

time <- 1:(365*24)/24
amplitude <- time * ((1+length(time))/24-time) / length(time)
set.seed(17)
x <- (x0 <- amplitude * cos(2*pi * (time + 6 * (time/365)^2))) + rnorm(length(time), sd=5)

The dark points are the underlying noiseless signal; the gray line is the signal with added Gaussian noise (subsampled at 7 hour intervals to make it clearer). The signal's amplitude obviously changes, but a little less obviously, the frequency increases with time, too.

Running a window over this series, computing the periodogram for each position, and extracting the value at the dominant frequency found in the first window, produces a series of powers corresponding to the centers of the windows. If there is no change in power over time, a plot of this series should be almost flat. Deviations can be interpreted as temporal changes in power. It is helpful to smooth this plot; a lowess smooth (with a narrow window of 1/4 the total width) is used here.

As a basis for comparison, the same procedure was applied to the underlying signal; it is plotted in gray. The smoothed windowed power function closely approximates the windowed power function for the underlying signal (and is very close at the largest powers).

(The powers shown here are computed using natural logarithms, not as decibels.)

One could make a plot for any desired frequency, in addition to this "dominant" or expected frequency.

Appendix

Here is the R code that produced the figures.

#
# Return the square of the periodogram (normalized by the length of x)
# evaluated at index n.
#
power <- function(x, n) {
  x.hat <- spec.pgram(x, plot=FALSE)
  (x.hat$spec[n] / length(x))^2       #$(TeX bug workaround)
}
#
# Window `power(*,n)` across array `x` using a weighted symmetric 
# window extending to `width` on either side (and therefore of length
# 2*width+1), subsampled every `skip` elements.
#
power.window <- function(x, n, width=1, weight=1, skip=1) {
  i <- seq.int(from=width+1, to=length(x)-width, by=skip)
  sapply(i, function(i) power(x[(i-width):(i+width)] * weight, n))
}
#
# Simulate and plot an interesting time series.
#
time <- 1:(365*24)/24
amplitude <- time * ((1+length(time))/24-time) / length(time)
set.seed(17)
x0 <- amplitude * cos(2*pi * (time + 6 * (time/365)^2))
x <- x0 + rnorm(length(time), sd=5)
i <- seq.int(from=1, to=length(time), by=7)
plot(time[i], x[i], type="l", xlab="Days", ylab="x", main="Data", col="Gray")
points(time, x0, pch=19, cex=0.2)
#
# Find the frequency where power is maximized within the initial window.
#
width <- 30*24
x.hat <- spec.pgram(x[1:(2*width+1)], plot=FALSE)
i.max <- which.max(x.hat$spec) #$ TeX bug workaround
#
# Compute and plot the power moving window for the underlying 
# series (`x0`) and its noisy version (`x`).  For plotting, subsample
# the series.
#
skip <- 7
x.power <- power.window(x, n=i.max, width=width, weight=1, skip=skip)
x0.power <- power.window(x0, n=i.max, width=width, weight=1, skip=skip)
x.power.smooth <- lowess(x.power, f=1/4)
i <- seq(from=(width+1)/24, by=skip/24, length.out=length(x.power))
main <- sprintf("Moving power at 1/%4.1f hours, %d hour window", 
                1/x.hat$freq[i.max], 2*width+1) #$ TeX bug
plot(i, log(x.power), type="l", lwd=2, xlab="Time (days)", main=main)
lines(i, log(x.power.smooth$y), lty=2, lwd=3, col="Blue")
lines(i, log(x0.power), lwd=2, col="Gray")
legend(30, -4, "Underlying", bty="n", box.col="White", box.lwd=0, lwd=2, col="Gray")
legend(30, -5, "Windowed", bty="n", box.col="White", box.lwd=0, lwd=2)
legend(30, -6, "Smooth", bty="n", box.col="White", box.lwd=0, lwd=3, col="Blue", lty=2)