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.
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.
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:
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.
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)\}\).
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):
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:
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:
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:
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
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.
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.