# Fourier Methods I

2 October 2018

$\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}\left[ #1 \right]} \newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]} \newcommand{\TrueRegFunc}{\mu} \newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator{\det}{det} \newcommand{\TrueNoise}{\epsilon} \newcommand{\EstNoise}{\widehat{\TrueNoise}} \newcommand{\SignalNoise}{N}$

# In our previous episodes: Decomposing data into components

• Eigenvectors of the smoother
• what’s kept by smoothing?
• what shows up in the residuals?
• Eigenvectors of the variance matrix for multivariate data
• what are the principal components?
• == what are the directions of maximum variance?
• == what are the uncorrelated components?
• == what components should we keep when reducing dimension?

# Today: Decomposing data into periodic components

• Helps with periodic trends/seasonality
• Also helps with covariance estimation
• And with any sort of linear operation

# Sine waves

We get data at times $$1, 2, \ldots n$$

curve(sin(2 * pi * x/n), xlim = c(0, n), ylim = c(-2, 2), type = "p", pch = 16,
cex = 0.5, col = "blue")

# Sine waves

curve(sin(2 * pi * x/n), xlim = c(0, n), ylim = c(-2, 2), type = "p", pch = 16,
cex = 0.5, col = "blue")
curve(sin(2 * pi * 3 * x/n), type = "b", pch = 16, cex = 0.5, col = "yellow3",
add = TRUE)

# Sine waves

curve(sin(2 * pi * x/n), xlim = c(0, n), ylim = c(-2, 2), type = "p", pch = 16,
cex = 0.5, col = "blue")
curve(sin(2 * pi * 3 * x/n), type = "p", pch = 16, cex = 0.5, col = "yellow3",
curve(sin(2 * pi * x/n) + sin(2 * pi * 3 * x/n), type = "p", pch = 16, cex = 0.5,
col = "green", add = TRUE)

# An ambition

• Examples like this suggest a new project: decompose arbitrary time series into a sum of sine waves
• i.e., write anything, periodic or not, as a sum of periodic components

# Some properties of sine waves

$\begin{eqnarray} \sin{(2\pi \omega \frac{t}{n})} & \text{is} & \text{periodic, with period} ~ n/\omega ~ \text{and frequency} ~\omega/n\\ \sin{-\theta} & = & -\sin{\theta}\\ \sin{0} & = & 0\\ \sin{\pi/2} & = & 1\\ \sin{\theta+\pi} & = & -\sin{\theta} \end{eqnarray}$

Consequences:

• If we just use a finite sum of sine waves, we’ll get periodic functions
• But we could match a non-periodic function on the domain where we see it
• If we just use $$\sin{2\pi \omega t/n}$$ functions, we can only fit data which goes through the origin as is mid-way between peak and trough at $$t=0$$

# Some properties of cosine waves

$\begin{eqnarray} \cos{(2\pi \omega \frac{t}{n})} & \text{is} & \text{periodic, with period} ~ n/\omega ~ \text{and frequency} ~\omega/n\\ \cos{-\theta} & = & \cos{\theta}\\ \cos{0} & = & 1\\ \cos{\pi/2} & = & 0\\ \cos{\theta+\pi} & = & -\cos{\theta} \end{eqnarray}$

Consequences:

• If we just use a finite sum of cosine waves, we’ll get periodic functions
• But we could match a non-periodic function on the domain where we see it
• Only using cosine waves only lets us fit data which is at a peak or trough at $$t=0$$

# Some properties of combining sines and cosines

$\begin{eqnarray} \sin{\theta+\phi} & = & \sin{\theta}\cos{\phi} + \sin{\phi}\cos{\theta}\\ \sin{\theta+\pi/2} & = & \cos{\theta}\\ \cos{\theta+\phi} & = & \cos{\theta}\cos{\phi} - \sin{\theta}\sin{\phi}\\ \cos{\theta-\pi/2} & = & \sin{\theta} \end{eqnarray}$

$$\Rightarrow$$ If we use $$\sin{(\phi_n + 2\pi \omega t/n)}$$ functions, we don’t need cosines, and vice versa

# Sines average out to zero

For any integer $$\omega$$, $\begin{eqnarray} \sum_{t=0}^{n-1}{\sin{2\pi\omega \frac{t}{n}}} & = & 0\\ \end{eqnarray}$ because $\sin{2\pi\omega \frac{t}{n}} = -\sin{2\pi\omega \frac{t+n/2\omega}{n}}$ i.e., time points come in pairs which cancel each other out in the sum

• ditto for averages of cosines $\sum_{t=0}^{n-1}{\cos{2\pi\omega\frac{t}{n}}} = 0 ~ \text{unless} ~ \omega =0$

# Sines of different frequencies are orthogonal

• The inner product of two sine waves is usually 0: $\begin{eqnarray} \sum_{t=1}^{n}{\sin{\left(2\pi\omega_1 \frac{t}{n}\right)} \sin{\left(2\pi\omega_2 \frac{t}{n}\right)}} & = & \frac{1}{2}\sum_{t=1}^{n}{\cos{\left(2\pi(\omega_1-\omega_2) \frac{t}{n}\right)} - \cos{\left(2\pi(\omega_1+\omega_2)\frac{t}{n}\right)}}\\ & = & 0 ~ \text{unless} ~ \omega_1=\omega_2 ~ \text{or}~ \omega_1 = -\omega_2 \end{eqnarray}$

• If $$\omega_1 = \omega_2 \neq 0$$, we get $$n/2$$
• If $$\omega_1 = -\omega_2 \neq 0$$, we get $$-n/2$$
• If $$\omega_1 = \omega_2 = 0$$, we get 0
• ditto for inner products of cosines

# We can’t resolve either very high or very low frequency sine waves

• High frequency: If $$\sin{\phi + 2\pi\omega t/n}$$ goes through a complete cycle between $$t$$ and $$t+1$$, then it looks exactly like $$\sin{\phi}$$ everywhere
• Low frequency: If $$\omega$$ is too small, we don’t see even half a cycle, and we can’t tell the difference between that and a constant plus higher-frequency terms

# Orthogonal vectors give us a basis

• Remember from linear algebra: vectors $$\vec{v}_1, \vec{v}_2, \ldots \vec{v}_n$$ are a basis when, for any vector $$\vec{x}$$, $\vec{x} = \sum_{i=1}^{n}{c_i \vec{v}_i}$
• In an $$n$$-dimensional space, any set of $$n$$ orthogonal vectors is a basis
• Also true of any $$n$$ linearly independent vectors
• So we can use the sines and cosines as a basis
• Not normalized, but that’s OK

# Re-writing the data in terms of sines and cosines

• Re-write any vector in terms of the basis: $x(t) = \sum_{\omega=0}^{n/2-1}{a_{\omega}\sin{\left(2\pi\omega \frac{t}{n}\right)} + b_{\omega}\cos{\left(2\pi\omega\frac{t}{n}\right)}}$

• Why only up to $$n/2-1$$?
• $$b_0 = \overline{x}$$ (why?) and $$a_0=0$$ (why?)
• The $$a$$ and $$b$$ coefficients are linked…

# Re-writing the data in terms of sines

Use the trig identities to eliminate cosines in favor of sines:

$x(t) = \sum_{\omega=0}^{n/2-1}{\alpha_{\omega} \sin{\left(2\pi\frac{t}{n} + \phi_{\omega}\right)}}$ with $\alpha_{\omega} = \sqrt{a_{\omega}^2+b_{\omega}^2}$

# Re-writing the data in terms of complex exponentials

• $$e^{i \theta} = \cos{\theta} + i \sin{\theta}$$
• Re-write to eliminate $$\sin$$ in favor of exponentials: $x(t) = \sum_{\omega=0}^{n-1}{c_{\omega} e^{2\pi i \omega t/n}}$

# The Fourier transform

• The Fourier transform of $$x(t)$$, $$t=0, \ldots n-1$$, is $\tilde{x}(\omega) \equiv \sum_{t=0}^{n-1}{x(t) e^{-2\pi i t \omega/n}}$
• These numbers are the Fourier coefficients
• $$\tilde{x}(\omega)= r e^{i\phi}$$ with $$r, \phi$$ real
• $$r =$$ amplitude
• $$\phi =$$ phase
• This is the frequency domain

# The inverse Fourier transform

• The inverse Fourier transform recovers $$x(t)$$ from $$\tilde{x}(\omega)$$: $x(t) = \frac{1}{n}\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t \omega/n}}$
• We are back in the time domain
• Always comes out real
• The un-normalized inverse transform is $\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t \omega/n}}$

# An example

n <- 60
x <- sin(2 * pi * (1:60)/n) + 0.5 * sin(2 * pi * 3 * (1:60)/n)
plot(x, xlab = "t", ylab = expression(x(t)))

# In R

• You could do everything directly from the definitions, but this is slow
• Use the fast Fourier transform (FFT), from the fft function
• The trick is cute, but it’s just a computing trick
• The inverse=TRUE option gives the unnormalized inverse
• Because R starts counting from 1, not 0, give vector x of length n,
fft(x)[h] = sum_{k=1}^n {x[k]*exp(-2*pi*1i*(k-1)*(h-1)/n)}

so entry h in the output of fft contains the Fourier coefficient for frequency (h-1)/n

# In R

plot(x, pch = 16)

# In R

plot(x, pch = 16)
points(Re(fft(fft(x), inverse = TRUE)/length(x)), col = "red", pch = 25)

# Why do we care?

• Periodic components show up!
• Differentiation and integration are easy: $\frac{d}{dt}\left( c_{\omega} e^{2\pi i \omega t/n}\right) = \frac{2\pi i \omega}{n} c_{\omega} e^{2\pi i \omega t/n}$ so $\widetilde{x^{\prime}}(\omega) = \frac{2\pi i \omega}{n} \tilde{x}(\omega)$
• (similarly for integration)
• Shifts in time are easy: if $$y(t) = x(t-t_0)$$ then $\tilde{y}(\omega) = \tilde{x}(\omega) e^{-2\pi i t_0 \omega/n}$
• Convolutions are easy

# Convolutions

• Fix any function $$h$$ of time
• The convolution of $$x(t)$$ with $$h$$ is a new function of time: $\begin{eqnarray} (x \star h)(t) &\equiv & \sum_{s=0}^{n-1}{x(s) h(t-s)} \end{eqnarray}$
• A weighted sum of $$x$$ values around $$t$$, with weights given by $$h$$
• We’ve seen this sort of thing before!
• Moving averages = convolution with indicator function
• Kernel smoothing = convolution with the kernel

# Convolutions and Fourier transforms

• Convolution in the time domain is multiplication in the frequency domain: $\widetilde{x \star h}(\omega) = \tilde{x}(\omega) \tilde{h}(\omega)$
• Reverse is also true: convolution in the frequency domain is multiplication in the time domain
• Doing convolution quickly:
• Take the Fourier transforms of $$x$$ and $$h$$ (using fft)
• Multiply the Fourier coefficients
• Take the inverse transform
• Implemented in the convolve function

# Autocovariance and Fourier transforms

• Suppose $$X(t)$$ stationary, with autocovariance $$\gamma(h)$$
• We can do a Fourier transform on $$\gamma$$: $\tilde{\gamma}(\omega) = \sum_{h=0}^{n-1}{\gamma(h) e^{-2\pi i h \omega/n}}$
• How does this relate to $$\tilde{x}$$?

# Power spectrum

• $$X(t)$$ is random, so $$\tilde{X}(\omega)$$ is also random
• Amplitude $$|\tilde{X}(\omega)|$$ is random
• Power $$|\tilde{X}(\omega)|^2$$ is random
• Power spectrum $S(\omega) = \Expect{|\tilde{X}(\omega)|^2}$ is not random

# The Wiener-Khinchin theorem:

• The power spectrum is the Fourier transform of the autocovariance: $\tilde{\gamma}(\omega) = S(\omega) = \Expect{|\tilde{X}(\omega)|^2}$
• Uses:
• If we can estimate one side, we can estimate the other
• If we can estimate $$\Expect{|\tilde{X}(\omega)|^2}$$, we can estimate $$\gamma(h)$$

# The periodogram

• The plot of $$|\tilde{x}(\omega)|^2$$ against $$\omega$$ is the periodogram
• As $$n \rightarrow \infty$$, $$\Expect{|\tilde{x}(\omega)|^2} \rightarrow S(\omega)$$
• but $$\Var{|\tilde{x}(\omega)|^2} \not\rightarrow 0$$
• Solution: smooth the periodogram
• Basic R command: spectrum
• Some sensible defaults about smoothing
• annoyingly, different normalization of frequencies than fft

# In R

spectrum(x)