I have a bunch of time series whose power spectra (FFT via

**R**'s

**spectrum()** function) I've been trying to visualize in an intuitive, aesthetically appealing way. At first, I just used lattice's

**bwplot**, but the spacing of the X-axis here really matters. The spectra's frequencies aren't regularly-spaced categories, which is the default of

**bwplot**. If all the series are of the same length, then the coefficients of the are all estimated at the same frequencies, and one can use

**plyr**'s split-apply-combine magic to compute the distribution of coefficients at each frequency.

I spent some time trying to faithfully reproduce the plotting layout of the figures produced by, e.g.

**spectrum(rnorm(1e3))** with respect to the axes. This was way more annoying than i expected...

To include a generic example, I started looking around for easily generated, interesting time series. The

logistic map is

so damned easy to compute and generally pretty interesting, capable of producing a range of frequencies. Still, I was rather surprised by the final output. Playing with these parameters produces a range of serious weirdness -- apparent ghosting, interesting asymmetries, etc..

Note that I've just used the logistic map for an interesting, self-contained example. You can use the exact same code to generate your own

**spectrum()** box-and-whisker plot by substituting your matrix of time series for

**mytimeseries** below. Do note here that series are by row, which is not exactly standard. For series by column, just change the margin from 1 to 2 in the call to

**adply** below.

## box-and-whisker plot of FFT power spectrum by frequency/period
## vary these parameters -- r as per logistic map
## see http://en.wikipedia.org/wiki/Logistic_map for details
logmap = function(n, r, x0=0.5) {
ret = rep(x0,n)
for (ii in 2:n) {
ret[ii] = r*ret[ii-1]*(1-ret[ii-1])
}
return(ret)
}
## modify these for interesting behavior differences
rfrom = 3.4
rto = 3.7
rsteps=200
seqlen=1e4
require(plyr)
require(ggplot2)
mytimeseries = aaply(seq(from=rfrom, to=rto, length.out=rsteps), 1,
function(x) {
logmap(seqlen, x)
})
## you can plug in any array here for mytimeseries
## each row is a timeseries
## for series by column, change the margin from 1 to 2, below
logspec = adply( mytimeseries, 1, function(x) {
## change spec.pgram parameters as per goals
ret = spectrum(x, taper=0, demean=T, pad=0, fast=T, plot=F)
return( data.frame(freq=ret$freq, spec=ret$spec))
})
## boxplot stats at each unique frequency
logspec.bw = ddply(logspec, 'freq', function(x) {
ret = boxplot.stats(log10(x$spec));
## sometimes $out is empty, bind NA to keep dimensions correct
return(data.frame(t(10^ret$stats), out=c(NA,10^ret$out)))
})
## plot -- outliers are dots
## be careful with spacing of hard returns -- ggplot has odd rules
## as close to default "spectrum" plot as possible
ggplot(logspec.bw, aes(x=1/(freq))) + geom_ribbon(aes(ymin=X1,
ymax=X5), alpha=0.35, fill='green') +
geom_ribbon(aes(ymin=X2, ymax=X4), alpha=0.5, fill='blue') +
geom_line(aes(y=X3)) +
geom_point(aes(y=out), color='darkred') +
scale_x_continuous( trans=c('inverse'), name='Period') +
scale_y_continuous(name='Power', trans='log10')

Following my nose further on the logistic map example, I also used the

**animate** package to make a movie that walks through the power spectrum of the map for a range of r values. Maybe one day I'll set this to music and post it to Youtube! Though I think some strategic speeding up and slowing down in important parameter regions is warranted for a truly epic tour of the logistic map's power spectrum.

require(animation)
saveVideo({
ani.options(interval = 1/5, ani.type='png', ani.width=1600, ani.height=1200)
for (ii in c(seq(from=3.2, to=3.5, length.out=200), seq(from=3.5, to=4, length.out=1200))) {
spectrum(logmap(1e4, ii), taper=0, demean=T, pad=0, fast=T, sub=round(ii, 3))
}
}, video.name = "logmap2.mp4", other.opts = "-b 600k")