---
title: Fourier Methods I
date: 2 October 2018
author: 36-467/36-667
output: slidy_presentation
bibliography: locusts.bib
---

```{r, include=FALSE}
# General set-up options
library(knitr)
opts_chunk$set(size="small", background="white", highlight=FALSE,
               cache=TRUE, autodep=TRUE,
               tidy=TRUE,  tidy.opts=list(comment=FALSE),
               warning=FALSE, message=FALSE,
               echo=FALSE)
```

\[
\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$

```{r}
n <- 60
```


```{r, echo=TRUE}
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

```{r, echo=TRUE}
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

```{r, echo=TRUE}
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", add=TRUE)
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 {.smaller}

\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 {.smaller}

\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 {.smaller}

\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

```{r, echo=TRUE}
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)))
```

## An example

```{r}
plot(fft(x), xlab="Real part of Fourier coefficient",
     ylab="Imaginary part of Fourier coefficient")
```

## An example

```{r}
plot(Mod(fft(x)), xlab=expression(omega),
         ylab="Magnitude of Fourier coefficient",
         xlim=c(1,30))
```


## 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

```{r, echo=TRUE}
plot(x, pch=16)
```

## In R


```{r, echo=TRUE}
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

```{r, echo=TRUE}
spectrum(x)
```

## In R

```{r, echo=TRUE}
library(gstat); data(wind)
wind <- wind[!(wind$month==2 & wind$day==29),]
spectrum(wind[,"DUB"], span=10)
```

## In R

```{r, echo=TRUE}
library(gstat); data(wind)
wind <- wind[!(wind$month==2 & wind$day==29),]
spectrum(wind[,"DUB"], span=10)
abline(v=1/365)
```


## Summary

## Some cautions:

- Lots of variants in notation/definitions
    + some people absorb the $2\pi$ into the definition of $\omega$
	+ some people absorb the $1/n$ into the definition of $\omega$
	+ some people write $k$ or $h$ or $f$ instead of $\omega$
	+ some people normalize the basis functions

## Backup: Eliminating sine+cosine in favor of sine

\begin{eqnarray}
a_{\omega} \sin{2\pi \omega t/n} + b_{\omega} \cos{2\pi \omega t/n}
& = & \sqrt{a_\omega^2+b_\omega^2}\left(\frac{a_{\omega}}{\sqrt{a_\omega^2+b\omega^2}}\sin{2\pi \omega t/n} +  \frac{b_{\omega}}{\sqrt{a_\omega^2+b\omega^2}}\cos{2\pi \omega t/n} \right) ~ \text{(algebra)}\\
& = & \sqrt{a_\omega^2+b_\omega^2}\left(\cos{\phi_\omega} \sin{2\pi \omega t/n} + \sin{\phi_\omega} \cos{2\pi \omega t/n}\right) ~ \text{(definition of sine and cosine)}\\
& = & \sqrt{a_\omega^2+b\omega^2} \sin{(2\pi \omega t/n + \phi_\omega)} ~ \text{(angle-adding trig identity)}\\
\end{eqnarray}

## Backup: why $x(t) e^{-2\pi i \omega t/n}$?

The hand-wavy explanation:

- When we take the inner product of two real vectors $\vec{u}$ and $\vec{v}$, it's the same as $\vec{u}^T \vec{v}$ ($[1\times p]$ $[p\times 1]$)
- But the inner product of anything with itself is a (squared) length, so $\geq 0$
- If $\vec{v} = a+bi$, then $\vec{v}^T \vec{v} = (a+bi)(a+bi) = a^2 - b^2 + 2iab$, which doesn't make sense as a _length_
- Solution: don't just transpose one vector, transpose it _and_ take its complex conjugate
    + Sometimes called finding the **adjoint** vector
	+ Often written $\vec{v}^*$ or $\vec{v}^\dagger$
- Now $\vec{v}^* \vec{v} = (a-bi)(a+bi) = a^2 + b^2$, which makes sense as a length
- We want to take the inner product of the vector of $x(t)$'s with the basis vector we get from $e^{2\pi i \omega t/n}$, so we need to complex-conjugate the basis vector

## Backup: why $x(t) e^{-2\pi i \omega t/n}$?

The "shut up and calculate" explanation:

- The (complex) numbers $\tilde{x}(\omega)$ are _defined_ to be the ones which give us back the
(real) numbers $x(t)$:
\[
x(t) = \frac{1}{n}\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t \omega/n}}
\]
- Pick your favorite frequency $\nu$, and multiply both sides by $e^{-2\pi i t \nu/n}$:
\[
x(t) e^{-2\pi i t \nu/n}  = \frac{1}{n}\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t (\omega-\nu)/n}}
\]
- Sum over all $t$:
\[
\sum_{t=1}^{n}{x(t) e^{-2\pi i t \nu/n}} = \frac{1}{n}\sum_{t=1}^{n}{\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega) e^{2\pi i t (\omega-\nu)/n}}}
\]
- Exchange the order of summation on the right-hand side:
\[
\sum_{t=1}^{n}{x(t) e^{-2\pi i t \nu/n}} = \frac{1}{n}\sum_{\omega=0}^{n/2-1}{\tilde{x}(\omega)\sum_{t=1}^{n}{ e^{2\pi i t (\omega-\nu)/n}}}
\]
- Now remember that sines and cosines average to zero, _except_ at zero frequency:
\[
\sum_{t=1}^{n}{x(t) e^{-2\pi i t \nu/n}} = \frac{1}{n}\tilde{x}(\nu)n = \tilde{x}(\nu)
\]