What is this all about?

Many real signals are/should be viewed as inherently “nonstationary” (and cannot be made “stationary” through a simple transformation). For example, the following audio signal of “How are you?” by its nature consists of different sounds at different times.

library(Rwave)

data(HOWAREYOU)
ts.how <- HOWAREYOU

library(audio)
#play(ts.how,rate=2800/.35)

plot(ts.how,type="l")

The usual “stationary” analysis of this signal, for example, through the raw periodogram given below, might be misleading, as it assumes that the signal looks alike in different time segments.

spec.pgram(ts.how,log='no',taper=0,pad=0,fast=FALSE,demean=TRUE,detrend=FALSE)

There are natural modifications of such spectral (frequency) analysis that consider frequency content locally in time. They make part of what is known as time-frequency analysis (sometimes time-frequency/time-scale analysis), to be discussed in this lecture.

Data/Theory goals of the lecture

The goal of the lecture is to provide the basics of time-frequency analysis. We will start with some general observations, which will then be specialized to the so-called windowed Fourier transform (WFT) and continuous wavelet transform (CWT). The purposes of WFT and CWT are slightly different, but they still fit naturally under the umbrella of time-frequency or time-frequency/time-scale analysis

The concepts will be illustrated using the \(\verb!R!\) package \(\verb!Rwave!\), associated with the monograph on the subject by Carmona et al. (1998). In fact, time-frequency analysis has historically been viewed as part of signal processing, where \(\verb!Matlab!\) is often preferred, and has seen relatively few developments on the \(\verb!R!\) front.

Basics on Fourier transforms

To present the theory of time-frequency analysis, it is more convenient to switch to continuous time \(t,u,\ldots\in \mathbb{R}\). The continuous analogue of the discrete Fourier transform behind the periodogram above is the (continuous) Fourier transform (FT) of function \(f:\mathbb{R} \to \mathbb{R}\ (\mbox{or}\ \mathbb{C})\), defined by \[ {\cal F}_f(w) = \widehat f(w) = \int_{\mathbb{R}} e^{-iwt} f(t) dt = \int_{\mathbb{R}} \cos(wt) f(t) dt - i\int_{\mathbb{R}} \sin(wt) f(t) dt,\quad w\in\mathbb{R}. \] This is defined for functions \(f\in L^2(\mathbb{R})\) and the resulting Fourier transform \(\widehat f \in L^2(\mathbb{R})\). But such “technical”" details will be excluded from the lecture, unless absolutely necessary.

Few important things to keep in mind about FT. The following formulas hold:

  • Inversion formula: \[ f(t) = \frac{1}{2\pi} \int_{\mathbb{R}} e^{iwt} \widehat f(w) dw. \]
  • Parseval formula: \[ \int_{\mathbb{R}} f(t) \overline{g(t)} dt = \frac{1}{2\pi} \int_{\mathbb{R}} \widehat f(w) \overline{\widehat g(w)} dw. \] (It is called Plancherel’s when \(f=g\).)

The inversion formula states that the signal \(f(t)\) is viewed as a weighted superposition of “atoms” \(g_{(w)}(t) = e^{iwt}\) across \(w\in \mathbb{R}\), where the weights are given by the FT of the signal. The atoms \(g_{(w)}(t) = e^{iwt}\) are “localized” at frequency \(w\) but not in time. We shall generalize this to consider atoms that are potentially “localized” in time as well.

Atomic decompositions (Continuous bases)

One thinks of a signal \(f(t)\) as consisting of a weighted superposition of atoms \(g_{(u,\xi)}(t)\) that are “localized”" in time \(u\) and frequency \(\xi\). In mathematical terms, this is expressed as \[ f(t) = \mathop{\int\int}_{\mathbb{R}^2} L_f(u,\xi) g_{(u,\xi)}(t) \mu(du,d\xi), \] where \(\mu(du,d\xi)\) is a suitable measure and the weights are given by \[ L_f(u,\xi) = \int_{\mathbb R} f(t) \overline{g_{(u,\xi)}(t)} dt = \frac{1}{2\pi} \int_{\mathbb R} \widehat f(w) \overline{\widehat g_{(u,\xi)}(w)} dw. \]

If \(g_{(u,\xi)}\) is localized at some time and some frequency, then \(L_f(u,\lambda)\) carries information about \(f\) and \(\widehat f\) simultaneously about that time and that frequency. This is conveniently represented in the so-called time frequency plane by using the Heisenberg boxes (or time-frequency spreads) of atoms \(g_{(u,\xi)}\).

Heisenberg (uncertainty) principle: It is not very difficult to show that supposing \(\|g\|_2=1\), \(\sigma_t \sigma_w\geq 1/2\) (with the equality for Gaussian density function \(g\)).

Closure relation: For the above to work, the atoms have to satisfy the closure relation (resolution of the identity) \[ \mathop{\int\int}_{\mathbb{R}^2} g_{(u,\xi)}(t) \overline{g_{(u,\xi)}(t')} \mu(du,d\xi) = \delta(t-t'), \] where \(\delta\) is the delta function at \(0\) (to be discussed in special cases below).

Another consequence of the atomic decomposition is that \[ L_f(u',\xi') = \mathop{\int\int}_{\mathbb{R}^2} K(u,u';\xi,\xi') L_f(u,\xi) \mu(du,d\xi), \] where \(K\) is the “reproducing kernel” given by \[ K(u,u';\xi,\xi') = \int_{\mathbb R} g_{(u,\xi)}(t) \overline{g_{(u',\xi')}(t)} dt = \frac{1}{2\pi} \int_{\mathbb R} \widehat g_{(u,\xi)}(w) \overline{\widehat g_{(u',\xi')}(w)} dw. \] This expresses the redundancy in the information provided by \(\{L_f(u,\xi)\}\).

Special case 1: Windowed Fourier transform (WFT)

The case \[ g_{(u,\xi)}(t) = e^{i\xi t} g(t-u) \] with real, symmetric \(g\) satisfying \(\|g\|_2=1\), or equivalently \[ \widehat g_{(u,\xi)}(w) = e^{-iu(w-\xi)} \widehat g(w-\xi), \] is known as the windowed Fourier transform (WFT); also short-term Fourier transform (STFT), time dependent Fourier transform, continuous Gabor transform, etc. The corresponding WFT is \[ S_f(u,\xi) = \int_{\mathbb R} f(t) e^{-i\xi t} g(t-u) dt = \frac{1}{2\pi} \int_{\mathbb R} \widehat f(w) e^{iu(w-\xi)} \widehat g(w-\xi) dw. \]

Some of the “windows” \(g\) used (without normalization):

  • Rectangular: \(g(t) = 1_{[-\frac{1}{2},\frac{1}{2}]}(t)\)
  • Hanning: \(g(t) = \cos^2(\pi t) 1_{[-\frac{1}{2},\frac{1}{2}]}(t)\)
  • Hamming: \(g(t) = (0.54 + 0.46 \cos(2\pi t)) 1_{[-\frac{1}{2},\frac{1}{2}]}(t)\)
  • Gaussian: \(g(t) = e^{-t^2/2\sigma^2}\)

Note: The atoms satisfy the closure relation (why?), inversion formula holds, there is a reproducing kernel, etc. The Heisenberg boxes look like:

Spectrogram: The spectrogram is defined as \(P_s(u,\xi) = |S_f(u,\xi)|^2\). The same term is also used for the corresponding plot. (For audio, it is also called sonogram.)

Note: In practice, discrete versions of continuous integrals are used.

Example 2.1 (Chirp): For \[ f(t) = e^{iat^2}\quad (a>0) \] and the Gaussian window, one can show that \[ P_s(u,\xi) = C(\sigma,a) e^{-\frac{\sigma^2(\xi-2au)^2}{1+4a^2\sigma^2}}. \] Note that it reaches its maximum at \(\xi(u) = 2au = (au^2)'\), which is called the instantaneous frequency.

x <- 1:512
chirp <- sin(2*pi*(x + 0.002*(x-256)*(x-256))/16)
plot (chirp,type="l")
title("Chirp signal")

#play(chirp,rate=2500)

nvoice <- 50
freqstep <- .005
scale <- 25
#scale <- 100
#scale <- 5
isize <- length(chirp)
cgtchirp <-cgt(chirp,nvoice,freqstep,scale,plot=FALSE)
image(1:isize, seq(0, nvoice * freqstep/2, length = nvoice), Mod(cgtchirp), xlab = "Time", ylab = "Frequency",col = grey(seq(0, 1, length = 256)))
        title("Gabor Transform Modulus")

Example 2.2 (Speech/Audio signal):

plot(ts.how,type="l")

nvoice <- 70
freqstep <- 1/70
scale <- 60
isize <- length(ts.how)
cgtH0WAREYOU<- cgt(ts.how,nvoice,freqstep,scale,plot=FALSE)
image(1:isize, seq(0, nvoice * freqstep/2, length = nvoice), Mod(cgtH0WAREYOU), xlab = "Time", ylab = "Frequency",col = grey(seq(0, 1, length = 256)))
        title("Gabor Transform Modulus")

Example 2.3 (Backscattered acoustic signals):

# Acoustic returns from an underwater metallic object
data("back1.000")
# Acoustic returns from natural underwater clutter
data("back1.220")

ts.back1 <- back1.000
ts.back2 <- back1.220

par(mfrow=c(1,2))
plot(ts.back1,type="l")
plot(ts.back2,type="l")

nvoice <- 50
freqstep <- 0.005
scale <- 50
isize <- length(ts.back1)
cgt.back1 <- cgt(ts.back1-mean(ts.back1),nvoice,freqstep,scale,plot=FALSE)
image(1:isize, seq(0, nvoice * freqstep/2, length = nvoice), Mod(cgt.back1), xlab = "Time", ylab = "Frequency",col = grey(seq(0, 1, length = 256)))
        title("Gabor Transform Modulus")
        
nvoice <- 50
freqstep <- 0.005
scale <- 50
isize <- length(ts.back2)
cgt.back2 <- cgt(ts.back2-mean(ts.back2),nvoice,freqstep,scale,plot=FALSE)
image(1:isize, seq(0, nvoice * freqstep/2, length = nvoice), Mod(cgt.back2), xlab = "Time", ylab = "Frequency",col = grey(seq(0, 1, length = 256)))
        title("Gabor Transform Modulus")

Other issues to discuss:

  • Shape and size of window
  • Phase plot
  • …

Special case 2: Continuous wavelet transform (CWT)

This special case is associated with the atoms \[ g_{(u,s)}(t) = \frac{1}{\sqrt{s}} \psi\Big( \frac{t-u}{s} \Big), \] where \(s>0\) is referred to as a scale, \(\psi\) is known as a wavelet (either real-valued or complex-valued). Equivalently, \[ \widehat g_{(u,s)}(w) = e^{-iuw} s^{1/2} \widehat \psi(sw). \] The corresponding transformation \[ W_f(u,s) = \int_{\mathbb R} f(t) \overline{\frac{1}{\sqrt{s}} \psi\Big( \frac{t-u}{s} \Big)} dt = \frac{1}{2\pi} \int_{\mathbb R} \widehat f(w) \overline{e^{-iuw} s^{1/2} \widehat \psi(sw)} dt. \] is called a continuous wavelet transform.

Some of the wavelets \(\psi\) used:

  • Mexican hat (Marr’s): \(\psi(t) = - \frac{\partial^2}{\partial t^2} e^{-t^2/\sigma^2}\)
  • Difference of Gaussians: \(\psi(t) = 2e^{-2t^2} - e^{-t^2/2}\)
  • Gabor/Morlet (complex-valued): \(\psi(t) = e^{-t^2/2\sigma^2} e^{i\eta t}\)
mor <- morlet(512,256,40)
par(mfrow=c (1,2))
plot(Re(mor),type="l")
title("Real part")
plot(Im(mor),type="l")
title ("Imaginary part")

Note: The atoms satisfy the closure relation (why?), but need (in the real-valued case) \(\int_{\mathbb R} \psi(t) dt =0\). The Heisenberg boxes look like:

Scalogram: The scalogram is defined as \(P_Wf(u,\xi) = |W_f(u,s)|^2 = |W_f(u,\frac{\mu_{w,\psi}}{\xi})|^2\) with \(\xi = \frac{\mu_{w,\psi}}{s}\). The same term is also used for the corresponding plot.

Note: In practice, discrete versions are used.

CWT is particularly suitable for point singularities. (Why?) (In contrast, WFT is for local periodicities.)

Example 2.4 (Some singularities):

sing <- numeric(256)
sing[30:60] <- sqrt(1:30)/sqrt(30)
sing [60:90] <- 1
sing[140] <- 0.5
sing[200] <- 0.5
sing[201] <- -0.5
par(mfrow=c(1,2))
plot(sing,type="l")
title("Some singularities")


isize <- length(sing)
noctave <- 5
nvoice <- 12
pp <- noctave * nvoice
cwtsing <-cwt(sing,noctave,nvoice,plot=FALSE)
#persp(Mod(cwtsing))
#image(Arg(cwtsing))
image(1:isize, seq(0, pp, length = pp), Mod(cwtsing), xlab = "Time", ylab = "log(scale)",main = "Wavelet Transform Modulus",col = grey(seq(0, 1, length = 256)))

Example 2.5 (Gravitational waves):

data(purwave)
ts.wave <- purwave

plot(ts.wave,type="l")

plot(ts.wave[7350:7510],type="l")

isize <- length(ts.wave)
noctave <- 5
nvoice <- 8
pp <- noctave * nvoice
cwtpure <- cwt(ts.wave, noctave, nvoice, plot=FALSE)
image(1:isize, seq(0, pp, length = pp), Mod(cwtpure), xlab = "Time", ylab = "log(scale)",main = "Wavelet Transform Modulus",col = grey(seq(0, 1, length = 256)))

data(noisywave)
ts.wave.noisy <- noisywave

plot(ts.wave.noisy,type="l")

cwtnoisy <- cwt(ts.wave.noisy, noctave, nvoice,  plot=FALSE)
image(1:isize, seq(0, pp, length = pp), Mod(cwtnoisy), xlab = "Time", ylab = "log(scale)",main = "Wavelet Transform Modulus",col = grey(seq(0, 1, length = 256)))

Other issues to discuss:

  • Other plots
  • …

Discussion on applications

  • Bird sonograms

  • Time frequency functional models (e.g. Holan et al. (2010), (2012))

  • Biophysical time series (EEG ,etc)

  • Geophysics (earthquakes, etc)

  • Next lecture (for wavelets): denoising, smoothness

Some final notes

  • Other classical textbooks: Daubechies (1992), Mallat (1998), Flandrin (1999), etc.
  • Extensions to \(2\)D or higher D.

Questions/Homework

Q0: Try some of the considered time-frequency/time-scale analysis methods on real time series.

Q1: Prove the Heisenberg principle.

Q2: Let \(f(t) = e^{i\phi(t)}\). Show that \(\int_{\mathbb R} \xi |S_f(u,\xi)|^2 d\xi = 2\pi \int_{\mathbb R} \phi'(t) |g(t-u)|^2 dt\) and interpret this result.

Q3: Let \(\psi\) be a real wavelet and \(\phi\) be a real function such that \[ |\widehat \phi(w)|^2 = \int_1^\infty |\widehat \psi(sw)|^2 \frac{ds}{s} = \int_w^\infty |\widehat \psi(\xi)|^2 \frac{d\xi}{\xi},\quad w>0. \]

Set \[ L_f(u,s) = \int_{\mathbb R} f(t) \overline{\frac{1}{\sqrt{s}}\, {\phi\left(\frac{t-u}{s}\right)}} dt,\quad u\in \mathbb{R},s>0. \] Without caring about technical issues (e.g. integrability), show that \[ f(t) = C_\psi^{-1} \left\{\frac{1}{a} \int_{\mathbb R} L_f(u,a) \frac{1}{\sqrt{a}}\, {\phi\left(\frac{t-u}{a}\right)} du+ \int_{\mathbb R} \int_0^a W_f(u,s) \frac{1}{\sqrt{s}}\, {\psi\left(\frac{t-u}{s}\right)} \frac{duds}{s^2} \right\} \] \[ = C_\psi^{-1} \lim_{a\to 0} \frac{1}{a}\int_{\mathbb R} L_f(u,a) \frac{1}{\sqrt{a}}\, {\phi\left(\frac{t-u}{a}\right)} du, \]

where \(W_f\) is a continuous wavelet transform. (Note: The function \(\phi\) is called a scaling function. The first term in the braces above is interpreted as an approximation term of \(f\) at a scale \(a\), and the second term as the details added over the scales smaller than \(a\).)

Q4: Explain what time frequency functional models are.

References

Carmona, R., W.-L. Hwang, and B. Torrésani. 1998. Practical Time-Frequency Analysis. Vol. 9. Wavelet Analysis and Its Applications. Academic Press, Inc., San Diego, CA.

Daubechies, I. 1992. Ten Lectures on Wavelets. Vol. 61. CBMS-Nsf Regional Conference Series in Applied Mathematics. Society for Industrial; Applied Mathematics (SIAM), Philadelphia, PA.

Flandrin, P. 1999. Time-Frequency/Time-Scale Analysis. Vol. 10. Wavelet Analysis and Its Applications. Academic Press, Inc., San Diego, CA.

Holan, S. H., C. K. Wikle, L. E. Sullivan-Beckers, and R. B. Cocroft. 2010. “Modeling Complex Phenotypes: Generalized Linear Models Using Spectrogram Predictors of Animal Communication Signals.” Biometrics 66 (3): 914–24.

Holan, S. H., W.-H. Yang, D. S. Matteson, and C. K. Wikle. 2012. “An Approach for Identifying and Predicting Economic Recessions in Real-Time Using Time-Frequency Functional Models.” Applied Stochastic Models in Business and Industry 28 (6): 485–99.

Mallat, S. 1998. A Wavelet Tour of Signal Processing. Academic Press, Inc., San Diego, CA.