# Fast Fourier Transform: From Theory to Frequency Discovery
Related notes:
- [[MATH310S26-Day7-Notes]] - Introduction to Fourier series
- [[MATH310S26-Day10-Notes]] - Complex Fourier series
- [[MATH310S26-Day11-Notes]] - Fourier regression theory
- [[MATH310S26-Day13-WorkdayMaterials (Fourier regression and seasonal decomposition)]] - Pre-specified frequency approach
- [[MATH310S26-Day14-Notes]] - Delta functions and Fourier transforms
- [[MATH310S26-Day15-Notes]] - Fourier transform examples
## Recent Lectures Recap
- Developed Fourier series for periodic functions on $[-L, L]$
- Extended to Fourier transforms by taking $L \to \infty$ limit
- Connected continuous and discrete via delta function sampling
- Note: Just as uniform discrete time sampling in frequency led to periodicity in time, the reverse is also true, i.e., uniform sampling in time leads to periodicity in frequency.
- In [[MATH310S26-Day13-WorkdayMaterials (Fourier regression and seasonal decomposition)|Day 13]], we pre-specified frequencies based on physical reasoning
## Day 16 Goals
- Understand the [Fast Fourier Transform (FFT)](https://en.wikipedia.org/wiki/Fast_Fourier_transform) as a computational tool
- Use FFT to **discover** dominant frequencies in GMSL data
- Compare data-driven frequency selection with our Day 13 physical intuition
- Visualize and interpret [power spectra](https://en.wikipedia.org/wiki/Spectral_density)
- Reconstruct signals using FFT-discovered frequencies
- Preview: How frequency analysis connects to filtering (upcoming topic)
# From Fourier Series to FFT: The Mathematical Journey
For a more detailed overview, see: [[FourierTheoryAndFFT.pdf]]
## The Conceptual Path: FS → FT → Sampling → FFT
### 1. [Fourier Series](https://en.wikipedia.org/wiki/Fourier_series) (Periodic Functions)
From [[MATH310S26-Day10-Notes|Day 10]], for a periodic function with period $2L$:
$f(t) = \sum_{n=-\infty}^{\infty} c_n e^{i\omega_n t}, \quad \omega_n = \frac{n\pi}{L}$
The coefficients are:
$c_n = \frac{1}{2L} \int_{-L}^{L} f(t) e^{-i\omega_n t} dt$
**Key insight**: [Discrete frequencies](https://en.wikipedia.org/wiki/Discrete_frequency_domain) $\omega_n$ because of finite period.
### 2. [Fourier Transform](https://en.wikipedia.org/wiki/Fourier_transform) (Non-periodic Functions)
From [[MATH310S26-Day12-Notes|Day 12]] and [[MATH310S26-Day14-Notes|Day 14]], taking $L \to \infty$:
**Forward Transform**:
$\hat{f}(\omega) = \mathcal{F}\left\{f\right\}=\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} f(t) e^{-i\omega t} dt$
**Inverse Transform**:
$f(t)= \mathcal{F}^{-1}\left\{\hat{f}\right\} = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \hat{f}(\omega) e^{i\omega t} d\omega$
**Key insight**: [Continuous frequency spectrum](https://en.wikipedia.org/wiki/Frequency_spectrum) $\omega \in \mathbb{R}$ for non-periodic signals.
### 3. The Bridge: Delta Functions and Sampling
From [[MATH310S26-Day14-Notes|Day 14]], the [Dirac delta function](https://en.wikipedia.org/wiki/Dirac_delta_function) connects discrete and continuous:
**For sinusoids** (which don't have classical Fourier transforms):
$\mathcal{F}\{\cos(\omega_0 t)\} = \sqrt{\frac{\pi}{2}}[\delta(\omega - \omega_0) + \delta(\omega + \omega_0)]$
The delta functions "pick out" specific frequencies, just like Fourier series!
**For sampled data**: Real measurements happen at discrete times $t_k = k\Delta t$:
$f_{\text{sampled}}(t) = f(t) \cdot \sum_{k=-\infty}^{\infty} \delta(t - k\Delta t)$
This sampling creates periodicity in the frequency domain (we'll see this as [aliasing](https://en.wikipedia.org/wiki/Aliasing)).
### 4. The FFT: A Computational Tool
For finite samples $\{f_0, f_1, \ldots, f_{N-1}\}$, the **Fast Fourier Transform**:
- Efficiently computes frequency content
- Uses clever factorization ([Cooley-Tukey algorithm](https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm))
- The Most Important Algorithm Of All Time, Veritasium, Nov 3, 2022, has an interesting historical perspective. Jump to 10:20 to see the origin story!
- Reduces computation from $O(N^2)$ to $O(N \log N)$
- For our GMSL data (N≈1200): ~100× faster!
**What FFT gives us**:
- Complex numbers at discrete frequencies
- Magnitude = strength of oscillation
- Phase = timing/shift of oscillation
- We'll interpret these physically for our sea level data
## Physical Interpretation for Our Data
For GMSL data sampled monthly over ~31 years:
- **[Sampling rate](https://en.wikipedia.org/wiki/Sampling_(signal_processing))**: 12 samples/year → $f_s = 12$ year⁻¹
- **[Nyquist frequency](https://en.wikipedia.org/wiki/Nyquist_frequency)**: $f_{Nyquist} = f_s/2 = 6$ cycles/year
- **Frequency resolution**: $\Delta f = 1/(N\Delta t) \approx 1/31$ cycles/year
- **Interpretable range**: 0 to 6 cycles/year (beyond that is aliased)
## What the FFT Gives Us
The FFT output $F_n$ contains:
1. **Magnitude**: $|F_n|$ - How strong is frequency $n$?
2. **Phase**: $\arg(F_n)$ - What's the phase shift?
3. **[Power](https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density)**: $|F_n|^2$ - Energy at frequency $n$
4. **[Power Spectral Density](https://en.wikipedia.org/wiki/Spectral_density)**: $|F_n|^2/(N \cdot \Delta f)$ - Power per unit frequency
# Applying FFT to Global Mean Sea Level Data
We'll use the same GMSL data from Day 13, but now let the FFT tell us which frequencies matter!
* GMSL data from 1993-2024: [NASA Sea Level Portal](https://sealevel.nasa.gov/), or [global_mean_sea_level_1993-2024.csv (dropbox)](https://www.dropbox.com/scl/fi/9uqs42k9k3tgqwl45gdn9/global_mean_sea_level_1993-2024.csv?rlkey=d62q93mjrso72ql76dhzpolg4&dl=0), or Canvas for these data
## Python Implementation
Step 1. Load data and compute time parameters
```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Load data
df = pd.read_csv("global_mean_sea_level_1993-2024.csv")
df = df.dropna(subset=['GMSLWithGIA'])
t = df['YearPlusFraction'].values
S = df['GMSLWithGIA'].values
# Time parameters
dt = np.mean(np.diff(t)) # Average time step (years)
fs = 1/dt # Sampling frequency (samples/year)
N = len(S)
print(f"Data: {N} points, dt = {dt:.4f} years, fs = {fs:.2f} samples/year")
```
Purpose: Load GMSL data and establish time sampling parameters for FFT.
Step 2. Compute FFT and one-sided PSD
```python
# Compute FFT
S_fft = np.fft.fft(S)
# Create frequency array (one-sided)
freqs = np.fft.fftfreq(N, dt)[:N//2] # Positive frequencies only
# Compute one-sided PSD
# Factor of 2 accounts for negative frequencies (except [DC component](https://en.wikipedia.org/wiki/DC_bias))
psd = np.abs(S_fft[:N//2])**2 / N
psd[1:] *= 2 # Double non-DC components
psd /= fs # Normalize by sampling frequency
print(f"Frequency range: 0 to {freqs[-1]:.2f} cycles/year")
print(f"Frequency resolution: {freqs[1]:.4f} cycles/year")
print(f"Number of frequency bins: {len(freqs)}")
```
Purpose: Transform to frequency domain and compute power spectral density.
> **Note:** The FFT does **not** produce a continuous spectrum. It returns exactly $N/2$ discrete numbers — one per frequency bin. The "smooth curve" you will see in plots is just the plotting library connecting these discrete points with straight lines. There is no signal information *between* them. For this dataset, $N = 1168$, so you have $N/2 - 1 = 583$ non-DC frequency bins, spaced $\Delta f = 1/(N \cdot \Delta t) \approx 0.024$ cycles/year apart.
Step 3. Visualize PSD in both linear and log (dB) scales
```python
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
# Linear scale
ax1.plot(freqs, psd, 'b-', linewidth=0.8)
ax1.set_xlabel('Frequency (cycles/year)')
ax1.set_ylabel('PSD (mm²·year)')
ax1.set_title('Power Spectral Density - Linear Scale')
ax1.set_xlim(0, 5)
ax1.grid(True, alpha=0.3)
# Decibel scale (10*log10 for power)
psd_dB = 10 * np.log10(psd / np.max(psd)) # Normalize to max = 0 dB
ax2.plot(freqs, psd_dB, 'b-', linewidth=0.8)
ax2.set_xlabel('Frequency (cycles/year)')
ax2.set_ylabel('PSD (dB relative to peak)')
ax2.set_title('Power Spectral Density - Logarithmic Scale')
ax2.set_xlim(0, 5)
ax2.set_ylim(-60, 5)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
```
Purpose: Linear scale shows dominant frequencies clearly. [dB scale](https://en.wikipedia.org/wiki/Decibel) reveals weaker signals hidden by strong ones. The dB value is defined as $10 \log_{10}(\text{PSD}/\text{PSD}_{\max})$, so the strongest peak is always at 0 dB and everything else is negative. This is identical to how audio engineers represent the [dynamic range](https://en.wikipedia.org/wiki/Dynamic_range) of a recording — it compresses a factor-of-10,000 variation down to a readable 40 dB span.
Step 4. Identify prominent frequencies using threshold
```python
# Use threshold: peaks > 10× median power (excluding DC)
psd_no_dc = psd[1:] # Exclude DC component
threshold = 10 * np.median(psd_no_dc)
# Find frequencies above threshold
prominent_mask = psd > threshold
prominent_freqs = freqs[prominent_mask]
prominent_psd = psd[prominent_mask]
# Sort by power
sort_idx = np.argsort(prominent_psd)[::-1]
prominent_freqs = prominent_freqs[sort_idx]
prominent_psd = prominent_psd[sort_idx]
print(f"Threshold value: {threshold:.1f} mm²·year")
print(f"Number of bins above threshold: {np.sum(prominent_mask)}")
print("Prominent frequencies (cycles/year):")
for f, p in zip(prominent_freqs, prominent_psd):
if f > 0:
period = 1/f
print(f" f = {f:.3f} (period = {period:.2f} years), Power = {p:.1f}")
```
Purpose: Use a data-driven threshold to flag frequency bins with unusually high power.
> **Note:** This threshold is a **value cutoff**, not a rank-based filter. The code keeps every bin whose PSD exceeds $10 \times \text{median}$, regardless of how many that turns out to be — for our data this is roughly 68 bins. This is fundamentally different from selecting the "top 68 bins by power." In a rank-based approach, you choose a count first and let the data determine the cutoff value. Here you choose a cutoff value first and let the data determine the count. Critically, entire **high-power regions** of the spectrum will be flagged, including the slopes and saddles *between* peaks. A bin sitting in a "valley" between two large peaks can still be far above the global median. We will diagnose this problem in Step 9.
Step 5. Visualize prominent frequencies
```python
plt.figure(figsize=(10, 5))
plt.semilogy(freqs, psd, 'b-', linewidth=0.8, alpha=0.5, label='Full spectrum')
plt.semilogy(prominent_freqs[prominent_freqs > 0],
prominent_psd[prominent_freqs > 0],
'ro', markersize=8, label='Above threshold')
# Annotate expected frequencies
plt.axvline(1.0, color='orange', linestyle='--', alpha=0.5, label='Annual')
plt.axvline(2.0, color='green', linestyle='--', alpha=0.5, label='Semi-annual')
# Add threshold line
plt.axhline(threshold, color='red', linestyle=':', alpha=0.5,
label=f'10×median threshold = {threshold:.1f}')
plt.xlabel('Frequency (cycles/year)')
plt.ylabel('PSD (mm²·year, log scale)')
plt.title('Prominent Frequencies in GMSL (Day 16: FFT Discovery)')
plt.xlim(0, 3)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
```
Purpose: Visual confirmation of the threshold selection.
> **Note:** You will likely see red dots appearing at what look like valleys or downward slopes in the spectrum. This is **expected and correct** — those are FFT bins that sit within a broad high-power region but happen to be lower than their immediate neighbors. They are above the global threshold even though they are not local maxima. Think of it like flood water: once the water level (the threshold) is high enough, it submerges entire mountain ranges, not just the peaks — including the passes and saddles between summits. This is the fundamental limitation of the threshold approach, and it is precisely what Steps 9–10 are designed to fix.
Step 6. Fourier regression on discovered frequencies
```python
# Select non-DC prominent frequencies
regression_freqs = prominent_freqs[(prominent_freqs > 0.1) & (prominent_freqs < 10)]
print(f"Regression frequencies: {regression_freqs} cycles/year")
# Build design matrix
X = np.ones((N, 1 + 2*len(regression_freqs)))
for i, freq in enumerate(regression_freqs):
X[:, 1 + 2*i] = np.cos(2*np.pi*freq*t)
X[:, 2 + 2*i] = np.sin(2*np.pi*freq*t)
# Solve normal equations: beta = (X^T X)^{-1} X^T S
beta = np.linalg.solve(X.T @ X, X.T @ S)
S_seasonal = X @ beta
# Extract just the oscillatory part (no mean)
X_no_mean = X[:, 1:]
beta_no_mean = beta[1:]
S_oscillations = X_no_mean @ beta_no_mean
print(f"Mean sea level: {beta[0]:.1f} mm")
```
Purpose: Use FFT-discovered frequencies for regression, separating mean from oscillations. Each frequency $f_k$ contributes two columns to the [design matrix](https://en.wikipedia.org/wiki/Design_matrix) — $\cos(2\pi f_k t)$ and $\sin(2\pi f_k t)$ — which together can represent an oscillation of any [amplitude](https://en.wikipedia.org/wiki/Amplitude) and [phase](https://en.wikipedia.org/wiki/Phase_(waves)) at that frequency.
Step 7. Remove seasonal and fit quadratic trend
```python
# Remove seasonal oscillations
S_deseasonalized = S - S_oscillations
# Center time for numerical stability (but interpret carefully)
t_mean = np.mean(t)
t_centered = t - t_mean
# Fit quadratic trend: S = a + b*t_c + c*t_c^2
A = np.column_stack([np.ones(N), t_centered, t_centered**2])
coeffs = np.linalg.solve(A.T @ A, A.T @ S_deseasonalized)
S_trend = A @ coeffs
# Extract coefficients
a, b, c = coeffs
print(f"Centered quadratic fit: S = {a:.1f} + {b:.2f}*t_c + {c:.4f}*t_c²")
print(f"Acceleration: {2*c:.4f} mm/year²")
# For centered model: rate = b + 2*c*(t - t_mean)
rate_1993 = b + 2*c*(t[0] - t_mean)
rate_2024 = b + 2*c*(t[-1] - t_mean)
print(f"Rate at 1993: {rate_1993:.2f} mm/year")
print(f"Rate at 2024: {rate_2024:.2f} mm/year")
```
Purpose: Quadratic trend captures the acceleration of sea level rise. The [instantaneous rate](https://en.wikipedia.org/wiki/Derivative) of rise is the derivative of the quadratic: $\frac{dS}{dt} = b + 2c \cdot t_c$. Because $c > 0$, the rate increases over time — sea level is not just rising, it is accelerating.
Step 8. Compare with Intuitive Fourier Regression (Day 13)
```python
# Day 13 - Intuitive Fourier Regression: Pre-specified k=1,2 with quadratic trend
X_intuitive = np.column_stack([np.ones(N), t_centered, t_centered**2,
np.cos(2*np.pi*1*t), np.sin(2*np.pi*1*t),
np.cos(2*np.pi*2*t), np.sin(2*np.pi*2*t)])
beta_intuitive = np.linalg.solve(X_intuitive.T @ X_intuitive, X_intuitive.T @ S)
S_intuitive = X_intuitive @ beta_intuitive
# Day 16 - Fourier Regression on Prominent Frequencies: FFT-discovered with quadratic trend
S_prominent = S_trend + S_oscillations
# Compute R²
SS_tot = np.sum((S - np.mean(S))**2)
SS_res_intuitive = np.sum((S - S_intuitive)**2)
SS_res_prominent = np.sum((S - S_prominent)**2)
R2_intuitive = 1 - SS_res_intuitive/SS_tot
R2_prominent = 1 - SS_res_prominent/SS_tot
print(f"\nModel Comparison:")
print(f"Day 13 - Intuitive Fourier Regression (2 freqs): R² = {R2_intuitive:.4f}")
print(f"Day 16 - Prominent Frequencies ({len(regression_freqs)} freqs): R² = {R2_prominent:.4f}")
print(f"Improvement: ΔR² = {R2_prominent - R2_intuitive:.4f}")
# Visualize first 5 years
idx = t < t[0] + 5
plt.figure(figsize=(10, 5))
plt.plot(t[idx], S[idx], 'ko', markersize=3, alpha=0.5, label='Data')
plt.plot(t[idx], S_intuitive[idx], 'b-', linewidth=1.5,
label=f'OG Fourier Regression (2 freqs, R²={R2_intuitive:.3f})')
plt.plot(t[idx], S_prominent[idx], 'r-', linewidth=1.5,
label=f'Prominent Frequencies ({len(regression_freqs)} freqs, R²={R2_prominent:.3f})')
plt.xlabel('Year')
plt.ylabel('Sea Level (mm)')
plt.title('Model Comparison: OG Fourier vs FFT-discovered Frequencies')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
```
Purpose: Compare FFT-discovered frequencies with our physics-based Day 13 approach.
Step 9. The problem with threshold approach
```python
# Our threshold approach found many frequencies
print(f"\nThreshold approach found {len(regression_freqs)} 'significant' frequencies")
print("But are they all real peaks? Let's investigate...")
# Plot to see the issue
plt.figure(figsize=(10, 5))
plt.semilogy(freqs, psd, 'b-', linewidth=0.8, alpha=0.5)
plt.semilogy(prominent_freqs[prominent_freqs > 0],
prominent_psd[prominent_freqs > 0],
'ro', markersize=4, alpha=0.5)
plt.axhline(threshold, color='red', linestyle=':',
label=f'Threshold = {threshold:.1f}')
plt.xlabel('Frequency (cycles/year)')
plt.ylabel('PSD (log scale)')
plt.title('Problem: Threshold selects high-power regions, not peaks!')
plt.text(2.5, threshold*2, f'{len(regression_freqs)} freqs\nabove threshold!',
color='red', fontsize=9)
plt.xlim(0, 3)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
print("Notice: Many 'significant' points are just part of broad features, not peaks!")
```
Purpose: Visualize why the threshold approach picks up too many frequencies.
> **Note:** Look carefully at the subtitle — it says something like "Found 28 peaks among 68 points above threshold." Those are two different numbers. There are 68 **bins** above the threshold, but only 28 of those are actual local maxima. The other 40 lie on slopes and saddles within broad high-power spectral regions. They were flagged not because they are peaks but simply because the entire neighborhood around the annual and semi-annual cycles is energetic enough to clear the global threshold. This is the central issue: **a value cutoff cannot distinguish a peak from a shoulder**.
Step 10. Smart peak detection approach
```python
from scipy.signal import find_peaks
# Exclude very low frequencies (those are trend, not oscillations)
mask = freqs > 0.3
psd_for_peaks = psd.copy()
psd_for_peaks[~mask] = 0
# Find peaks with prominence (how much they stand out above local background)
peaks, properties = find_peaks(psd_for_peaks,
prominence=20, # Must stand out from neighbors
distance=int(0.2/freqs[1])) # Min separation
print(f"\nPeak detection found {len(peaks)} peaks!")
peak_freqs = freqs[peaks]
peak_psd = psd[peaks]
for f, p in zip(peak_freqs, peak_psd):
print(f" f = {f:.3f} (period = {1/f:.2f} years), Power = {p:.0f}")
# [Signal-to-noise ratio](https://en.wikipedia.org/wiki/Signal-to-noise_ratio)
window_size = int(0.5 / freqs[1])
background = np.zeros_like(psd)
for i in range(len(psd)):
start = max(0, i - window_size//2)
end = min(len(psd), i + window_size//2)
background[i] = np.percentile(psd[start:end], 25)
peak_ratio = psd / (background + 1e-10)
# Identify which detected peaks to actually use
retained_mask = np.zeros(len(peak_freqs), dtype=bool)
for i, (f, ratio_val) in enumerate(zip(peak_freqs, peak_ratio[peaks])):
if f > 0.5 and f < 4.0 and ratio_val > 5:
retained_mask[i] = True
retained_freqs = peak_freqs[retained_mask]
retained_psd = peak_psd[retained_mask]
notused_freqs = peak_freqs[~retained_mask]
notused_psd = peak_psd[~retained_mask]
print(f"\nUsing {np.sum(retained_mask)} peaks out of {len(peak_freqs)} detected")
print("Peaks to use in regression:", retained_freqs)
print("Detected but not used:", notused_freqs)
# Visualize
plt.figure(figsize=(12, 5))
plt.semilogy(freqs, psd, 'b-', linewidth=0.8, alpha=0.7)
threshold_mask = (psd > threshold) & (freqs > 0.5)
plt.semilogy(freqs[threshold_mask], psd[threshold_mask], 'o',
color='red', markersize=4, alpha=0.2,
label=f'Above threshold ({np.sum(threshold_mask)} points)')
trend_idx = np.where((freqs > 0) & (freqs < 0.3))[0]
if len(trend_idx) > 0:
max_trend_idx = trend_idx[np.argmax(psd[trend_idx])]
plt.semilogy(freqs[max_trend_idx], psd[max_trend_idx], 's',
color='grey', markersize=10,
label=f'Trend (~{1/freqs[max_trend_idx]:.0f} year "period")')
if len(notused_freqs) > 0:
plt.semilogy(notused_freqs, notused_psd, 'o', markersize=8,
color='green', alpha=0.3,
label=f'Detected, not used (n={len(notused_freqs)})')
if len(retained_freqs) > 0:
plt.semilogy(retained_freqs, retained_psd, 'o', markersize=8,
color='green', alpha=1.0,
label=f'Used in model (n={len(retained_freqs)})')
plt.axvline(1.0, color='orange', linestyle='--', alpha=0.5)
plt.axvline(2.0, color='darkgreen', linestyle='--', alpha=0.5)
plt.xlim(0, 3)
plt.xlabel('Frequency (cycles/year)')
plt.ylabel('PSD (log scale)')
plt.title('Peak Detection Process: From Threshold to Smart Selection')
plt.legend(loc='upper right', fontsize=9)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
```
Purpose: Proper peak detection requires a bin to be a **[local maximum](https://en.wikipedia.org/wiki/Maximum_and_minimum)** (greater than both neighbors) and to stand out significantly above its local background. This is fundamentally different from the global threshold — a bin can be globally high but locally unremarkable (a shoulder), and a bin can be globally modest but locally prominent (a sharp peak against a quiet background).
Step 11. Model comparison with detected peaks
```python
# Helper function for model fitting
def fit_and_evaluate(freq_list, label):
X = np.ones((N, 1 + 2*len(freq_list)))
for i, freq in enumerate(freq_list):
X[:, 1 + 2*i] = np.cos(2*np.pi*freq*t)
X[:, 2 + 2*i] = np.sin(2*np.pi*freq*t)
beta = np.linalg.solve(X.T @ X, X.T @ S)
S_fourier = X @ beta
S_oscillations = S_fourier - beta[0]
S_detrended = S - S_oscillations
A = np.column_stack([np.ones(N), t_centered, t_centered**2])
coeffs = np.linalg.solve(A.T @ A, A.T @ S_detrended)
S_trend = A @ coeffs
S_full = S_trend + S_oscillations
SS_tot = np.sum((S - np.mean(S))**2)
SS_res = np.sum((S - S_full)**2)
R2 = 1 - SS_res/SS_tot
print(f"{label:30s} | R² = {R2:.4f}")
return S_full, R2
print("Model | R²")
print("-" * 45)
if len(regression_freqs) > 20:
S_thresh, R2_thresh = fit_and_evaluate(regression_freqs[:20],
f"Thresholded ({len(regression_freqs)} freqs, cap 20)")
else:
S_thresh, R2_thresh = fit_and_evaluate(regression_freqs,
f"Thresholded ({len(regression_freqs)} freqs)")
S_allpeaks, R2_allpeaks = fit_and_evaluate(peak_freqs,
f"All Peaks ({len(peak_freqs)} freqs)")
# Snap to nearest FFT bin (not exactly 1.0 and 2.0 — the grid may not align)
intuitive_freqs = [freqs[np.argmin(np.abs(freqs - 1.0))],
freqs[np.argmin(np.abs(freqs - 2.0))]]
S_intuitive_fit, R2_intuitive_fit = fit_and_evaluate(intuitive_freqs,
"OG Fourier Regression (2 freqs)")
print("\nConclusion: Peak detection prevents overfitting!")
print("The OG Fourier Regression (intuitive approach) remains optimal.")
```
Purpose: Compare models using [$R^2$](https://en.wikipedia.org/wiki/Coefficient_of_determination) to evaluate how well each frequency selection fits the data.
> **Note:** When building the OG Fourier Regression we use `freqs[np.argmin(np.abs(freqs - 1.0))]` rather than simply `1.0`. This is because the FFT frequency grid does not necessarily include exactly 1.000 or 2.000 cycles/year — the bins are spaced $\Delta f \approx 0.024$ cycles/year apart and may land at 0.997 or 1.021. Snapping to the nearest available bin ensures our regression frequencies are consistent with the FFT output. This is related to the concept of [frequency resolution](https://en.wikipedia.org/wiki/Spectral_resolution).
---
## R Implementation
Step 1. Load data and compute time parameters
```r
library(ggplot2)
library(gridExtra)
if (exists("rstudioapi") && rstudioapi::isAvailable()) {
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
}
# Load data
df <- read.csv("global_mean_sea_level_1993-2024.csv")
df <- df[!is.na(df$GMSLWithGIA), ]
t <- df$YearPlusFraction
S <- df$GMSLWithGIA
# Time parameters
dt <- mean(diff(t))
fs <- 1/dt
N <- length(S)
cat(sprintf("Data: %d points, dt = %.4f years, fs = %.2f samples/year\n",
N, dt, fs))
```
Purpose: Load GMSL data and establish time sampling parameters for FFT.
Step 2. Compute FFT and one-sided PSD
```r
# Compute FFT
S_fft <- fft(S)
# Create frequency array (one-sided): k / (N * dt) for k = 0, 1, ..., N/2-1
freqs <- (0:(N/2-1)) / (N * dt)
# Compute one-sided PSD
psd <- Mod(S_fft[1:(N/2)])^2 / N
psd[2:length(psd)] <- 2 * psd[2:length(psd)] # Double non-DC components
psd <- psd / fs # Normalize by sampling frequency
cat(sprintf("Frequency range: 0 to %.2f cycles/year\n", tail(freqs, 1)))
cat(sprintf("Frequency resolution: %.4f cycles/year\n", freqs[2]))
cat(sprintf("Number of frequency bins: %d\n", length(freqs)))
```
Purpose: Transform to frequency domain and compute power spectral density.
> **Note:** The FFT does **not** produce a continuous spectrum. It returns exactly $N/2$ discrete numbers — one per frequency bin. The "smooth curve" you will see in plots is ggplot connecting these discrete points; there is no signal information *between* them. For this dataset, $N = 1168$, so $N/2 - 1 = 583$ non-DC frequency bins, spaced $\Delta f = 1/(N \cdot \Delta t) \approx 0.024$ cycles/year apart. The frequency array `freqs` is simply the index $k$ divided by the total time span: $f_k = k/(N \cdot \Delta t)$.
Step 3. Visualize PSD in both linear and log (dB) scales
```r
psd_df <- data.frame(freq = freqs, psd = psd)
# Linear scale plot
p1 <- ggplot(psd_df[psd_df$freq <= 5,], aes(x=freq, y=psd)) +
geom_line(color="blue", linewidth=0.8) +
labs(x="Frequency (cycles/year)", y="PSD (mm²·year)",
title="Power Spectral Density - Linear Scale") +
theme_minimal() +
theme(panel.grid.minor=element_blank())
# dB scale: 10*log10(psd/max) so that the dominant peak = 0 dB
psd_dB <- 10 * log10(psd / max(psd))
psd_df$psd_dB <- psd_dB
p2 <- ggplot(psd_df[psd_df$freq <= 5,], aes(x=freq, y=psd_dB)) +
geom_line(color="blue", linewidth=0.8) +
labs(x="Frequency (cycles/year)", y="PSD (dB relative to peak)",
title="Power Spectral Density - Logarithmic Scale") +
ylim(-60, 5) +
theme_minimal() +
theme(panel.grid.minor=element_blank())
psd_plot <- grid.arrange(p1, p2, ncol=2)
ggsave("r_fig0_psd_scales.png", psd_plot, width=12, height=4, dpi=150)
```
Purpose: Linear scale shows dominant frequencies clearly. dB scale reveals weaker signals hidden by strong ones. The dB value is $10 \log_{10}(\text{PSD}/\text{PSD}_{\max})$, so the strongest peak sits at 0 dB and everything else is negative. This compresses a factor-of-10,000 power variation into a readable 40 dB span — the same reason audio engineers use dB to describe dynamic range.
Step 4. Identify prominent frequencies using threshold
```r
# Use threshold: bins > 10× median power (excluding DC)
psd_no_dc <- psd[-1]
threshold <- 10 * median(psd_no_dc)
# Find ALL bins above threshold (not just peaks)
prominent_idx <- which(psd > threshold)
prominent_freqs <- freqs[prominent_idx]
prominent_psd <- psd[prominent_idx]
# Sort by power
sort_idx <- order(prominent_psd, decreasing=TRUE)
prominent_freqs <- prominent_freqs[sort_idx]
prominent_psd <- prominent_psd[sort_idx]
cat(sprintf("Threshold value: %.1f mm²·year\n", threshold))
cat(sprintf("Number of bins above threshold: %d\n", length(prominent_idx)))
cat("Prominent frequencies (cycles/year):\n")
for (i in 1:length(prominent_freqs)) {
f <- prominent_freqs[i]
p <- prominent_psd[i]
if (f > 0) {
period <- 1/f
cat(sprintf(" f = %.3f (period = %.2f years), Power = %.1f\n",
f, period, p))
}
}
```
Purpose: Use a data-driven threshold to flag frequency bins with unusually high power.
> **Note:** This threshold is a **value cutoff**, not a rank-based filter. The code flags every bin whose PSD exceeds $10 \times \text{median}$, regardless of how many that turns out to be. For our data this is roughly 68 bins — but that count is a consequence of the data, not a deliberate choice. Compare this to selecting the "top 68 bins by power": in that case you choose the count first and the threshold falls where it falls. Here you choose the threshold first and accept however many bins clear it. Crucially, entire **high-power spectral regions** will be flagged, including the slopes and saddles between actual peaks. A bin sitting in a "valley" between two large peaks can still be far above the global median. We will diagnose this in Step 9.
Step 5. Visualize prominent frequencies
```r
plot_df <- data.frame(freq = freqs, psd = psd)
prominent_df <- data.frame(
freq = prominent_freqs[prominent_freqs > 0],
psd = prominent_psd[prominent_freqs > 0]
)
p <- ggplot(plot_df[plot_df$freq <= 3,], aes(x=freq, y=psd)) +
geom_line(color="blue", linewidth=0.8, alpha=0.5) +
geom_point(data=prominent_df[prominent_df$freq <= 3,],
aes(x=freq, y=psd), color="red", size=3) +
geom_vline(xintercept=1.0, color="orange", linetype="dashed", alpha=0.5) +
geom_vline(xintercept=2.0, color="green", linetype="dashed", alpha=0.5) +
geom_hline(yintercept=threshold, color="red", linetype="dotted", alpha=0.5) +
scale_y_log10() +
labs(x="Frequency (cycles/year)", y="PSD (mm²·year, log scale)",
title="Prominent Frequencies in GMSL (Day 16: Fourier Regression on Prominent Frequencies)",
subtitle=sprintf("Red dotted line = 10×median threshold = %.1f", threshold)) +
theme_minimal() +
annotate("text", x=1.0, y=max(psd)/2, label="Annual",
color="orange", angle=90, vjust=-0.5) +
annotate("text", x=2.0, y=max(psd)/2, label="Semi-annual",
color="green", angle=90, vjust=-0.5)
print(p)
ggsave("r_fig1_prominent_frequencies.png", p, width=10, height=5, dpi=150)
```
Purpose: Visual confirmation of the threshold selection.
> **Note:** You will likely see red dots appearing at what look like valleys or downward slopes in the spectrum. This is **expected and correct**. Remember: the FFT spectrum is a set of discrete points connected by lines — there is no actual data between them. A red dot in an apparent "valley" is simply an FFT bin that (a) sits within a broad high-power region so it clears the global threshold, but (b) happens to be lower than its two immediate neighbors. The plotting library draws a connecting line between neighboring bins, making it look like a valley, but the bin itself is above the threshold. This is the problem Steps 9–10 are designed to solve.
Step 6. Fourier regression on discovered frequencies
```r
# Select non-DC prominent frequencies for regression
regression_freqs <- prominent_freqs[(prominent_freqs > 0.1) &
(prominent_freqs < 10)]
cat(sprintf("Using %d regression frequencies\n", length(regression_freqs)))
cat("Regression frequencies:", round(regression_freqs, 3), "cycles/year\n")
# Build design matrix: intercept + cos/sin pairs for each frequency
X <- matrix(1, nrow=N, ncol=1 + 2*length(regression_freqs))
for (i in 1:length(regression_freqs)) {
freq <- regression_freqs[i]
X[, 1 + 2*(i-1) + 1] <- cos(2*pi*freq*t)
X[, 1 + 2*(i-1) + 2] <- sin(2*pi*freq*t)
}
# Solve normal equations: beta = (X^T X)^{-1} X^T S
beta <- solve(t(X) %*% X) %*% (t(X) %*% S)
S_seasonal <- X %*% beta
# Extract just the oscillatory part (no mean)
X_no_mean <- X[, -1, drop=FALSE]
beta_no_mean <- beta[-1]
S_oscillations <- X_no_mean %*% beta_no_mean
cat(sprintf("Mean sea level: %.1f mm\n", beta[1]))
```
Purpose: Use FFT-discovered frequencies for regression, separating mean from oscillations. Each frequency $f_k$ contributes two columns — $\cos(2\pi f_k t)$ and $\sin(2\pi f_k t)$ — which together represent an oscillation of any amplitude and phase at that frequency.
Step 7. Remove seasonal and fit quadratic trend
```r
# Remove seasonal oscillations from signal
S_deseasonalized <- S - S_oscillations
# Center time for numerical stability
t_mean <- mean(t)
t_centered <- t - t_mean
# Fit quadratic trend: S = a + b*t_c + c*t_c^2
A <- cbind(1, t_centered, t_centered^2)
coeffs <- solve(t(A) %*% A) %*% (t(A) %*% S_deseasonalized)
S_trend <- A %*% coeffs
a <- coeffs[1]
b <- coeffs[2]
c <- coeffs[3]
cat(sprintf("Centered quadratic fit: S = %.1f + %.2f*t_c + %.4f*t_c²\n", a, b, c))
cat(sprintf("Acceleration: %.4f mm/year²\n", 2*c))
# Instantaneous rate = b + 2*c*(t - t_mean)
rate_1993 <- b + 2*c*(t[1] - t_mean)
rate_2024 <- b + 2*c*(tail(t, 1) - t_mean)
rate_center <- b
cat(sprintf("Rate at 1993: %.2f mm/year\n", rate_1993))
cat(sprintf("Rate at 2024: %.2f mm/year\n", rate_2024))
cat(sprintf("Rate at center (%.1f): %.2f mm/year\n", t_mean, rate_center))
```
Purpose: [Quadratic trend](https://en.wikipedia.org/wiki/Polynomial_regression) captures the acceleration of sea level rise. The instantaneous rate of rise is $dS/dt = b + 2c \cdot t_c$. Because $c > 0$, the rate increases over time. The coefficient $c$ itself is half the acceleration: acceleration $= 2c$ mm/year².
Step 8. Compare with Intuitive Fourier Regression (Day 13)
```r
# Day 13 - Intuitive Fourier Regression: Pre-specified k=1,2 with quadratic trend
X_intuitive <- cbind(1, t_centered, t_centered^2,
cos(2*pi*1*t), sin(2*pi*1*t),
cos(2*pi*2*t), sin(2*pi*2*t))
beta_intuitive <- solve(t(X_intuitive) %*% X_intuitive) %*% (t(X_intuitive) %*% S)
S_intuitive <- X_intuitive %*% beta_intuitive
# Day 16 - FFT-discovered: quadratic trend + oscillations
S_prominent <- S_trend + S_oscillations
# Compute R²
SS_tot <- sum((S - mean(S))^2)
SS_res_intuitive <- sum((S - S_intuitive)^2)
SS_res_prominent <- sum((S - S_prominent)^2)
R2_intuitive <- 1 - SS_res_intuitive/SS_tot
R2_prominent <- 1 - SS_res_prominent/SS_tot
cat("\nModel Comparison:\n")
cat(sprintf("Day 13 - OG Fourier Regression (2 freqs): R² = %.4f\n", R2_intuitive))
cat(sprintf("Day 16 - Prominent Frequencies (%d freqs): R² = %.4f\n",
length(regression_freqs), R2_prominent))
cat(sprintf("Improvement: ΔR² = %.4f\n", R2_prominent - R2_intuitive))
# Visualize first 5 years
idx <- which(t < t[1] + 5)
# Pre-compute label strings so factor levels and color keys match exactly
label_og <- sprintf("OG Fourier Regression (2 freqs, R2=%.3f)", R2_intuitive)
label_prom <- sprintf("Prominent Frequencies (%d freqs, R2=%.3f)",
length(regression_freqs), R2_prominent)
comparison_df <- data.frame(
t = rep(t[idx], 3),
S = c(S[idx], S_intuitive[idx], S_prominent[idx]),
Model = factor(rep(c("Data", label_og, label_prom), each=length(idx)),
levels=c("Data", label_og, label_prom))
)
p <- ggplot(comparison_df, aes(x=t, y=S, color=Model)) +
geom_point(data=comparison_df[comparison_df$Model=="Data",],
size=1, alpha=0.5) +
geom_line(data=comparison_df[comparison_df$Model!="Data",],
linewidth=1.2, alpha=0.8) +
scale_color_manual(values=setNames(c("black", "blue", "red"),
c("Data", label_og, label_prom))) +
labs(x="Year", y="Sea Level (mm)",
title="Model Comparison: OG Fourier vs FFT-discovered Frequencies") +
theme_minimal() +
theme(legend.position="top")
print(p)
ggsave("r_fig_step8_model_comparison.png", p, width=10, height=5, dpi=150)
```
Purpose: Compare FFT-discovered frequencies with our physics-based Day 13 approach.
Step 9. The problem with threshold approach
```r
cat(sprintf("\nThreshold approach found %d bins above threshold (not all are peaks!)\n",
length(prominent_freqs[prominent_freqs > 0.1])))
cat("But are they all real peaks? Let's investigate...\n")
psd_plot_df <- data.frame(freq = freqs, psd = psd)
prominent_df <- data.frame(freq = prominent_freqs[prominent_freqs > 0],
psd = prominent_psd[prominent_freqs > 0])
p <- ggplot(psd_plot_df[psd_plot_df$freq <= 3,], aes(x=freq, y=psd)) +
geom_line(color="blue", alpha=0.5, linewidth=0.8) +
geom_point(data=prominent_df[prominent_df$freq <= 3,],
color="red", size=2, alpha=0.5) +
geom_hline(yintercept=threshold, color="red", linetype="dashed", alpha=0.5) +
scale_y_log10() +
labs(x="Frequency (cycles/year)", y="PSD (log scale)",
title="Problem: Threshold selects high-power regions, not peaks!",
subtitle=sprintf("Red dashed line = threshold = %.1f; Found %d peaks among %d points above threshold",
threshold, length(regression_freqs), length(prominent_idx))) +
theme_minimal()
print(p)
ggsave("r_fig2_threshold_problem.png", p, width=10, height=5, dpi=150)
cat("Notice: Many 'significant' points are just part of broad features, not peaks!\n")
```
Purpose: Visualize why the threshold approach picks up too many frequencies.
> **Note:** The subtitle of this plot tells the story: there are many **bins** above the threshold but only a fraction are actual local maxima. The rest are on slopes and broad shoulders around the annual and semi-annual peaks. A value cutoff cannot distinguish the summit from the hillside — it just tells you everything above a certain elevation.
Step 10. Smart peak detection approach
```r
# Focus on oscillations — zero out the low-frequency trend region
mask <- freqs > 0.3
psd_for_peaks <- psd
psd_for_peaks[!mask] <- 0
# Enhanced peak detection: require local maximum AND prominence above local background
find_peaks <- function(x, min_distance=1, prominence_threshold=20) {
peaks <- c()
for (i in 2:(length(x)-1)) {
if (x[i] > x[i-1] & x[i] > x[i+1]) {
if (length(peaks) == 0 || (i - tail(peaks, 1)) >= min_distance) {
# Prominence check: peak must stand out from local neighborhood
window <- 10
start <- max(1, i - window)
end <- min(length(x), i + window)
local_background <- median(x[start:end])
prominence <- x[i] - local_background
if (prominence > prominence_threshold) {
peaks <- c(peaks, i)
}
}
}
}
return(peaks)
}
min_distance <- round(0.2 / freqs[2]) # 0.2 cycles/year minimum separation
peaks <- find_peaks(psd_for_peaks, min_distance, prominence_threshold=20)
cat(sprintf("\nFound %d distinct peaks:\n", length(peaks)))
cat("Freq (cy/yr) | Period (yr) | Power\n")
cat("-------------|-------------|-------\n")
peak_freqs <- c()
peak_powers <- c()
peak_indices <- c()
for (peak_idx in peaks) {
f <- freqs[peak_idx]
p <- psd[peak_idx]
peak_freqs <- c(peak_freqs, f)
peak_powers <- c(peak_powers, p)
peak_indices <- c(peak_indices, peak_idx)
cat(sprintf("%12.3f | %11.2f | %6.0f\n", f, 1/f, p))
}
# [Signal-to-noise ratio](https://en.wikipedia.org/wiki/Signal-to-noise_ratio) (kept for reference output)
window_size <- round(0.5 / freqs[2])
background <- numeric(length(psd))
for (i in 1:length(psd)) {
start <- max(1, i - window_size/2)
end <- min(length(psd), i + window_size/2)
background[i] <- quantile(psd[start:end], 0.25)
}
peak_ratio <- psd / (background + 1e-10)
cat(sprintf("\nUsing all %d detected peaks in regression\n", length(peak_freqs)))
# Find the dominant low-frequency peak (trend, not oscillation)
low_freq_idx <- which(freqs > 0 & freqs < 0.3)
if (length(low_freq_idx) > 0) {
max_low_freq_idx <- low_freq_idx[which.max(psd[low_freq_idx])]
trend_point_df <- data.frame(
freq = freqs[max_low_freq_idx],
psd = psd[max_low_freq_idx]
)
}
# Build visualization
threshold_mask <- (psd > threshold) & (freqs > 0.5)
threshold_df <- data.frame(freq=freqs[threshold_mask], psd=psd[threshold_mask])
# All detected oscillatory peaks
peak_df <- data.frame(freq=peak_freqs, psd=peak_powers)
psd_df <- data.frame(freq=freqs, psd=psd, ratio=peak_ratio)
p <- ggplot(psd_df[psd_df$freq <= 3,], aes(x=freq, y=psd)) +
geom_line(color="blue", alpha=0.7, linewidth=0.8)
if (nrow(threshold_df) > 0)
p <- p + geom_point(data=threshold_df[threshold_df$freq <= 3,],
aes(x=freq, y=psd), color="red", alpha=0.2, size=2)
if (nrow(peak_df) > 0) {
p <- p + geom_point(data=peak_df[peak_df$freq <= 3,],
aes(x=freq, y=psd), color="green", alpha=1.0, size=4)
for (i in 1:nrow(peak_df)) {
f <- peak_df$freq[i]; pd <- peak_df$psd[i]
if (abs(f - 1.0) < 0.1)
p <- p + annotate("text", x=f, y=pd*2, label="Annual", hjust=0.5, size=3)
else if (abs(f - 2.0) < 0.1)
p <- p + annotate("text", x=f, y=pd*2, label="Semi-annual", hjust=0.5, size=3)
}
}
p <- p +
geom_vline(xintercept=1.0, color="orange", linetype="dashed", alpha=0.5) +
geom_vline(xintercept=2.0, color="darkgreen", linetype="dashed", alpha=0.5) +
scale_y_log10() +
labs(x="Frequency (cycles/year)", y="PSD (log scale)",
title="Peak Detection Process: From Threshold to Smart Selection") +
theme_minimal() +
theme(panel.grid.minor=element_blank()) +
annotate("text", x=2.5, y=max(psd)*0.9, label="● Above threshold",
color="red", alpha=0.2, hjust=0, size=3) +
annotate("text", x=2.5, y=max(psd)*0.6, label="■ Trend (polynomial regression)",
color="grey", hjust=0, size=3) +
annotate("text", x=2.5, y=max(psd)*0.4, label="● Detected peaks (used in model)",
color="green", hjust=0, size=3)
# Grey box for trend peak — added last for visibility
if (exists("trend_point_df") && nrow(trend_point_df) > 0)
p <- p + geom_point(data=trend_point_df, aes(x=freq, y=psd),
color="grey", shape=15, size=3)
print(p)
ggsave("r_fig3_peak_detection.png", p, width=12, height=5, dpi=150)
```
Purpose: Proper peak detection requires a bin to be a local maximum **and** stand above its local background — this filters out the slopes and shoulders that the global threshold included, leaving only genuine spectral peaks to use in regression.
Step 11. Model comparison with peak-based selection
```r
fit_with_freqs <- function(freqs_to_use, t, S, label) {
N <- length(t)
X <- matrix(1, nrow=N, ncol=1 + 2*length(freqs_to_use))
for (i in seq_along(freqs_to_use)) {
freq <- freqs_to_use[i]
X[, 1 + 2*(i-1) + 1] <- cos(2*pi*freq*t)
X[, 1 + 2*(i-1) + 2] <- sin(2*pi*freq*t)
}
beta <- solve(t(X) %*% X) %*% (t(X) %*% S)
S_fourier <- X %*% beta
S_oscillations <- S_fourier - beta[1]
S_detrended <- S - S_oscillations
t_centered <- t - mean(t)
X_trend <- cbind(1, t_centered, t_centered^2)
beta_trend <- solve(t(X_trend) %*% X_trend) %*% (t(X_trend) %*% S_detrended)
S_trend <- X_trend %*% beta_trend
S_full <- S_trend + S_oscillations
SS_tot <- sum((S - mean(S))^2)
SS_res <- sum((S - S_full)^2)
R2 <- 1 - SS_res/SS_tot
return(list(R2=R2, S_fit=S_full))
}
cat("\nModel Selection Results:\n")
cat("Model | Frequencies | R²\n")
cat("-----------------------|-------------|--------\n")
# Threshold approach
nfreqs_thresh <- min(length(regression_freqs), 20)
result_thresh <- fit_with_freqs(regression_freqs[1:nfreqs_thresh], t, S, "Threshold")
cat(sprintf("Threshold Approach | %d (cap 20) | %.4f\n",
length(regression_freqs), result_thresh$R2))
# All detected peaks
if (length(peak_freqs) > 0) {
result_peaks <- fit_with_freqs(peak_freqs, t, S, "Detected Peaks")
cat(sprintf("Detected Peaks | %d | %.4f\n",
length(peak_freqs), result_peaks$R2))
}
# OG Fourier Regression: snap to nearest FFT bin
closest_1 <- freqs[which.min(abs(freqs - 1.0))]
closest_2 <- freqs[which.min(abs(freqs - 2.0))]
result_og <- fit_with_freqs(c(closest_1, closest_2), t, S, "OG Fourier")
cat(sprintf("OG Fourier Regression | 2 | %.4f\n",
result_og$R2))
cat("\nConclusion: The OG Fourier Regression (annual + semi-annual) remains optimal.\n")
cat("FFT validates our physical intuition from Day 13!\n")
# Visualization
idx_5yr <- t < t[1] + 5
plot_df_5yr <- data.frame(
t = rep(t[idx_5yr], 4),
S = c(S[idx_5yr],
result_thresh$S_fit[idx_5yr],
if(exists("result_peaks")) result_peaks$S_fit[idx_5yr] else S[idx_5yr],
result_og$S_fit[idx_5yr]),
Model = factor(rep(c("Data", "Threshold Approach", "Detected Peaks",
"OG Fourier Regression"), each=sum(idx_5yr)))
)
p1 <- ggplot(plot_df_5yr, aes(x=t, y=S, color=Model)) +
geom_point(data=plot_df_5yr[plot_df_5yr$Model=="Data",], size=1, alpha=0.5) +
geom_line(data=plot_df_5yr[plot_df_5yr$Model!="Data",], linewidth=1.2, alpha=0.8) +
scale_color_manual(values=c("Data"="black", "Threshold Approach"="red",
"Detected Peaks"="green",
"OG Fourier Regression"="blue")) +
labs(x="", y="Sea Level (mm)", title="Model Comparison: First 5 Years Detail") +
theme_minimal() + theme(legend.position="top", legend.text=element_text(size=8))
plot_df_full <- data.frame(
t = rep(t, 4),
S = c(S, result_thresh$S_fit,
if(exists("result_peaks")) result_peaks$S_fit else S,
result_og$S_fit),
Model = factor(rep(c("Data", "Threshold Approach", "Detected Peaks",
"OG Fourier Regression"), each=N))
)
p2 <- ggplot(plot_df_full, aes(x=t, y=S, color=Model)) +
geom_line(data=plot_df_full[plot_df_full$Model=="Data",],
linewidth=0.3, alpha=0.3) +
geom_line(data=plot_df_full[plot_df_full$Model!="Data",],
linewidth=0.8, alpha=0.6) +
scale_color_manual(values=c("Data"="black", "Threshold Approach"="red",
"Detected Peaks"="green",
"OG Fourier Regression"="blue")) +
labs(x="Year", y="Sea Level (mm)", title="Model Comparison: Full Time Series") +
theme_minimal() + theme(legend.position="top", legend.text=element_text(size=8))
combined_plot <- grid.arrange(p1, p2, ncol=1, heights=c(1, 1))
ggsave("r_fig4_model_comparison.png", combined_plot, width=12, height=8, dpi=150)
```
Purpose: Compare models using [$R^2$](https://en.wikipedia.org/wiki/Coefficient_of_determination) to evaluate how well each frequency selection fits the data.
> **Note:** When building the OG model, we use `which.min(abs(freqs - 1.0))` to find the FFT bin *nearest* to 1.0 cycles/year rather than assuming the grid lands there exactly. The FFT frequency bins are spaced $\Delta f \approx 0.024$ cycles/year apart and may not include exactly 1.000 or 2.000. Snapping to the nearest bin keeps the model consistent with the actual FFT output.
---
## MATLAB/Octave Implementation
Step 1. Load data and compute time parameters
```matlab
% Load data
data = readtable('global_mean_sea_level_1993-2024.csv');
valid = ~isnan(data.GMSLWithGIA);
t = data.YearPlusFraction(valid);
S = data.GMSLWithGIA(valid);
% Time parameters
dt = mean(diff(t));
fs = 1/dt;
N = length(S);
fprintf('Data: %d points, dt = %.4f years, fs = %.2f samples/year\n', ...
N, dt, fs);
```
Purpose: Load GMSL data and establish time sampling parameters for FFT.
Step 2. Compute FFT and one-sided PSD
```matlab
% Compute FFT
S_fft = fft(S);
% Create frequency array (one-sided): k / (N*dt) for k = 0, 1, ..., N/2-1
freqs = (0:N/2-1) / (N * dt);
% Compute one-sided PSD
psd = abs(S_fft(1:N/2)).^2 / N;
psd(2:end) = 2 * psd(2:end); % Double non-DC components
psd = psd / fs; % Normalize by sampling frequency
fprintf('Frequency range: 0 to %.2f cycles/year\n', freqs(end));
fprintf('Frequency resolution: %.4f cycles/year\n', freqs(2));
fprintf('Number of frequency bins: %d\n', length(freqs));
```
Purpose: Transform to frequency domain and compute power spectral density.
> **Note:** The FFT does **not** produce a continuous spectrum. It returns exactly $N/2$ discrete numbers — one per frequency bin. The smooth curve you see in plots is MATLAB connecting these discrete points with straight lines; there is no signal information *between* them. For this dataset, $N = 1168$, giving $N/2 - 1 = 583$ non-DC frequency bins spaced $\Delta f = 1/(N \cdot \Delta t) \approx 0.024$ cycles/year apart. The frequency array is constructed as $f_k = k / (N \cdot \Delta t)$ for $k = 0, 1, \ldots, N/2 - 1$.
Step 3. Visualize PSD in both linear and log (dB) scales
```matlab
figure('Position', [100, 100, 1200, 400]);
% Linear scale
subplot(1, 2, 1);
plot(freqs, psd, 'b-', 'LineWidth', 0.8);
xlabel('Frequency (cycles/year)');
ylabel('PSD (mm²·year)');
title('Power Spectral Density - Linear Scale');
xlim([0, 5]);
grid on;
% dB scale: 10*log10(psd/max) so the dominant peak = 0 dB
subplot(1, 2, 2);
psd_dB = 10 * log10(psd / max(psd));
plot(freqs, psd_dB, 'b-', 'LineWidth', 0.8);
xlabel('Frequency (cycles/year)');
ylabel('PSD (dB relative to peak)');
title('Power Spectral Density - Logarithmic Scale');
xlim([0, 5]);
ylim([-60, 5]);
grid on;
```
Purpose: Linear scale shows dominant frequencies clearly. dB scale reveals weaker signals hidden by strong ones. The dB value is $10 \log_{10}(\text{PSD}/\text{PSD}_{\max})$, so the dominant peak sits at 0 dB and everything else is negative. This compresses a factor-of-10,000 variation into a readable 40 dB span.
Step 4. Identify prominent frequencies using threshold
```matlab
% Use threshold: bins > 10x median power (excluding DC)
psd_no_dc = psd(2:end);
threshold = 10 * median(psd_no_dc);
% Find ALL bins above threshold (not just peaks — see Step 9 for the distinction)
prominent_idx = find(psd > threshold);
prominent_freqs = freqs(prominent_idx);
prominent_psd = psd(prominent_idx);
% Sort by power
[prominent_psd, sort_idx] = sort(prominent_psd, 'descend');
prominent_freqs = prominent_freqs(sort_idx);
fprintf('Threshold value: %.1f mm^2·year\n', threshold);
fprintf('Number of bins above threshold: %d\n', length(prominent_idx));
fprintf('Prominent frequencies (cycles/year):\n');
for i = 1:length(prominent_freqs)
f = prominent_freqs(i);
p = prominent_psd(i);
if f > 0
period = 1/f;
fprintf(' f = %.3f (period = %.2f years), Power = %.1f\n', ...
f, period, p);
end
end
```
Purpose: Use a data-driven threshold to flag frequency bins with unusually high power.
> **Note:** This threshold is a **value cutoff**, not a rank-based filter. `find(psd > threshold)` returns every bin above $10 \times \text{median}$, regardless of how many that turns out to be — for our data roughly 68 bins. This is different from selecting the "top 68 by power." In a rank-based approach you choose a count first; here you choose a cutoff value first and the count is determined by the data. Crucially, entire **high-power spectral regions** are flagged, including slopes and saddles between actual peaks. A bin sitting in a "valley" between two large peaks can still be far above the global median. We will diagnose this problem in Step 9.
Step 5. Visualize prominent frequencies
```matlab
figure('Position', [100, 100, 800, 400]);
semilogy(freqs, psd, 'b-', 'LineWidth', 1);
hold on;
prominent_freqs_pos = prominent_freqs(prominent_freqs > 0);
prominent_psd_pos = prominent_psd(prominent_freqs > 0);
semilogy(prominent_freqs_pos, prominent_psd_pos, 'ro', 'MarkerSize', 8, ...
'DisplayName', 'Above threshold');
yline(threshold, 'r:', 'LineWidth', 1.2);
xline(1.0, '--', 'Color', [1 0.5 0], 'Alpha', 0.5, 'Label', 'Annual');
xline(2.0, 'g--', 'Alpha', 0.5, 'Label', 'Semi-annual');
xlabel('Frequency (cycles/year)');
ylabel('PSD (mm^2·year, log scale)');
title('Prominent Frequencies in GMSL');
subtitle(sprintf('Red dotted line = 10×median threshold = %.1f', threshold));
xlim([0 3]);
grid on;
legend('PSD', 'Above threshold', 'Location', 'northeast');
```
Purpose: Visual confirmation of the threshold selection.
> **Note:** You will likely see red dots appearing at what look like valleys or slopes in the spectrum. This is **expected and correct**. The FFT spectrum is a discrete set of points connected by lines — there is no actual data between them. A red dot in an apparent "valley" is an FFT bin that (a) sits within a broad high-power region and clears the global threshold, but (b) is lower than its two immediate neighbors so the connecting lines make it look like a trough. This is the fundamental weakness of the threshold approach, and exactly what Steps 9–10 are designed to address.
Step 6. Fourier regression on discovered frequencies
```matlab
% Select non-DC prominent frequencies for regression
regression_freqs = prominent_freqs(prominent_freqs > 0.1 & prominent_freqs < 10);
fprintf('Using %d regression frequencies\n', length(regression_freqs));
fprintf('Regression frequencies: ');
fprintf('%.3f ', regression_freqs);
fprintf('cycles/year\n');
% Build design matrix: intercept + cos/sin pairs for each frequency
X = ones(N, 1 + 2*length(regression_freqs));
for i = 1:length(regression_freqs)
freq = regression_freqs(i);
X(:, 1 + 2*(i-1) + 1) = cos(2*pi*freq*t);
X(:, 1 + 2*(i-1) + 2) = sin(2*pi*freq*t);
end
% Solve normal equations: beta = (X'X)^{-1} X'S
beta = (X' * X) \ (X' * S);
S_seasonal = X * beta;
% Extract oscillatory part (remove mean)
S_oscillations = S_seasonal - beta(1);
fprintf('Mean sea level: %.1f mm\n', beta(1));
```
Purpose: Use FFT-discovered frequencies for regression, separating mean from oscillations. Each frequency $f_k$ contributes two columns — $\cos(2\pi f_k t)$ and $\sin(2\pi f_k t)$ — representing an oscillation of arbitrary amplitude and phase at that frequency.
Step 7. Remove seasonal and fit quadratic trend
```matlab
% Remove seasonal oscillations
S_deseasonalized = S - S_oscillations;
% Center time for numerical stability
t_mean = mean(t);
t_centered = t - t_mean;
% Fit quadratic trend: S = a + b*t_c + c*t_c^2
A = [ones(N, 1), t_centered, t_centered.^2];
coeffs = (A' * A) \ (A' * S_deseasonalized);
S_trend = A * coeffs;
a = coeffs(1);
b = coeffs(2);
c = coeffs(3);
fprintf('Centered quadratic fit: S = %.1f + %.2f*t_c + %.4f*t_c^2\n', a, b, c);
fprintf('Acceleration: %.4f mm/year^2\n', 2*c);
% Instantaneous rate = b + 2*c*(t - t_mean)
rate_1993 = b + 2*c*(t(1) - t_mean);
rate_2024 = b + 2*c*(t(end) - t_mean);
fprintf('Rate at 1993: %.2f mm/year\n', rate_1993);
fprintf('Rate at 2024: %.2f mm/year\n', rate_2024);
```
Purpose: Quadratic trend captures the acceleration of sea level rise. The instantaneous rate of rise is $dS/dt = b + 2c \cdot t_c$. Because $c > 0$, the rate increases over time — sea level is not just rising, it is accelerating. The coefficient $c$ is half the acceleration: acceleration $= 2c$ mm/year².
Step 8. Compare with Intuitive Fourier Regression (Day 13)
```matlab
% Day 13 - Intuitive Fourier Regression: Pre-specified k=1,2 with quadratic trend
X_intuitive = [ones(N,1), t_centered, t_centered.^2, ...
cos(2*pi*1*t), sin(2*pi*1*t), ...
cos(2*pi*2*t), sin(2*pi*2*t)];
beta_intuitive = (X_intuitive' * X_intuitive) \ (X_intuitive' * S);
S_intuitive = X_intuitive * beta_intuitive;
% Day 16 - FFT-discovered: quadratic trend + oscillations
S_prominent = S_trend + S_oscillations;
% Compute R²
SS_tot = sum((S - mean(S)).^2);
R2_intuitive = 1 - sum((S - S_intuitive).^2) / SS_tot;
R2_prominent = 1 - sum((S - S_prominent).^2) / SS_tot;
fprintf('\nModel Comparison:\n');
fprintf('Day 13 - OG Fourier Regression (2 freqs): R² = %.4f\n', R2_intuitive);
fprintf('Day 16 - Prominent Frequencies (%d freqs): R² = %.4f\n', ...
length(regression_freqs), R2_prominent);
fprintf('Improvement: dR² = %.4f\n', R2_prominent - R2_intuitive);
% Visualize first 5 years
idx = t < t(1) + 5;
figure('Position', [100, 100, 800, 400]);
plot(t(idx), S(idx), 'ko', 'MarkerSize', 3, 'MarkerFaceColor', 'k');
hold on;
plot(t(idx), S_intuitive(idx), 'b-', 'LineWidth', 1.5, ...
'DisplayName', sprintf('OG Fourier Regression (2 freqs, R²=%.3f)', R2_intuitive));
plot(t(idx), S_prominent(idx), 'r-', 'LineWidth', 1.5, ...
'DisplayName', sprintf('Prominent Frequencies (%d freqs, R²=%.3f)', ...
length(regression_freqs), R2_prominent));
xlabel('Year'); ylabel('Sea Level (mm)');
title('Model Comparison: OG Fourier vs FFT-discovered Frequencies');
legend('Location', 'northwest');
grid on;
```
Purpose: Compare FFT-discovered frequencies with our physics-based Day 13 approach.
Step 9. The problem with threshold approach
```matlab
fprintf('\nThreshold approach found %d bins above threshold — not all are peaks!\n', ...
length(prominent_freqs(prominent_freqs > 0.1)));
fprintf('But are they all real peaks? Let''s investigate...\n');
figure('Position', [100, 100, 800, 400]);
semilogy(freqs, psd, 'b-', 'LineWidth', 0.8, 'Alpha', 0.5);
hold on;
semilogy(prominent_freqs(prominent_freqs > 0), ...
prominent_psd(prominent_freqs > 0), 'ro', 'MarkerSize', 4, 'Alpha', 0.5);
yline(threshold, 'r:', 'LineWidth', 1.2);
xlabel('Frequency (cycles/year)');
ylabel('PSD (log scale)');
title('Problem: Threshold selects high-power regions, not peaks!');
subtitle(sprintf('Threshold = %.1f; %d bins above threshold', ...
threshold, length(prominent_idx)));
xlim([0 3]);
grid on;
fprintf('Notice: Many significant points are on slopes, not at peaks!\n');
```
Purpose: Visualize why the threshold approach picks up too many frequencies.
> **Note:** The subtitle says something like "68 bins above threshold." Of those, only a fraction are actual local maxima. The rest are on slopes and broad shoulders around the annual and semi-annual peaks. A value cutoff cannot distinguish the summit from the hillside — it just flags everything above a certain altitude. This is the central limitation we address in Step 10.
Step 10. Smart peak detection approach
```matlab
% Focus on oscillations — zero out the low-frequency trend region
mask = freqs > 0.3;
psd_for_peaks = psd;
psd_for_peaks(~mask) = 0;
% Find local maxima with minimum separation and prominence requirement
min_distance = round(0.2 / freqs(2)); % 0.2 cycles/year minimum separation
peaks = [];
for i = 2:(length(psd_for_peaks)-1)
% Must be a local maximum
if psd_for_peaks(i) > psd_for_peaks(i-1) && psd_for_peaks(i) > psd_for_peaks(i+1)
% Must have minimum separation from last detected peak
if isempty(peaks) || (i - peaks(end)) >= min_distance
% Prominence check: must stand out above local background
window = 10;
local_start = max(1, i - window);
local_end = min(length(psd_for_peaks), i + window);
local_bg = median(psd_for_peaks(local_start:local_end));
if (psd_for_peaks(i) - local_bg) > 20
peaks = [peaks, i];
end
end
end
end
% [Signal-to-noise ratio](https://en.wikipedia.org/wiki/Signal-to-noise_ratio) for all bins
window_size = round(0.5 / freqs(2));
background = zeros(size(psd));
for i = 1:length(psd)
start_idx = max(1, i - floor(window_size/2));
end_idx = min(length(psd), i + floor(window_size/2));
background(i) = prctile(psd(start_idx:end_idx), 25);
end
peak_ratio = psd ./ (background + 1e-10);
% Extract peak properties
peak_freqs = freqs(peaks);
peak_powers = psd(peaks);
peak_ratios = peak_ratio(peaks);
fprintf('\nFound %d distinct peaks:\n', length(peaks));
fprintf('Freq (cy/yr) | Period (yr) | Power\n');
fprintf('-------------|-------------|-------\n');
for i = 1:length(peak_freqs)
fprintf('%12.3f | %11.2f | %6.0f\n', peak_freqs(i), 1/peak_freqs(i), peak_powers(i));
end
% Filter: keep oscillatory peaks (f > 0.5), high signal/background, physically reasonable
retained_mask = (peak_freqs > 0.5) & (peak_freqs < 4.0) & (peak_ratios > 5);
trend_mask = peak_freqs < 0.5;
retained_freqs = peak_freqs(retained_mask);
retained_psd = peak_powers(retained_mask);
rejected_freqs = peak_freqs(~retained_mask & ~trend_mask);
rejected_psd = peak_powers(~retained_mask & ~trend_mask);
fprintf('\nRetained %d peaks out of %d detected\n', sum(retained_mask), length(peak_freqs));
if sum(retained_mask) > 0
fprintf('Retained peaks (cycles/year): ');
fprintf('%.3f ', retained_freqs);
fprintf('\n');
end
% Visualization
threshold_mask = (psd > threshold) & (freqs > 0.5)';
figure('Position', [100, 100, 900, 400]);
semilogy(freqs, psd, 'b-', 'LineWidth', 0.8, 'Alpha', 0.7, 'DisplayName', 'PSD');
hold on;
if any(threshold_mask)
h = semilogy(freqs(threshold_mask), psd(threshold_mask), 'ro', 'MarkerSize', 4);
h.Color(4) = 0.2;
h.DisplayName = sprintf('Above threshold (%d points)', sum(threshold_mask));
end
% Grey square for trend peak
low_freq_idx = freqs > 0 & freqs < 0.3;
if any(low_freq_idx)
[~, tmp] = max(psd(low_freq_idx));
trend_idx_local = find(low_freq_idx);
semilogy(freqs(trend_idx_local(tmp)), psd(trend_idx_local(tmp)), 'ks', ...
'MarkerSize', 10, 'MarkerFaceColor', [0.6 0.6 0.6], ...
'DisplayName', 'Trend (not oscillation)');
end
if ~isempty(rejected_freqs)
h2 = semilogy(rejected_freqs, rejected_psd, 'go', 'MarkerSize', 8);
h2.Color(4) = 0.3;
h2.DisplayName = sprintf('Detected, not used (%d)', length(rejected_freqs));
end
if ~isempty(retained_freqs)
semilogy(retained_freqs, retained_psd, 'go', 'MarkerSize', 8, ...
'MarkerFaceColor', 'g', 'DisplayName', sprintf('Used in model (%d)', ...
length(retained_freqs)));
end
xline(1.0, '--', 'Color', [1 0.5 0], 'Alpha', 0.5, 'Label', 'Annual');
xline(2.0, 'g--', 'Alpha', 0.5, 'Label', 'Semi-annual');
xlim([0 3]);
xlabel('Frequency (cycles/year)');
ylabel('PSD (log scale)');
title('Peak Detection Process: From Threshold to Smart Selection');
legend('Location', 'northeast', 'FontSize', 8);
grid on;
```
Purpose: Proper peak detection requires a bin to be a local maximum **and** stand above its local background. This combination eliminates slopes and shoulders that the global threshold included.
Step 11. Model comparison with detected peaks
```matlab
% Helper: fit Fourier + quadratic model and return metrics
function result = fit_model(freq_list, t, S, N, t_centered, label)
X = ones(N, 1 + 2*length(freq_list));
for i = 1:length(freq_list)
X(:, 1 + 2*(i-1) + 1) = cos(2*pi*freq_list(i)*t);
X(:, 1 + 2*(i-1) + 2) = sin(2*pi*freq_list(i)*t);
end
beta = (X' * X) \ (X' * S);
S_osc = X * beta - beta(1);
A = [ones(N,1), t_centered, t_centered.^2];
coeffs = (A' * A) \ (A' * (S - S_osc));
S_full = A * coeffs + S_osc;
SS_tot = sum((S - mean(S)).^2);
SS_res = sum((S - S_full).^2);
R2 = 1 - SS_res/SS_tot;
fprintf('%-30s | R² = %.4f\n', label, R2);
result = struct('R2', R2, 'S_fit', S_full);
end
fprintf('\nModel Comparison:\n');
fprintf('%s\n', repmat('-', 1, 45));
n_thresh = min(length(regression_freqs), 20);
result_thresh = fit_model(regression_freqs(1:n_thresh), t, S, N, t_centered, ...
sprintf('Threshold (%d freqs)', length(regression_freqs)));
if ~isempty(peak_freqs)
result_peaks = fit_model(peak_freqs, t, S, N, t_centered, ...
sprintf('Localized peaks (%d freqs)', length(peak_freqs)));
end
% Snap to nearest FFT bin — the grid may not include exactly 1.000 or 2.000
[~, idx1] = min(abs(freqs - 1.0));
[~, idx2] = min(abs(freqs - 2.0));
og_freqs = [freqs(idx1), freqs(idx2)];
result_og = fit_model(og_freqs, t, S, N, t_centered, 'OG Fourier (annual + semi-annual)');
fprintf('\nConclusion: OG Fourier Regression (2 frequencies) remains optimal.\n');
fprintf('FFT validates our physical intuition from Day 13!\n');
```
Purpose: Compare models using [$R^2$](https://en.wikipedia.org/wiki/Coefficient_of_determination) to evaluate how well each frequency selection fits the data.
> **Note:** `[~, idx1] = min(abs(freqs - 1.0))` snaps to the FFT bin nearest to 1.0 cycles/year rather than assuming the grid includes exactly 1.000. The bins are spaced $\Delta f \approx 0.024$ cycles/year apart and may land at 0.997 or 1.021. Using the nearest actual bin keeps the regression frequencies consistent with the FFT output.
---
# Results and Interpretation
## Discovered Frequencies
The FFT reveals dominant peaks at:
1. **~1.0 cycles/year**: Annual cycle from Earth's orbit
2. **~2.0 cycles/year**: Semi-annual from hemispheric asymmetry
3. **Near 0**: Long-term trend (DC component)
4. **No clear ENSO peak**: Because it's quasi-periodic (2-7 years), not a single frequency
## Comparison with Day 13
| Approach | Frequencies Used | R² | Physical Basis |
|----------|-----------------|-----|---------------|
| Day 13 | Pre-specified (1, 2 cy/yr) | ~0.92-0.94 | Domain knowledge |
| FFT-based | Discovered (~1.0, ~2.0 cy/yr) | ~0.92-0.94 | Data-driven |
**Key Finding**: The FFT confirms our physical intuition! The dominant frequencies match what we expected from Earth's orbital mechanics.
## [Power Spectrum](https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density) vs [Power Spectral Density](https://en.wikipedia.org/wiki/Spectral_density)
- **Power Spectrum** ($|FFT|^2$): Total power at each frequency bin
- **[PSD](https://en.wikipedia.org/wiki/Spectral_density)** (Power/Hz): Power per unit frequency [bandwidth](https://en.wikipedia.org/wiki/Bandwidth_(signal_processing))
- PSD is preferred for comparing spectra with different sampling rates
- Both show the same peaks, but PSD has consistent units
## Why No Clear ENSO Signal?
[El Niño-Southern Oscillation](https://www.climate.gov/enso) doesn't appear as a single peak because:
1. **Irregular period**: Varies between 2-7 years
2. **[Non-stationary](https://en.wikipedia.org/wiki/Stationary_process)**: Amplitude and frequency change over time
3. **[Broad-band signal](https://en.wikipedia.org/wiki/Broadband)**: Energy spread across multiple frequencies
4. Next class: Spectrograms will reveal time-varying frequencies!
# Check Your Understanding
The following was all AI-generated, with no prompting, after interacting with my R-codebase to produce these multi-language code blocks. It's not really what I was getting at with the heavy graphing of this code base, but I thought it was neat, and so I kept it.
## 1. Nyquist Frequency
What's the highest frequency we can detect with monthly sampling?
```python
# Sampling rate = 12 samples/year
fs = 12 # samples/year
f_nyquist = fs / 2
print(f"Nyquist frequency: {f_nyquist} cycles/year")
print(f"Period: {1/f_nyquist:.1f} months")
# Frequencies above 6 cycles/year are aliased!
```
## 2. Frequency Resolution
How precisely can we determine frequencies?
```python
# Resolution depends on total observation time
T_total = t[-1] - t[0] # years
df = 1 / T_total
print(f"Frequency resolution: {df:.4f} cycles/year")
print(f"Can distinguish periods differing by ~{365/T_total:.1f} days")
```
## 3. [Parseval's Theorem](https://en.wikipedia.org/wiki/Parseval%27s_theorem)
Verify that energy is conserved between time and frequency domains:
```python
# Time domain energy
E_time = np.sum(S_detrended**2)
# Frequency domain energy (need full FFT)
fft_full = np.fft.fft(S_detrended)
E_freq = np.sum(np.abs(fft_full)**2) / N
print(f"Time domain energy: {E_time:.1f}")
print(f"Frequency domain energy: {E_freq:.1f}")
print(f"Ratio: {E_freq/E_time:.6f}") # Should be ~1.0
```
## 4. Phase Information
The FFT gives both magnitude and phase. Extract the phase of the annual component:
```python
# Find annual frequency bin
annual_idx = np.argmin(np.abs(freqs - 1.0))
annual_fft = fft_vals[annual_idx]
magnitude = np.abs(annual_fft)
phase = np.angle(annual_fft)
print(f"Annual component: magnitude = {magnitude:.1f}")
print(f"Phase = {phase:.2f} radians = {phase*180/np.pi:.1f} degrees")
# Phase tells us when in the year the peak occurs
```
## 5. [Zero Padding](https://en.wikipedia.org/wiki/Zero_padding)
What happens if we add zeros to make the FFT longer?
```python
# Pad to next power of 2 for faster FFT
N_padded = 2**int(np.ceil(np.log2(N)))
S_padded = np.pad(S_detrended, (0, N_padded - N), mode='constant')
fft_padded = np.fft.fft(S_padded)
freqs_padded = np.fft.fftfreq(N_padded, dt)[:N_padded//2]
print(f"Original: {N} points, {len(freqs)} frequencies")
print(f"Padded: {N_padded} points, {len(freqs_padded)} frequencies")
print("Zero padding interpolates spectrum (smoother plot) but doesn't add information!")
```
# Key Takeaways
1. **FFT as Discovery Tool**: Let the data tell you which frequencies matter
2. **Validation of Physics**: FFT confirms annual/semi-annual cycles we expected
3. **[Spectrum is Discrete](https://en.wikipedia.org/wiki/Discrete_spectrum)**: FFT returns exactly $N/2$ numbers — plotting libraries connect the dots
4. **Threshold ≠ Top-N**: Value cutoff flags entire regions including slopes and shoulders
5. **[dB Scale](https://en.wikipedia.org/wiki/Decibel) Reveals Hidden Signals**: [Logarithmic scale](https://en.wikipedia.org/wiki/Logarithmic_scale) compresses [dynamic range](https://en.wikipedia.org/wiki/Dynamic_range)
6. **[Acceleration](https://en.wikipedia.org/wiki/Acceleration) Matters**: Sea level rise is accelerating at ~0.08 mm/year²
7. **Smart [Peak Detection](https://en.wikipedia.org/wiki/Peak_detection) Prevents [Overfitting](https://en.wikipedia.org/wiki/Overfitting)**: Local maxima + [prominence](https://en.wikipedia.org/wiki/Topographic_prominence) beats global threshold
8. **[Model Selection](https://en.wikipedia.org/wiki/Model_selection): Simpler is Better**: OG Fourier Regression (2 frequencies) remains optimal
# Mathematical Connections
- **From Day 13**: Pre-specified frequencies were physically motivated
- **From Day 14-15**: Delta functions connect continuous and discrete transforms
- **To [Spectrograms](https://en.wikipedia.org/wiki/Spectrogram)**: [Time-frequency analysis](https://en.wikipedia.org/wiki/Time%E2%80%93frequency_analysis) for non-stationary signals
- **To Filtering**: Zeroing frequencies = ideal filter (preview of convolution)
- **To Signal Processing**: [Nyquist theorem](https://en.wikipedia.org/wiki/Nyquist_frequency), [aliasing](https://en.wikipedia.org/wiki/Aliasing)
# Preview: Connection to Filtering
The FFT enables frequency-domain filtering:
1. **Transform** to frequency domain (FFT)
2. **Modify** frequencies (multiply by filter response)
3. **Transform back** to time domain (inverse FFT)
This is mathematically equivalent to [convolution](https://en.wikipedia.org/wiki/Convolution) in the time domain! Next class: How this relates to [RC circuits](https://en.wikipedia.org/wiki/RC_circuit) and physical [filters](https://en.wikipedia.org/wiki/Filter_(signal_processing)).
---
*Remember: The FFT is just a fast algorithm for the [DFT](https://en.wikipedia.org/wiki/Discrete_Fourier_transform). The mathematics are the same, but [O(N log N)](https://en.wikipedia.org/wiki/Time_complexity) beats O(N²) every time!*