Wavelet analysis; significance levels; confidence intervals;

Recently, there is an upsurge of the use of wavelet analysis for reaction time data. Hilde Geurts wanted to apply this technique to follow up a spectral analysis RT time series, for comparisson between different clinical groups. A good introduction, including the practical statistical issues of wavelet analysis, can be found on this web site: Wavelet analysis; significance levels; confidence intervals;. (In particular: download the pdf "A practical guide to wavelet analysis", published in the Bulletin of the American Meteorological Assocation.) There are several contributed packages for wavelet analysis in R. One of them, Rwave, provides the function cwt for the calculation of the continuous wavelet transform using Morlet wavelets, as discussed on web site indicated above. A difficulty with wavelets analysis is that, although wavelets are often explained as providing a time-frequency decomposistion of a time series, the appropriateness of this assertion realy depends on the wavelet basis function used (the so called Mother Wavelet). This interpretation can be safely assumed when Morlet basis wavelet is used—it should be kept in mind though, that on the lower end of the frequency spectrum the continuos wavelet transform provides little time resolution but high frequency resolution, whereas on the higher end of the spectrum the cwt provides little frequency resolution but high time resolution. Because the 'frequency spectrum' interpretation is not appropriate for all types of wavelets, and because of the different scales of resolution, it is custom to display the wavelet transform along a time and a 'scale' axis. However, for those accustomed to frequency spectra, this is not so easy to interpret, and it is desirable to display an appropriate frequency axis in stead of the scale access. Unfortunately this is not implemented in the Rwave package. Here is a function that will display the CWT (computed with cwt) with a frequency axis:

freqwave <-
function (x, a, noctave = log2(nextn(length(x))), nvoice = 12,
    xlab = 'time', ylab = 'frequency (Hz)', ...)
{
    if (any(is.na(a)))
        a <- cwt(x, noctave, nvoice, plot = FALSE)
    T = deltat(x) * (length(x) - 1)
    f = pretty(log2(nextn(length(x))) - 1
        - (0:(nvoice * noctave))/nvoice)
    lt = pretty(time(x))
    lf = (1/T) * 2^pretty(rev(log2(nextn(length(x)))
        - 1 - (0:(nvoice * noctave))/nvoice))
    f = pretty(rev(log2(nextn(length(x))) - 1
        - (0:(nvoice * noctave))/nvoice))
    bt = 1/diff(range(time(x)))
    at = -bt * min(time(x))
    bf = -1/diff(range(rev(log2(nextn(length(x)))
         - 1 - (0:(nvoice * noctave))/nvoice)))
    af = -bf * max(rev(log2(nextn(length(x))) - 1
         - (0:(nvoice * noctave))/nvoice))
    filled.contour(Mod(a), xlab = xlab, ylab = ylab,
        plot.axes = { axis(1, at + bt * lt, lab = lt)
        axis(2, af + bf * f, lab = format(lf, dig = 1))
    }, ...)
    invisible(list(t = time(x), f = lf, cwt = a))
}

It accepts a time series object as it first argument, and a second argument which is either NA or the complex matrix returned by cwt. If the second argument is NA, the CWT will be computed. Note that the frequency axis only makes sense if the sampling frequency, or deltat, of the time series object is set. As an example,

library(Rwave)
t = seq(0,1,len=512)
x = 2 * sin(2*pi*16*t)*exp(-(t-.25)^2/.001)
x = x + sin(2*pi*64*t)*exp(-(t-.75)^2/.001)
x = ts(x,deltat=1/512) # sampling frequency is 512 Hz
freqwave(x,NA)

The result is this: