---
title: Fourier Analysis
date: 3 September 2026 (Lecture 4)
author: "36-740, Fall 2026"
output:
  slidy_presentation:
    math_method:
      engine: mathjax
      url: https://bactra.org/mathjax/tex-svg.js
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, warning=FALSE, message=FALSE,
               echo=FALSE)
```


\[
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]}
\newcommand{\TrueRegFunc}{\mu}
\newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}}
\newcommand{\TrueNoise}{\epsilon}
\newcommand{\EstNoise}{\widehat{\TrueNoise}}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator{\dof}{DoF}
\]


## Previously

- Decompose or analyze data into additive components:
    + Eigenvectors of the smoothing matrix (linear smoothers)
    + Eigenvectors of the variance matrix (principal compnents analysis)
- Today: decompose into _periodic_ additive components
    + Even when the data is _not_ periodic
	
## Start with the Karhunen-Loeve transform {.smaller}

- **Given**: $X(t)$ defined on a finite domain (e.g., interval) $D$, $\Expect{X(t)}=0$ (centering)
    + Substitute $\mathbf{r}$ for $t$, or $(\mathbf{r}, t)$, etc.
- **Desired** $X(t) = \sum_{i=1}^{\infty}{S_i \phi_i(t)}$ with random coefficients $S_i$
    + Implies $\Expect{S_i}=0$
- **Require** orthonormality of the functions: $\int_{D}{\phi_i(t)  \phi_j(t) dt} = \delta_{ij}$
- **Implies**  $\Cov{X(t), X(s)} = \sum_{i=1}^{\infty}{\sum_{j=1}^{\infty}}{\phi_i(t) \phi_j(s) \Cov{S_i, S_j}}$
- **Require** minimal mean squared error, so we pick the $\phi_i$ to minimize $\Expect{\int_{D}{\left(X(t) - \sum_{i=1}^{q}{S_i \phi_i(t)}\right)^2 dt}}$ at each $q$, subject to the orthonormality constraint
- **Implies**  $\int_{D}{\Cov{X(t), X(s)} \phi_i(s) ds} = \lambda_i \phi_i(t)$
    + An eigenproblem again, but now solved by an eigenfunction, not an eigenvector
	+ Eigenproblems for integral operators are classic early-20th-century math [@Courant-Hilbert-methods-of-math-phys, ch. III]
         * $\Cov{X(t), X(s)}$ is the **kernel** of the integral operator (and yes, stuff like this is [where our "kernel" comes from](https://bactra.org/weblog/two-kinds-of-kernel-methods.html))
- **Implies** $\Cov{S_i, S_j} = \lambda_i \delta_{ij}$
    + $\therefore$ $\Cov{X(t), X(s)} = \sum_{i}^{\infty}{\phi_i(t) \phi_i(s) \lambda_i}$
- We need some sort of estimate of the covariance function...

## Even more abstract Karhunen-Loeve

- **Given**: function $X$ in some space with an inner product $\langle X, Y \rangle$, with $\Expect{X} = 0$ (centering)
    + Associated norm $\|X\|^2 = \langle X, X \rangle$
- **Desired**: $X = \sum_{i=1}^{\infty}{S_i \phi_i}$ with random coefficients $S_i$
- **Require** orthonormality, $\langle \phi_i, \phi_j \rangle = \delta_{ij}$
- **Implies** $\Cov{X(t), X(s)} = \sum_{i}{\sum_{j}{\phi_i(t) \phi_j(s) \Cov{S_i, S_j}}}$
- **Require** minimal error, minimize $\Expect{\| X - \sum_{i=1}^{q}{S_i \phi_i}\|^2}$ over choice of $\phi_i$
- Define $\gamma_t(s) \equiv \Cov{X(t), X(s)}$
- Doing the minimization, with the Lagrange multiplier to enforce orthonormality, gives us $\langle \gamma_t, \phi_i \rangle = \lambda_i \phi_i(t)$ which is an eigenvalue equation
    + $\langle \gamma_t, af + bg \rangle = a \langle \gamma_t, f \rangle + b \langle \gamma_t, g \rangle$ so this is a linear operator
- And we get the implication that $\Cov{S_i, S_j} = \lambda_i \delta_{ij}$

## Some things to notice about Karhunen-Loeve / PCA {.smaller}

- We are transforming to new coordinates
- The expected value of the transform is zero
    + Because we centered $X$; otherwise, add the constant function to soak up the mean
- The variance of the transformed coordinates is _not_ zero
    + We're dealing with a random process
- The variance of transformed coordinates tells us a _lot_ about the underlying process
    + $\Var{S_i} = \lambda_i$ and $\lambda_i$ is an eigenvalue
- The covariance across the transformed coordinates is zero
- The transformation is different for every process

## Let's use sine waves! {.smaller}

- Let's _fix_ the basis functions, independent of the random process
- There are many, many choices
- Sine waves are particularly use for _stationary_ processes
- **Stationary** = the distribution does not change over time
- **Weakly stationary** = **second-order stationary** = mean, variance and covariance do not change over time
- This plays nicely with sine waves because $e^{2\pi i \omega (t-t_0)} = e^{-2\pi i \omega t_0} e^{2\pi i \omega t}$
   + Shifting in time changes the phase of the wave, in a trivial way, but nothing else
   + Most fixed choices of basis function are not so nice!
- "Decompose an arbitrary function into a linear combination of sine waves" = **Fourier analysis**
- Pros: the math is nice, and well-developed
- Cons: 
    + We'll need a _lot_ of coefficients / transformed coordinates
    + The new coefficients won't be uncorrelated (in general)
    + We need complex-valued coefficients

## A few notational conventions for complex-valued vectors

- $c = a+bi$ $\Leftrightarrow$ $\overline{c} = a - bi$ (complex conjugation)
- $|c|^2 = \overline{c}c = c\overline{c} =$ modulus squared of $c$
- $\overline{\mathbf{c}} =$ complex conjugate each element of $\mathbf{c}$
- $\mathbf{c}^{\dagger} = \overline{\mathbf{c}^T} =$ transpose and complex conjugate
   + Read as "$\mathbf{c}$ adjoint"
   + Doesn't matter what order you transpose and conjugate
- For a complex-valued vector $\vec{c}$, $\|\vec{c}\|^2 = \vec{c}^{\dagger} \vec{c}$
    + think of $\vec{c}$ as a $p\times 1$ matrix

## The Fourier transform (in 1D) {.smaller}

- For continuous time,
\[
\tilde{x}(\omega) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}{x(t) e^{2\pi i \omega t} dt}
\]
    + Note: a continuous infinity of basis functions
- For discrete time, $j=1, \ldots n$,
\[
\tilde{x}_k = \frac{1}{\sqrt{n}}\sum_{j=1}^{n}{x_j e^{2\pi i k j / n}}
\]
    + Only $n$ (orthonormal) basis functions
    + If the physical time between measurements ($\equiv$ **sampling interval**) is $\Delta t$, then the corresponding physical frequencies $\omega_k = \frac{k}{n \Delta t}$
- The Fourier transform is a linear operator say $\mathcal{F}$
    + Maps functions on the **time domain** into functions on the **frequency domain**
        * Strictly speaking $\mathcal{F}$ maps complex-valued functions of the real line to complex-valued functions of the real line, $\mathcal{F}: (\mathbf{C}^{\mathbf{R}}) \mapsto (\mathbf{C}^{\mathbf{R}})$, but domain and the range mean different things


## The Fourier transform (in 1D) {.smaller}

- The **inverse Fourier transform** is
\[
\tilde{\tilde{x}}(t) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}{\tilde{x}(\omega) e^{-2\pi i \omega t} d\omega} = x(t)
\]
or
\[
\tilde{\tilde{x}}_j = \frac{1}{\sqrt{n}}\sum_{k=-n/2}^{n/2}{\tilde{x}_k e^{-2\pi i k j/n}} = x_j
\]
    + This is, again, a linear operator
    + The inverse transform takes frequency-domain functions to time-domain functions
    + The inverse transform looks _almost_ just like the forward transform
    + WARNING: Some people use slightly different conventions for forward/inverse FT; some people absorb the $2\pi$ into $\omega$ or $k$



## Convolution {.smaller}

- The **convolution** of $f$ and $g$ is a new function, $f*g$:
\[
(f*g)(t) \equiv \int{f(s) g(t-s) ds} = (g*f)(t)
\]
    + **EXERCISE**: Prove $f*g = g*f$
    + similarly in discrete time
- If $g$ is a pdf, then $(f*g)$ is a weighted average of $f$ at each point
- Convolving $f$ with a pdf gives a new, smoother function
    + Kernel (= Nadaraya-Watson) smoothing is convolving with a pdf (= the kernel)
- If $\mathrm{pdf}(X) = f$ and $\mathrm{pdf}(Z) = g$ and $X$ and $Z$ are independent, then $Y=X+Z$ has pdf $f*g$
- The Fourier transform turns convolution into multiplication:
\[
\widetilde{f*g}(\omega) = \tilde{f}(\omega) \tilde{g}(\omega)
\]
   + Conversely, it turns multiplication into convolution: $\widetilde{fg} = \tilde{f} * \tilde{g}$
    + Proving this is left as an **EXERCISE**
        * See the reading if you get stuck, but try it on your own


## Spectrum

- The Fourier coefficients are called the **spectrum** of the function
- Scalar multiplication multiplies the spectrum
    + The Fourier transform is a linear operator
- Time shifts change the spectrum _only_ by phase, not in magnitude
- The Fourier coefficients are random variables
- The Fourier coefficients _remain_ random even in the infinite-time limit
   + Because they need to reproduce that particular random trajectory exactly
   
## Some examples

- A (Dirac) delta function: $\delta(t)$ is defined by
\[
\int{f(t) \delta(t) dt} = f(0)
\]
    + Therefore $\int{f(t) \delta(t-t_0) dt} = f(t_0)$
- The Fourier transform of the delta function is
\[
\tilde{\delta}(\omega) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}{e^{2\pi i \omega t} \delta(t) dt} = \frac{e^{2\pi i \omega 0}}{\sqrt{2\pi}} = \frac{1}{\sqrt{2\pi}}
\]
    + Equal Fourier coefficients at all frequencies
- The transform of a Fourier basis function of frequency $\omega_0$ is a delta function at $\omega_0$ (**exercise**)
- The Fourier transform of a Gaussian is another Gaussian (**exercise**)


## The power spectrum

- The **power spectrum** is the modulus squared:
\[
S(\omega)= |\tilde{X}(\omega)|^2
\]
    + The power (= energy per unit time) transmitted by a physical wave is $\propto$ amplitude squared
        * At least for small amplitudes...
- The power spectrum does not care about the phase of the Fourier coefficients, just their magnitude
- $\therefore$ the power spectrum is invariant under translation

## Wiener-Khinchin

- Suppose $X(t)$ is second-order stationary, so that
\[
\Cov{X(t), X(s)} = \gamma(|t-s|)
\]
the **autocovariance function**
- $\gamma$ has some Fourier transform
- **Wiener-Khinchin theorem**:
\[
\tilde{\gamma}(\omega) = \Expect{S(\omega)}
\]

## Estimating the power spectrum {.smaller}

#### The **periodogram**

- Take the discrete Fourier transform of $X_j$, get $\tilde{X}_k$
- Set $\hat{S}(\omega_k) = |\tilde{X}_{k}|^2$
- This is unbiased
- And $\Cov{\hat{S}(\omega), \hat{S}(\nu)} \rightarrow 0$ as $n\rightarrow \infty$
- _But_ we have limited frequency resolution $\Rightarrow$ "leakage" of power from nearby frequencies
- _And_ $\Var{\hat{S}(\omega)} \not\rightarrow 0$ as $n\rightarrow \infty$

#### Practically

- Smooth the periodogram by convolving it with a pdf
    + _Deliberately_ add bias, to reduce variance
- This is actually where statisticians first used kernel methods, in the 1950s
    + The people who did this went on to invent kernel density estimation [@Rosenblatt-introduces-kernel-density-estimation; @Parzen-introduces-kernel-density-estimation]
- Smoothing less and less as we get more data gives consistent estimates of the power spectrum



## Power spectrum interpretation

- If there's a periodic component with period $\tau$, there should be a peak in the power spectrum at frequency $1/\tau$
    + Also at multiples $2/\tau$, $3/\tau$, etc. (**harmonics**)
- There _can_ be peaks in the power spectrum for other reasons as well




## Covariance $\Leftrightarrow$ Power spectrum

- The autocovariance function and the power spectrum are a Fourier-transform pair
- $\therefore$ either one can be used to find the other
- $\therefore$ either one puts constraints on the other
- Every (valid) covariance function is the Fourier transform of a (valid) power spectrum
- Less obvious: for every (valid) power spectrum, there is some stationary process with that autocovariance function
    + Useful for simulation
- Practically: estimate whichever one is easier to get, and then Fourier transform to get the other


## Spatial and spatio-temporal Fourier transforms

- Replace the scalar frequency $\omega$ with a vector of appropriate dimension
   + Physicists calls the resulting basis functions "plane waves"
      * In 3D, the surfaces of constant phase are planes, perpendicular to the direction the wave is traveling in
- Otherwise, everything works very much as before

## Filtering, smoothing, trends {.smaller}

- The inverse Fourier transform recovers the original time series _exactly_
- We can deliberately change some of the Fourier coefficients, and then the inverse transform will give us something else
- One common trick: $\tilde{x}(\omega) \rightarrow \tilde{x}(\omega) \kappa(\omega)$ where $\kappa(u) \geq 0$, $\kappa(u) \rightarrow 0$ as $|u|\rightarrow \infty$
    + $\kappa$ is sometimes, but not always, a pdf
- Then $\hat{x} = \mathcal{F}^{-1} (\tilde{x} \kappa)$
- The fact that $\kappa(u) \rightarrow 0$ as $|u| \rightarrow \infty$ makes this a **low-pass filter**
    + Low-frequency components are preserved, high-frequency components are attenuated
    + What would a **high-pass** filter be?
- Multiplication in the frequency domain is convolution in the time domain: $\hat{x} = x * \tilde{\kappa}$
    + This is a smoother again
- One way to estimate a trend



## Non-Stationarity {.smaller}

1. We can de-trend, and then apply Fourier analysis to the fluctuations.
2. We can break the data up into segments, and hope it's (approximately) stationary within each segment.
    a. A **spectrogram** is a plot showing (local) Fourier transforms for successive time segments
    b. This leads to a test of stationarity: we should get the same Fourier spectrum from each segment, plus or minus sampling error.
       i. @Fuentes-Guttorp-Sampson-transforms, sec. 1.6.2 gives one version of this idea, with conditions under which we can make it into an ANOVA test and use the standard machinery
       ii. @Kantz-Schreiber-2nd, sec. 2.2--2.3 and 13.2, give what is (essentially) a bootstrap version of the same idea
3. Give up on the Fourier basis in favor of **wavelets**
    - Fourier basis functions have equal amplitude everywhere
    - Come up with wiggly, wave-like functions that are localized
	- It's just another set of basis functions, so we just take inner products again


## Summing up

- We can always decompose our data into an additive combination of sine waves (=Fourier basis functions)
- This lets us express many things like convolution very simply
- Because sine waves just change phase under translation, this is extra useful for stationary processes
- The autocovariance function is the Fourier transform of the power spectrum


## Backup: Fourier transforms and calculus

- Inverse Fourier transform: $x(t) = \frac{1}{\sqrt{2\pi}}\int{\tilde{x}(\omega) e^{-2\pi i \omega t} d\omega}$
- Say $\dot{x}(t) = \frac{dx}{dt}(t)$
- Exchanging integral and derivative, $\dot{x}(t) = \frac{1}{2\pi}\int{\tilde{x}(\omega) \left( \frac{d}{dt} e^{-2\pi i \omega t} \right) dt}$
- Therefore $\widetilde{\dot{x}}(\omega) = - 2\pi i \omega \tilde{x}(\omega)$
    + Exercise: Find the parallel result for the time integral of $x$
- Similarly for higher derivatives of $x(t)$: the Fourier transform of the $k^{\mathrm{th}}$ derivative is $(-2\pi i \omega)^k \tilde{x}(\omega)$
    + Some people use this to _define_ fractional derivatives and fractional integrals
- Differentiating is one form of high-pass filtering; integrating is a form of low-pass filtering


## Backup: Measure Theory

- Something like the Dirac delta "function" actually makes a lot more sense as a measure
    + Specifically $\delta$ is the probability measure that puts probability 1 on the origin.
- Every measure $\mu$ defines a linear map from functions to the real or complex numbers, via $\int{f(x) \mu(dx)}$
    + Hence the physicists' folklore about how $\delta$ only makes sense in integrals
    + Hence also the **de Finetti notation**, abbreviating $\int{f(x) \mu(dx)} = \mu f$.
    + Question: is every linear map from functions to numbers an integral w.r.t. some measure?
- In the more general theory of Fourier analysis, the spectrum is really a measure $\tilde{F}$, and we should write $x(t) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}{e^{-2\pi \omega t} \tilde{F}(d\omega)}$
    + The integrals I wrote are for the well-behaved case where $\tilde{F}$ is absolutely continuous w.r.t. Lebesgue measure, with density $\tilde{x}(\omega)$
- There are good accounts of all this in lots of books on stochastic processes, e.g. @Loeve-probability, @Bartlett-stochastic, etc.
    + @CRS-almost-none, ch. 19, is my attempt to make this as comprehensible as possible
	
## Backup: Simulating from a power spectrum {.smaller}

> **Theorem**: Suppose $g(\omega) \geq 0$, and define $G(\omega) = \int_{-\infty}^{\omega}{g(\nu) d\nu}$.  Further suppose that $G(-\infty) = 0$ and $G(\infty)=\sigma^2/\sqrt{2\pi} > 0$.  Then there exists a centered, weakly-stationary stochastic process $X$ with autocovariance function $\tilde{g}$.  In particular, $\Var{X(t)} = \sigma^2$ for all $t$, and $\Cov{X(t_1), X(t_2)} = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}{ e^{-2\pi i (t_2-t_1) \omega} g(\omega) d\omega}$.

**Sketch proof**: Draw $W \sim G/(\sigma^2 \sqrt{2\pi})$.  That is, the pdf of $W$ is $g(\omega)/(\sigma^2 \sqrt{2 \pi})$, and $\Pr{\left(W leq \omega\right)} = G(\omega)/(\sigma^2 \sqrt{2\pi})$.  Draw $A \sim \mathrm{Unif}(0, 2\pi)$, independent of $N$.  Set $X(t) \equiv \frac{\sigma} e^{iA} e^{-iNt}$.  It's easy to check that $\Expect{e^{iA}} = 0$, so $\Expect{X(t)} = 0$ too.  And
\[
\Expect{\overline{X(t_1)} X(t_2)} = \sigma^2 \Expect{e^{iN(t_2 -t_1)}} = \int_{-\infty}^{\infty}{e^{-2\pi i \omega (t_2-t_1)} \frac{g(\omega)}{\sqrt{2\pi}} d\omega}
\]
as desired.

Notes:

- This is one of those places where things would be cleaner if I didn't have the factor of $1/\sqrt{2\pi}$ in the inverse Fourier transform
- We need to complex-conjugate either $X(t_1)$ or $X(t_2)$ to get a covariance function for complex-valued processes.
- You can modify this to handle dis-continuous power spectra, too.


## Backup: Simulating from a power spectrum (cont'd.) {.smaller}

- The construction on the previous slide might not be very satisfying, because every realization is just a sine wave, and the stochasticity only becomes apparent across realizations
- Here's a related idea: for each $\omega$, $A(\omega) \mathrm{Unif}(0, 2\pi)$, IIDly, except that $A(-\omega) = \overline{A(\omega)}$
\[
X(t) =  (2\pi)^{-1/4} \int{e^{iA(\omega)} e^{-2\pi i \omega t} \sqrt{g(\omega)} d\omega}
\]
Now
\[
\Expect{\overline{X(t_1)} X(t_2)} = \frac{1}{\sqrt{2\pi}} \int{\left( \int{\Expect{e^{-iA(\omega)} e^{iA(\nu)}} e^{-2\pi i (\omega t_i - \nu t_2)} \sqrt{g(\nu)} d\nu} \right) \sqrt{g(\omega)} d\omega} = \frac{1}{\sqrt{2\pi}} \int{e^{-2\pi i \omega (t_2 - t_1)} g(\omega) d\omega}
\]
(How did I get rid of the double integral?)









## References