What is the frequency domain?
Fourier transforms
Spectral analysis
Wavelets
23 May 2023
Topics for today
Time domain
We having been examining changes in \(x_t\) over time
Time domain
We can think of this as comparing changes in amplitude (displacement) with time
Frequency domain
Today we’ll consider how amplitude changes with frequency
Jean-Baptiste Fourier (1768 - 1830)
French mathematician & physicist best known for his studies of heat transfer
First described what we now call the “greenhouse effect”
Solving hard problems
Solving the heat equation involves solving partial differential equations conditional on some boundary conditions
\[ \begin{matrix} \large \text{Problem} \\ really \Bigg \downarrow hard \\ \large \text{Solution} \end{matrix} \]
Fourier’s approach
Find \(f(t)\) and \(\hat{f}(t)\), such that
\[ \begin{matrix} \large \text{Problem} && \xrightarrow{f(t)} && \large \text{Transformed problem} \\ really \Bigg \downarrow hard && && much \Bigg \downarrow easier \\ \large \text{Solution} && \xleftarrow{\hat{f}(t)} && \large \text{Transformed solution} \\ \end{matrix} \]
Fourier series
Complex periodic functions can be written as infinite sums of sine waves
\[ f(t) = a_0 + \sum_{k = 1}^\infty a_k \sin(2 \pi f_0 k t + p_k) \]
where
\(k\) is the wave number (index)
\(a_k\) is the amplitude of wave \(k\)
\(f_0\) is the fundamental frequency
\(p_k\) is the phase shift
Fourier series
A finite example
\[ f(t) = \sum_{k = 1}^5 \frac{1}{k} \sin(2 \pi k t + k^2) \]
Fourier series
Fourier series
Fourier series
Here’s an animated example from Wikipedia
Fourier transform
We can make use of Euler’s formula
\[ \cos(2 \pi k) + i \sin(2 \pi k) = e^{i 2 \pi k} \]
and write the Fourier transform of \(f(t)\) as
\[ f(t) = \int_{-\infty}^\infty \hat{f}(k) ~ e^{i 2 \pi t k} ~ d k \]
where \(k\) is the frequency
Discrete Fourier transform
Fourier transform
\[ f_k = \sum_{n=0}^{N-1} x_t ~ e^{-i 2 \pi n k} \]
Fourier transforms in R
R uses what’s known as Fast Fourier transform via fft(), which returns the amplitude at each frequency
ft <- fft(xt) ## often normalize by the length ft <- fft(xt) / length(xt)
Fourier represention of our \(\{x_t\}\)
Discrete Inverse Fourier transform
Fourier transform
\[ f_k = \sum_{n=0}^{N-1} x_t ~ e^{-i 2 \pi n k} \]
Inverse
\[ x_t = \sum_{k=0}^{N-1} f_k ~ e^{i 2 \pi n k} \]
Inverse Fourier transforms in R
i <- complex(1, re = 0, im = 1)
xx <- rep(NA, TT)
kk <- seq(TT) - 1
## Inverse Fourier transform
## ft <- fft(xt)
for(t in kk) {
xx[t+1] <- sum(ft * exp(i*2*pi*kk*t/TT))
}
Original \(\{x_t\}\) & our inverse transform
Inverse Fourier transforms in R
ift <- fft(ft, inverse = TRUE)
Original \(\{x_t\}\) & R’s inverse transform
Spectral analysis
Spectral analysis
Spectral analysis refers to a general way of decomposing time series into their constituent frequencies
Spectral analysis
Consider a linear regression model for \(\{x_t\}\) with various sines and cosines as predictors
\[ x_t = a_0 + \sum_{k = 1}^{n/2-1} a_k \cos(2 \pi f_0 k t/n) + b_k \sin(2 \pi f_0 k t/n) \]
Periodogram
The periodogram measures the contributions of each frequency \(k\) to \(\{x_t\}\)
\[ P_k = a^2_k + b^2_k \]
Estimate the periodogram in R
spectrum(xt, log = "on") spectrum(xt, log = "off") spectrum(xt, log = "dB")
Periodogram for our \(\{x_t\}\)
spectrum(xt, log = "dB")
Periodogram for our \(\{x_t\}\)
Density on natural scale & frequency in cycles per time
Spectral density estimation via AR(p)
For an AR(p) process
\[ x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \cdots + \phi_p x_{t-p} + e_t \]
The spectral density is
\[ S(f,\phi_1, \dots, \phi_p, \sigma^2) = \frac{\sigma^2 \Delta t}{\vert 1 - \sum_{k = 1}^p \phi_k e^{-i 2 \pi f k \Delta t} \vert ^2} \]
Limits to spectral analysis
Spectral analysis works well for
stationary time series
identifying periodic signals corrupted by noise
Limits to spectral analysis
Spectral analysis works well for
stationary time series
identifying periodic signals corrupted by noise
But…
it’s an inconsistent estimator for most real data sets
it’s generally biased
Wavelets
Shifting frequencies
What if the frequency changes over time?
Wavelets
For non-stationary time series we can use so-called wavelets
A wavelet is a function that is localized in time & frequency
Graphical forms for decomposition
Graphical form for decomposition
What is a wavelet?
Formally, a wavelet \(\psi\) is defined as
\[ \psi_{\sigma, \tau}(t) = \frac{1}{\sqrt{\vert \sigma \vert}} \psi \left( \frac{t - \tau}{\sigma} \right) \]
where \(\tau\) determines its position & \(\sigma\) determines its frequency
Graphical form for decomposition
Properties of wavelets
It goes up and down
\[ \int_{-\infty}^{\infty} \psi(t) ~ dt = 0 \]
It has a finite sum
\[ \int_{-\infty}^{\infty} \vert \psi(t) \vert ~ dt < \infty \]
How are wavelets defined?
In terms of scaling functions that describe
Dilations \(~~~~~~~~~~~~~~ \psi(t) \rightarrow \psi(2t)\)
Translations \(~~~~~~~~~ \psi(t) \rightarrow \psi(t - 1)\)
How are wavelets defined?
More generally,
\[ \psi_{j,k}(t) = 2^{j/2} \psi(2^j t - k) \]
where
\(j\) is the dilation index
\(k\) is the translation index
and
\(2^{j/2}\) is a normalization constant
Wavelets in practice
There are many options for \(\psi(t)\), but we’ll use scaling functions and define
\[ \psi(t) = \sum_{k=0}^K c_k \psi(2x - k) \]
where the \(c_k\) are filter coefficients*
*Note that \(\psi(t)\) gets “smoother” as \(K\) increases
Haar’s scaling function
Simple, but commonly used, where \(\small K = 1; ~ c_0 = 1; ~ c_1 = 1\)
\[ \psi(t) = \sum_{k=0}^K c_k \psi(2t - k) \\ \big \Downarrow \\ \psi(t) = \psi(2t) + \psi(2t - 1) \]
The only function that satisfies this is:
\[ \begin{align} \psi(t) &= 1 ~ \text{if} ~ 0 \leq t \leq 1 \\ \psi(t) &= 0 ~ \text{otherwise} \end{align} \]
Haar’s scaling function
Haar’s scaling function
In terms of the dilation
\[ \begin{align} \psi(2t) &= 1 ~ \text{if} ~ 0 \leq t \leq 0.5 \\ \psi(2t) &= 0 ~ \text{otherwise} \end{align} \]
and translation
\[ \begin{align} \psi(2t - 1) &= 1 ~ \text{if} ~ 0.5 \leq t \leq 1 \\ \psi(2t - 1) &= 0 ~ \text{otherwise} \end{align} \]
Haar’s scaling function (father)
Haar’s mother wavelet
Wavelets are created via differencing of scaling functions
\[ \psi(t) = \sum_{k=0}^1 (-1)^k c_k \psi(2t - k) \]
where \((-1)^k\) creates the difference
Haar’s mother wavelet
Family of Haar’s wavelets
So-called “child” wavelets are created via dilation & translation
\[ \psi_{j,k}(t) = 2^{j/2} \psi(2^j t - k) \]
The mother Haar wavelet has \(j = 0\)
Family of Haar’s wavelets
So-called “child” wavelets are created via dilation & translation
\[ \psi_{j,k}(t) = 2^{j/2} \psi(2^j t - k) \]
The basic Haar wavelet has \(j = 0\)
Setting \(j = 1\) yields a daughter
\[ \psi_{j,k}(t) = \sqrt{2} \psi(2 t - k) \]
Haar’s daughter wavelet
\[ \psi(t) = \sum_{k=0}^1 (-1)^k c_k \sqrt{2} \psi(2 t - k) \]
Recall that \((-1)^k\) creates the difference
A daughter wavelet of Haar’s
Other wavelets
There are many forms of wavelets, many of which were developed in the past 50 years
Morlet
Mexican Hat
Who does this?
Wavelet analysis is used widely in audio & video compression
JPEG
Estimating wavelet transforms in R
We’ll use the WaveletComp package, which uses the Morlet wavelet
We’ll also use the L Washington temperature data from the MARSS package
library(WaveletComp) ## L WA temperature data tmp <- MARSS::lakeWAplanktonTrans[,"Temp"] ## WaveletComp needs data as df dat <- data.frame(tmp = tmp)
Estimating wavelet transforms in R
Use analyze.wavelet() to estimate the wavelet transform
w_est <- analyze.wavelet(dat, "tmp", ## need both df & colname
loess.span = 0, ## no de-trending
dt = 1/12, ## monthly sampling
lowerPeriod = 1/6, ## default = 2*dt
n.sim = 100,
verbose = FALSE)
## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
Estimating wavelets in R
Use wt.image() to plot the spectrum
Inverse wavelet transforms
Involves integral calculus
\[ f(t) = \frac{1}{C_{\psi}} \int_a \int_b < f(t), \psi_{a,b}(t) > \psi_{a,b}(t) db \frac{da}{a^2} \]
Inverse wavelet transforms in R
Use reconstruct() to get estimate of original time series