Spectral Analysis
Spectral Analysis¶
Learning Objectives¶
- Understand the definition and properties of the power spectral density (PSD)
- Master the Wiener-Khinchin theorem relating autocorrelation and PSD
- Learn non-parametric spectral estimation methods (periodogram, Bartlett, Welch, Blackman-Tukey)
- Understand the fundamental resolution-variance tradeoff in spectral estimation
- Learn parametric spectral estimation methods (AR, MA, ARMA models)
- Apply model order selection criteria (AIC, BIC)
- Perform cross-spectral analysis and compute coherence
- Create and interpret spectrograms for time-varying spectral analysis
- Implement spectral analysis techniques using Python's
scipy.signalandmatplotlib
Table of Contents¶
- Introduction to Spectral Analysis
- Power Spectral Density
- Wiener-Khinchin Theorem
- Periodogram
- Bartlett's Method
- Welch's Method
- Blackman-Tukey Method
- Resolution-Variance Tradeoff
- Parametric Methods: AR Models
- Parametric Methods: MA and ARMA
- Model Order Selection
- Cross-Spectral Analysis and Coherence
- Spectrogram
- Python Implementation
- Exercises
1. Introduction to Spectral Analysis¶
1.1 The Problem¶
Given a finite observation of a random or deterministic signal, we want to estimate how the signal's power (or energy) is distributed across frequencies. This is the spectral estimation problem.
Unlike the DFT, which computes the exact spectral representation of a finite deterministic signal, spectral analysis deals with: - Random processes: signals with statistical properties rather than exact values - Finite data: we only observe $N$ samples of a potentially infinite-duration process - Estimation uncertainty: any estimate has bias and variance
1.2 Applications¶
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
โ Spectral Analysis Applications โ
โโโโโโโโโโโโโโโโโโโโโโฌโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโค
โ Audio/Speech โ Pitch detection, formant analysis, โ
โ โ noise characterization โ
โโโโโโโโโโโโโโโโโโโโโโผโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโค
โ Communications โ Channel characterization, interference โ
โ โ detection, modulation classification โ
โโโโโโโโโโโโโโโโโโโโโโผโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโค
โ Biomedical โ EEG/ECG analysis, sleep staging, โ
โ โ heart rate variability โ
โโโโโโโโโโโโโโโโโโโโโโผโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโค
โ Radar/Sonar โ Target detection, Doppler estimation โ
โโโโโโโโโโโโโโโโโโโโโโผโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโค
โ Vibration Analysis โ Machine health monitoring, modal analysis โ
โโโโโโโโโโโโโโโโโโโโโโผโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโค
โ Geophysics โ Seismic analysis, ocean wave spectra โ
โโโโโโโโโโโโโโโโโโโโโโผโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโค
โ Astrophysics โ Pulsar timing, gravitational wave detection โ
โโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
1.3 Two Approaches¶
Non-parametric methods: Make minimal assumptions about the signal. Estimate the PSD directly from the data (e.g., periodogram, Welch's method).
Parametric methods: Assume the signal is generated by a model (e.g., AR, MA, ARMA). Estimate model parameters, then compute the PSD from the model.
2. Power Spectral Density¶
2.1 Definition for Random Processes¶
For a wide-sense stationary (WSS) random process $x[n]$, the power spectral density (PSD) is defined as:
$$S_x(e^{j\omega}) = \sum_{k=-\infty}^{\infty} r_{xx}[k] e^{-j\omega k}$$
where $r_{xx}[k] = E\{x[n] x^*[n-k]\}$ is the autocorrelation sequence.
2.2 Properties of the PSD¶
- Real-valued: $S_x(e^{j\omega}) \in \mathbb{R}$ for all $\omega$
- Non-negative: $S_x(e^{j\omega}) \geq 0$ for all $\omega$
- Periodic: $S_x(e^{j\omega}) = S_x(e^{j(\omega + 2\pi)})$
- Even symmetry (for real $x[n]$): $S_x(e^{j\omega}) = S_x(e^{-j\omega})$
- Total power: $r_{xx}[0] = E\{|x[n]|^2\} = \frac{1}{2\pi} \int_{-\pi}^{\pi} S_x(e^{j\omega}) d\omega$
2.3 PSD of Filtered Processes¶
If $y[n] = h[n] * x[n]$ (LTI filtering), then:
$$S_y(e^{j\omega}) = |H(e^{j\omega})|^2 S_x(e^{j\omega})$$
This is the fundamental relationship between input and output power spectra.
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fft import fft, fftfreq
def psd_of_filtered_noise():
"""Demonstrate PSD of filtered white noise."""
np.random.seed(42)
N = 10000
fs = 1000
# White noise: flat PSD
x = np.random.randn(N)
# Design bandpass filter
sos = signal.butter(4, [100, 200], btype='bandpass', fs=fs, output='sos')
y = signal.sosfilt(sos, x)
# Estimate PSDs using Welch's method
f_x, Pxx = signal.welch(x, fs=fs, nperseg=512)
f_y, Pyy = signal.welch(y, fs=fs, nperseg=512)
# Theoretical: |H(f)|^2 * Pxx
w, H = signal.sosfreqz(sos, worN=len(f_x), fs=fs)
Pyy_theory = np.abs(H)**2 * np.mean(Pxx) # Pxx โ constant for white noise
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Input PSD
axes[0, 0].semilogy(f_x, Pxx, 'b-', linewidth=1.5)
axes[0, 0].set_title('Input PSD (White Noise)')
axes[0, 0].set_xlabel('Frequency (Hz)')
axes[0, 0].set_ylabel('PSD (Vยฒ/Hz)')
axes[0, 0].grid(True, alpha=0.3)
# Filter response
axes[0, 1].plot(w, 20*np.log10(np.abs(H) + 1e-15), 'g-', linewidth=1.5)
axes[0, 1].set_title('|H(f)|ยฒ (Bandpass Filter)')
axes[0, 1].set_xlabel('Frequency (Hz)')
axes[0, 1].set_ylabel('Magnitude (dB)')
axes[0, 1].grid(True, alpha=0.3)
# Output PSD: measured vs theoretical
axes[1, 0].semilogy(f_y, Pyy, 'r-', linewidth=1.5, label='Estimated')
axes[1, 0].semilogy(w, Pyy_theory, 'k--', linewidth=1.5, label='Theoretical')
axes[1, 0].set_title('Output PSD: Sy = |H|ยฒ ยท Sx')
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('PSD (Vยฒ/Hz)')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# Time domain
t = np.arange(N) / fs
axes[1, 1].plot(t[:500], x[:500], 'b-', alpha=0.5, label='Input')
axes[1, 1].plot(t[:500], y[:500], 'r-', linewidth=1.5, label='Output')
axes[1, 1].set_title('Time Domain')
axes[1, 1].set_xlabel('Time (s)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.suptitle('PSD of Filtered Process: Sy(f) = |H(f)|ยฒ Sx(f)', fontsize=14)
plt.tight_layout()
plt.savefig('psd_filtered_noise.png', dpi=150)
plt.close()
psd_of_filtered_noise()
3. Wiener-Khinchin Theorem¶
3.1 Statement¶
The Wiener-Khinchin theorem states that for a wide-sense stationary process, the power spectral density and the autocorrelation function form a Fourier transform pair:
$$S_x(e^{j\omega}) = \sum_{k=-\infty}^{\infty} r_{xx}[k] e^{-j\omega k} \quad \text{(DTFT)}$$
$$r_{xx}[k] = \frac{1}{2\pi} \int_{-\pi}^{\pi} S_x(e^{j\omega}) e^{j\omega k} d\omega \quad \text{(inverse DTFT)}$$
3.2 Significance¶
This theorem provides two fundamental approaches to spectral estimation:
- Direct approach: Estimate the PSD directly from the data (periodogram methods)
- Indirect approach: First estimate the autocorrelation $\hat{r}_{xx}[k]$, then take its Fourier transform (Blackman-Tukey method)
3.3 Verification¶
def verify_wiener_khinchin():
"""Verify the Wiener-Khinchin theorem numerically."""
np.random.seed(42)
N = 50000
fs = 1000
# Generate colored noise (AR process)
# x[n] = 0.9*x[n-1] + w[n]
a = 0.9
w = np.random.randn(N)
x = np.zeros(N)
x[0] = w[0]
for n in range(1, N):
x[n] = a * x[n-1] + w[n]
# Method 1: Direct PSD estimation (periodogram)
f_direct, Pxx_direct = signal.periodogram(x, fs=fs)
# Method 2: Via autocorrelation (Wiener-Khinchin)
max_lag = 500
r_xx = np.correlate(x, x, mode='full')[N-1-max_lag:N+max_lag] / N
lags = np.arange(-max_lag, max_lag + 1)
# DTFT of autocorrelation
omega = np.linspace(0, np.pi, len(f_direct))
Pxx_wk = np.zeros(len(omega))
for i, w_val in enumerate(omega):
Pxx_wk[i] = np.real(np.sum(r_xx * np.exp(-1j * w_val * lags))) / fs
# Theoretical PSD for AR(1): sigma_w^2 / |1 - a*e^(-jw)|^2
sigma_w = 1.0
Pxx_theory = sigma_w**2 / (np.abs(1 - a * np.exp(-1j * omega))**2) / fs
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Autocorrelation
axes[0, 0].plot(lags, r_xx, 'b-', linewidth=1.5)
axes[0, 0].set_title('Estimated Autocorrelation r_xx[k]')
axes[0, 0].set_xlabel('Lag k')
axes[0, 0].set_ylabel('r_xx[k]')
axes[0, 0].grid(True, alpha=0.3)
# Theoretical autocorrelation: r_xx[k] = sigma_w^2 * a^|k| / (1-a^2)
r_theory = sigma_w**2 * a**np.abs(lags) / (1 - a**2)
axes[0, 1].plot(lags, r_xx, 'b-', alpha=0.5, label='Estimated')
axes[0, 1].plot(lags, r_theory, 'r--', linewidth=1.5, label='Theoretical')
axes[0, 1].set_title('Autocorrelation: Estimated vs Theoretical')
axes[0, 1].set_xlabel('Lag k')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# PSD comparison
f_hz = omega / (2 * np.pi) * fs
axes[1, 0].semilogy(f_direct, Pxx_direct, 'b-', alpha=0.3, label='Periodogram')
axes[1, 0].semilogy(f_hz, Pxx_wk, 'g-', linewidth=2, label='Wiener-Khinchin')
axes[1, 0].semilogy(f_hz, Pxx_theory, 'r--', linewidth=2, label='Theoretical')
axes[1, 0].set_title('PSD Comparison')
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('PSD (Vยฒ/Hz)')
axes[1, 0].legend()
axes[1, 0].grid(True, which='both', alpha=0.3)
# PSD in dB
axes[1, 1].plot(f_hz, 10*np.log10(Pxx_wk + 1e-15), 'g-', linewidth=2,
label='Wiener-Khinchin')
axes[1, 1].plot(f_hz, 10*np.log10(Pxx_theory + 1e-15), 'r--', linewidth=2,
label='Theoretical')
axes[1, 1].set_title('PSD (dB Scale)')
axes[1, 1].set_xlabel('Frequency (Hz)')
axes[1, 1].set_ylabel('PSD (dB/Hz)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.suptitle('Wiener-Khinchin Theorem Verification', fontsize=14)
plt.tight_layout()
plt.savefig('wiener_khinchin.png', dpi=150)
plt.close()
verify_wiener_khinchin()
4. Periodogram¶
4.1 Definition¶
The periodogram is the simplest PSD estimator, based on the squared magnitude of the DFT:
$$\hat{S}_x^{(\text{per})}(f_k) = \frac{1}{N f_s} |X[k]|^2 = \frac{1}{N f_s} \left| \sum_{n=0}^{N-1} x[n] e^{-j 2\pi k n / N} \right|^2$$
where $f_k = k f_s / N$ and $k = 0, 1, \ldots, N-1$.
4.2 Statistical Properties¶
Bias: The expected value of the periodogram is the true PSD convolved with the Fejer kernel:
$$E\left\{\hat{S}_x^{(\text{per})}(e^{j\omega})\right\} = \frac{1}{2\pi} S_x(e^{j\omega}) * W_N(e^{j\omega})$$
where $W_N(e^{j\omega}) = \frac{1}{N}\left|\frac{\sin(N\omega/2)}{\sin(\omega/2)}\right|^2$ is the Fejer kernel (squared Dirichlet kernel).
- As $N \to \infty$, the bias decreases (asymptotically unbiased)
- For finite $N$, there is spectral leakage due to the finite observation window
Variance: The periodogram is not consistent -- its variance does not decrease with increasing $N$:
$$\text{Var}\left\{\hat{S}_x^{(\text{per})}(e^{j\omega})\right\} \approx S_x^2(e^{j\omega})$$
This means the periodogram fluctuates wildly regardless of data length, making it unreliable for smooth PSD estimation.
4.3 Implementation and Demonstration¶
def periodogram_analysis():
"""Analyze the properties of the periodogram estimator."""
np.random.seed(42)
fs = 1000
# True signal: two sinusoids in noise
f1, f2 = 100, 120 # Hz (close frequencies)
A1, A2 = 1.0, 0.8
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Different data lengths
for idx, N in enumerate([64, 256, 1024, 4096]):
ax = axes[idx // 2, idx % 2]
t = np.arange(N) / fs
x = (A1 * np.sin(2*np.pi*f1*t) +
A2 * np.sin(2*np.pi*f2*t) +
np.random.randn(N) * 0.5)
# Periodogram
f, Pxx = signal.periodogram(x, fs=fs)
ax.semilogy(f, Pxx, 'b-', linewidth=0.8, alpha=0.8)
# Mark true frequencies
ax.axvline(f1, color='r', linestyle='--', alpha=0.5, label=f'{f1} Hz')
ax.axvline(f2, color='g', linestyle='--', alpha=0.5, label=f'{f2} Hz')
# Frequency resolution
delta_f = fs / N
ax.set_title(f'N={N}, ฮf = {delta_f:.1f} Hz '
f'({"resolved" if delta_f < abs(f2-f1) else "NOT resolved"})')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD (Vยฒ/Hz)')
ax.set_xlim(0, 300)
ax.legend(fontsize=8)
ax.grid(True, which='both', alpha=0.3)
plt.suptitle('Periodogram: Effect of Data Length on Resolution', fontsize=14)
plt.tight_layout()
plt.savefig('periodogram_analysis.png', dpi=150)
plt.close()
periodogram_analysis()
4.4 Periodogram Variance Demonstration¶
def periodogram_variance():
"""Show that periodogram variance doesn't decrease with N."""
np.random.seed(42)
fs = 1000
# AR(1) process (known PSD)
a = 0.9
sigma_w = 1.0
num_realizations = 50
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for ax, N in zip(axes, [256, 1024, 4096]):
for trial in range(num_realizations):
w = np.random.randn(N) * sigma_w
x = np.zeros(N)
x[0] = w[0]
for n in range(1, N):
x[n] = a * x[n-1] + w[n]
f, Pxx = signal.periodogram(x, fs=fs)
ax.semilogy(f, Pxx, 'b-', alpha=0.1, linewidth=0.5)
# Theoretical PSD
f_theory = np.linspace(0, fs/2, 500)
omega = 2 * np.pi * f_theory / fs
Pxx_theory = sigma_w**2 / (np.abs(1 - a*np.exp(-1j*omega))**2) / fs
ax.semilogy(f_theory, Pxx_theory, 'r-', linewidth=2, label='True PSD')
ax.set_title(f'N = {N} ({num_realizations} realizations)')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD (Vยฒ/Hz)')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
plt.suptitle('Periodogram Variance Does NOT Decrease with N', fontsize=14)
plt.tight_layout()
plt.savefig('periodogram_variance.png', dpi=150)
plt.close()
periodogram_variance()
5. Bartlett's Method¶
5.1 Principle¶
Bartlett's method (1948) reduces variance by averaging multiple periodograms:
- Divide the data into $K$ non-overlapping segments of length $L$ ($N = KL$)
- Compute the periodogram of each segment
- Average the $K$ periodograms
$$\hat{S}_x^{(\text{Bartlett})}(f) = \frac{1}{K} \sum_{i=0}^{K-1} \hat{S}_{x_i}^{(\text{per})}(f)$$
5.2 Properties¶
Variance reduction: By averaging $K$ independent periodograms:
$$\text{Var}\left\{\hat{S}_x^{(\text{Bartlett})}\right\} \approx \frac{1}{K} S_x^2(e^{j\omega})$$
Resolution loss: Each segment is shorter ($L = N/K$), so the frequency resolution is:
$$\Delta f = \frac{f_s}{L} = \frac{K f_s}{N}$$
This is the resolution-variance tradeoff: reducing variance by factor $K$ degrades resolution by factor $K$.
5.3 Implementation¶
def bartlett_method(x, fs, K):
"""
Bartlett's method for PSD estimation.
Parameters:
x: input signal
fs: sampling frequency
K: number of non-overlapping segments
Returns:
f: frequency vector
Pxx: PSD estimate
"""
N = len(x)
L = N // K # Segment length
# Compute periodogram of each segment
f = np.fft.rfftfreq(L, 1/fs)
Pxx_sum = np.zeros(len(f))
for i in range(K):
segment = x[i*L : (i+1)*L]
X_seg = np.fft.rfft(segment)
Pxx_seg = np.abs(X_seg)**2 / (L * fs)
Pxx_sum += Pxx_seg
Pxx = Pxx_sum / K
return f, Pxx
# Compare periodogram vs Bartlett with different K
np.random.seed(42)
N = 4096
fs = 1000
# AR(1) process
a = 0.9
w = np.random.randn(N)
x = np.zeros(N)
x[0] = w[0]
for n in range(1, N):
x[n] = a * x[n-1] + w[n]
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
configs = [
(1, 'Periodogram (K=1)'),
(4, 'Bartlett (K=4)'),
(16, 'Bartlett (K=16)'),
(64, 'Bartlett (K=64)'),
]
# Theoretical PSD
f_theory = np.linspace(0, fs/2, 500)
omega = 2 * np.pi * f_theory / fs
Pxx_theory = 1.0 / (np.abs(1 - a*np.exp(-1j*omega))**2) / fs
for ax, (K, title) in zip(axes.flat, configs):
f, Pxx = bartlett_method(x, fs, K)
L = N // K
ax.semilogy(f, Pxx, 'b-', linewidth=1, alpha=0.8, label='Estimate')
ax.semilogy(f_theory, Pxx_theory, 'r-', linewidth=2, label='True PSD')
ax.set_title(f'{title}, L={L}, ฮf={fs/L:.1f} Hz')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD (Vยฒ/Hz)')
ax.legend(fontsize=8)
ax.grid(True, which='both', alpha=0.3)
plt.suptitle("Bartlett's Method: Resolution-Variance Tradeoff", fontsize=14)
plt.tight_layout()
plt.savefig('bartlett_method.png', dpi=150)
plt.close()
6. Welch's Method¶
6.1 Improvements over Bartlett¶
Welch's method (1967) extends Bartlett's method with two key improvements:
- Overlapping segments: Segments overlap by a fraction (typically 50%), yielding more segments from the same data
- Windowing: Each segment is multiplied by a window function (e.g., Hanning) to reduce spectral leakage
$$\hat{S}_x^{(\text{Welch})}(f) = \frac{1}{K} \sum_{i=0}^{K-1} \frac{1}{L U} \left| \sum_{n=0}^{L-1} w[n] x_i[n] e^{-j2\pi f n / f_s} \right|^2$$
where $U = \frac{1}{L}\sum_{n=0}^{L-1} w^2[n]$ is the window power normalization.
6.2 Properties¶
- More segments: With 50% overlap and $N$ total samples, we get $K \approx 2N/L - 1$ segments (roughly double Bartlett)
- Lower variance: More averaging plus windowed leakage reduction
- Same resolution: Determined by segment length $L$
6.3 Implementation and Analysis¶
def welch_analysis():
"""Comprehensive Welch's method analysis."""
np.random.seed(42)
N = 8192
fs = 1000
# Signal: two closely spaced tones + broadband noise
t = np.arange(N) / fs
f1, f2 = 100, 108 # 8 Hz apart
x = (np.sin(2*np.pi*f1*t) + 0.7*np.sin(2*np.pi*f2*t) +
0.5*np.random.randn(N))
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
# 1. Effect of segment length (nperseg)
for i, nperseg in enumerate([128, 512, 2048]):
f, Pxx = signal.welch(x, fs=fs, nperseg=nperseg, noverlap=nperseg//2)
axes[0, i].semilogy(f, Pxx, 'b-', linewidth=1.5)
axes[0, i].axvline(f1, color='r', linestyle='--', alpha=0.5)
axes[0, i].axvline(f2, color='g', linestyle='--', alpha=0.5)
axes[0, i].set_title(f'nperseg = {nperseg}, ฮf = {fs/nperseg:.1f} Hz')
axes[0, i].set_xlabel('Frequency (Hz)')
axes[0, i].set_ylabel('PSD')
axes[0, i].set_xlim(50, 200)
axes[0, i].grid(True, which='both', alpha=0.3)
# 2. Effect of overlap
nperseg = 512
for i, overlap_frac in enumerate([0, 0.5, 0.75]):
noverlap = int(nperseg * overlap_frac)
f, Pxx = signal.welch(x, fs=fs, nperseg=nperseg, noverlap=noverlap)
n_segs = (N - noverlap) // (nperseg - noverlap) if noverlap < nperseg else 1
axes[1, i].semilogy(f, Pxx, 'b-', linewidth=1.5)
axes[1, i].axvline(f1, color='r', linestyle='--', alpha=0.5)
axes[1, i].axvline(f2, color='g', linestyle='--', alpha=0.5)
axes[1, i].set_title(f'Overlap = {overlap_frac*100:.0f}%, '
f'K โ {n_segs} segments')
axes[1, i].set_xlabel('Frequency (Hz)')
axes[1, i].set_ylabel('PSD')
axes[1, i].set_xlim(50, 200)
axes[1, i].grid(True, which='both', alpha=0.3)
plt.suptitle("Welch's Method: Segment Length and Overlap Effects", fontsize=14)
plt.tight_layout()
plt.savefig('welch_analysis.png', dpi=150)
plt.close()
welch_analysis()
6.4 Window Effects¶
def welch_window_comparison():
"""Compare different windows in Welch's method."""
np.random.seed(42)
N = 4096
fs = 1000
# Signal with a weak tone near a strong one
t = np.arange(N) / fs
x = (np.sin(2*np.pi*100*t) +
0.01*np.sin(2*np.pi*130*t) + # 40 dB weaker
0.1*np.random.randn(N))
windows = ['boxcar', 'hann', 'hamming', 'blackman', 'blackmanharris']
nperseg = 512
fig, ax = plt.subplots(figsize=(14, 6))
for window in windows:
f, Pxx = signal.welch(x, fs=fs, nperseg=nperseg,
noverlap=nperseg//2, window=window)
ax.plot(f, 10*np.log10(Pxx + 1e-15), linewidth=1.5, label=window)
ax.axvline(100, color='gray', linestyle='--', alpha=0.5)
ax.axvline(130, color='gray', linestyle='--', alpha=0.5)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD (dB/Hz)')
ax.set_title('Window Function Effect: Detecting a Weak Tone (-40 dB)')
ax.set_xlim(50, 200)
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('welch_windows.png', dpi=150)
plt.close()
welch_window_comparison()
7. Blackman-Tukey Method¶
7.1 Principle¶
The Blackman-Tukey method (1958) is an indirect approach based on the Wiener-Khinchin theorem:
- Estimate the autocorrelation sequence $\hat{r}_{xx}[k]$
- Apply a window $w[k]$ to the autocorrelation (to limit the lag range)
- Compute the DTFT of the windowed autocorrelation
$$\hat{S}_x^{(\text{BT})}(e^{j\omega}) = \sum_{k=-(M-1)}^{M-1} w[k] \hat{r}_{xx}[k] e^{-j\omega k}$$
where $M$ is the maximum lag.
7.2 Autocorrelation Estimation¶
The biased autocorrelation estimate:
$$\hat{r}_{xx}[k] = \frac{1}{N} \sum_{n=0}^{N-1-|k|} x[n+|k|] x^*[n], \quad |k| \leq N-1$$
The unbiased estimate:
$$\hat{r}_{xx}^{(\text{unbiased})}[k] = \frac{1}{N - |k|} \sum_{n=0}^{N-1-|k|} x[n+|k|] x^*[n]$$
The biased estimate is preferred because it guarantees a non-negative PSD, while the unbiased estimate can yield negative PSD values.
7.3 Lag Window¶
The lag window $w[k]$ controls the smoothing applied to the PSD estimate: - Narrow window ($M$ small): smooth PSD, low variance, poor resolution - Wide window ($M$ large): rough PSD, high variance, good resolution
Common lag windows: Bartlett (triangular), Parzen, Tukey (cosine-tapered).
7.4 Implementation¶
def blackman_tukey_method(x, fs, max_lag, window='bartlett'):
"""
Blackman-Tukey PSD estimation.
Parameters:
x: input signal
fs: sampling frequency
max_lag: maximum lag for autocorrelation
window: lag window type
Returns:
f: frequency vector
Pxx: PSD estimate
"""
N = len(x)
# Biased autocorrelation estimate
r_xx = np.correlate(x, x, mode='full') / N
center = N - 1
r_xx_lags = r_xx[center - max_lag:center + max_lag + 1]
lags = np.arange(-max_lag, max_lag + 1)
# Apply lag window
if window == 'bartlett':
w = 1 - np.abs(lags) / max_lag
elif window == 'parzen':
u = np.abs(lags) / max_lag
w = np.where(u <= 0.5,
1 - 6*u**2 + 6*u**3,
2*(1 - u)**3)
elif window == 'tukey':
w = 0.5 * (1 + np.cos(np.pi * lags / max_lag))
else:
w = np.ones(len(lags))
r_windowed = r_xx_lags * w
# Compute PSD via DTFT
nfft = 2048
f = np.linspace(0, fs/2, nfft//2 + 1)
omega = 2 * np.pi * f / fs
Pxx = np.zeros(len(f))
for i, w_val in enumerate(omega):
Pxx[i] = np.real(np.sum(r_windowed * np.exp(-1j * w_val * lags))) / fs
# Ensure non-negative
Pxx = np.maximum(Pxx, 0)
return f, Pxx
# Compare with different max_lag values
np.random.seed(42)
N = 4096
fs = 1000
a = 0.9
w = np.random.randn(N)
x = np.zeros(N)
x[0] = w[0]
for n in range(1, N):
x[n] = a * x[n-1] + w[n]
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Theoretical
f_theory = np.linspace(0, fs/2, 500)
omega_t = 2 * np.pi * f_theory / fs
Pxx_theory = 1.0 / (np.abs(1 - a*np.exp(-1j*omega_t))**2) / fs
for ax, max_lag in zip(axes.flat, [16, 64, 256, 1024]):
f_bt, Pxx_bt = blackman_tukey_method(x, fs, max_lag, window='bartlett')
ax.semilogy(f_bt, Pxx_bt, 'b-', linewidth=1.5, label='BT estimate')
ax.semilogy(f_theory, Pxx_theory, 'r--', linewidth=2, label='True PSD')
ax.set_title(f'Blackman-Tukey (max_lag = {max_lag})')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD (Vยฒ/Hz)')
ax.legend(fontsize=8)
ax.grid(True, which='both', alpha=0.3)
plt.suptitle('Blackman-Tukey Method: Effect of Maximum Lag', fontsize=14)
plt.tight_layout()
plt.savefig('blackman_tukey.png', dpi=150)
plt.close()
8. Resolution-Variance Tradeoff¶
8.1 Fundamental Limitation¶
All non-parametric spectral estimators face a fundamental tradeoff:
$$\text{Resolution} \times \text{Variance} \approx \text{constant}$$
- Better frequency resolution requires longer observation windows, which means fewer averages and higher variance
- Lower variance requires more averaging, which means shorter segments and coarser resolution
8.2 Quantitative Analysis¶
For Welch's method with $K$ segments of length $L$:
- Resolution: $\Delta f = f_s / L$
- Normalized variance: $\text{Var}\{\hat{S}\} / S^2 \approx 1/K$
- Product: $\Delta f \times \text{Var}\{\hat{S}\} / S^2 \approx f_s / (KL) = f_s / N$
The resolution bandwidth (RBW) and the equivalent degrees of freedom (EDOF) characterize this tradeoff:
$$\text{EDOF} \approx 2K \quad \text{(for non-overlapping segments)}$$
8.3 Visualization¶
def resolution_variance_tradeoff():
"""Visualize the resolution-variance tradeoff."""
np.random.seed(42)
N = 8192
fs = 1000
# Two tones: 100 Hz and 110 Hz (10 Hz apart)
t = np.arange(N) / fs
x = np.sin(2*np.pi*100*t) + np.sin(2*np.pi*110*t) + 0.5*np.random.randn(N)
# Sweep segment length
nperseg_values = [64, 128, 256, 512, 1024, 2048, 4096]
resolutions = []
variances = []
fig, axes = plt.subplots(2, 4, figsize=(18, 8))
for i, nperseg in enumerate(nperseg_values):
# Multiple realizations to estimate variance
psds = []
for trial in range(50):
x_trial = (np.sin(2*np.pi*100*t) + np.sin(2*np.pi*110*t) +
0.5*np.random.randn(N))
f, Pxx = signal.welch(x_trial, fs=fs, nperseg=nperseg,
noverlap=nperseg//2)
psds.append(Pxx)
psds = np.array(psds)
mean_psd = np.mean(psds, axis=0)
std_psd = np.std(psds, axis=0)
resolutions.append(fs / nperseg)
variances.append(np.mean(std_psd**2 / (mean_psd**2 + 1e-15)))
if i < 7:
ax = axes[i // 4, i % 4]
ax.semilogy(f, mean_psd, 'b-', linewidth=1.5)
ax.fill_between(f, mean_psd - std_psd, mean_psd + std_psd,
alpha=0.2, color='blue')
ax.axvline(100, color='r', linestyle='--', alpha=0.5)
ax.axvline(110, color='g', linestyle='--', alpha=0.5)
ax.set_title(f'L={nperseg}, ฮf={fs/nperseg:.1f} Hz')
ax.set_xlim(50, 200)
ax.set_xlabel('Freq (Hz)')
ax.grid(True, which='both', alpha=0.3)
# Summary plot
ax_summary = axes[1, 3]
ax_summary.loglog(resolutions, variances, 'ro-', markersize=8, linewidth=2)
ax_summary.set_xlabel('Resolution ฮf (Hz)')
ax_summary.set_ylabel('Normalized Variance')
ax_summary.set_title('Resolution-Variance Tradeoff')
ax_summary.grid(True, which='both', alpha=0.3)
plt.suptitle('Resolution vs Variance Tradeoff in Welch\'s Method', fontsize=14)
plt.tight_layout()
plt.savefig('resolution_variance.png', dpi=150)
plt.close()
resolution_variance_tradeoff()
9. Parametric Methods: AR Models¶
9.1 Autoregressive (AR) Model¶
An AR process of order $p$ satisfies:
$$x[n] = -\sum_{k=1}^{p} a_k x[n-k] + w[n]$$
where $w[n]$ is white noise with variance $\sigma_w^2$.
The PSD of an AR($p$) process is:
$$S_x(e^{j\omega}) = \frac{\sigma_w^2}{|A(e^{j\omega})|^2} = \frac{\sigma_w^2}{\left|1 + \sum_{k=1}^{p} a_k e^{-j\omega k}\right|^2}$$
9.2 Why AR Models?¶
AR models are popular because: - The PSD has a smooth, peaked shape good for modeling spectral peaks - Parameter estimation reduces to solving linear equations (Yule-Walker) - They can achieve excellent frequency resolution for short data - Many natural signals (speech, vibration) are well-modeled by AR processes
9.3 Yule-Walker Equations¶
The AR parameters are determined by the Yule-Walker equations:
$$\begin{bmatrix} r_{xx}[0] & r_{xx}[1] & \cdots & r_{xx}[p-1] \\ r_{xx}[1] & r_{xx}[0] & \cdots & r_{xx}[p-2] \\ \vdots & \vdots & \ddots & \vdots \\ r_{xx}[p-1] & r_{xx}[p-2] & \cdots & r_{xx}[0] \end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_p \end{bmatrix} = -\begin{bmatrix} r_{xx}[1] \\ r_{xx}[2] \\ \vdots \\ r_{xx}[p] \end{bmatrix}$$
In matrix form: $\mathbf{R} \mathbf{a} = -\mathbf{r}$
9.4 Levinson-Durbin Algorithm¶
The Levinson-Durbin algorithm solves the Yule-Walker equations efficiently in $O(p^2)$ operations (vs. $O(p^3)$ for direct matrix inversion):
def levinson_durbin(r, order):
"""
Levinson-Durbin recursion for AR parameter estimation.
Parameters:
r: autocorrelation sequence r[0], r[1], ..., r[order]
order: AR model order p
Returns:
a: AR coefficients [a1, a2, ..., ap]
sigma2: driving noise variance
reflection_coeffs: reflection coefficients (PARCOR)
"""
# Initialize
a = np.zeros(order)
reflection_coeffs = np.zeros(order)
# Order 1
a[0] = -r[1] / r[0]
reflection_coeffs[0] = a[0]
sigma2 = r[0] * (1 - a[0]**2)
# Recursion
for m in range(1, order):
# Compute reflection coefficient
numerator = r[m + 1] + np.sum(a[:m] * r[m:0:-1])
k_m = -numerator / sigma2
reflection_coeffs[m] = k_m
# Update AR coefficients
a_new = np.zeros(order)
a_new[m] = k_m
for j in range(m):
a_new[j] = a[j] + k_m * a[m - 1 - j]
a = a_new
# Update variance
sigma2 = sigma2 * (1 - k_m**2)
return a, sigma2, reflection_coeffs
def ar_psd(a, sigma2, fs, nfft=1024):
"""Compute PSD from AR parameters."""
f = np.linspace(0, fs/2, nfft//2 + 1)
omega = 2 * np.pi * f / fs
# A(e^jw) = 1 + a1*e^(-jw) + a2*e^(-j2w) + ...
A_freq = np.ones(len(f), dtype=complex)
for k in range(len(a)):
A_freq += a[k] * np.exp(-1j * (k+1) * omega)
Pxx = sigma2 / (np.abs(A_freq)**2 * fs)
return f, Pxx
9.5 AR Spectral Estimation¶
def ar_spectral_estimation():
"""Demonstrate AR spectral estimation."""
np.random.seed(42)
N = 256 # Short data record
fs = 1000
# True signal: sum of two narrow peaks
t = np.arange(N) / fs
x = (np.sin(2*np.pi*100*t) + 0.7*np.sin(2*np.pi*120*t) +
0.3*np.random.randn(N))
# Autocorrelation estimate (biased)
r_full = np.correlate(x, x, mode='full') / N
center = N - 1
max_order = 50
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
# Welch's method (for reference)
f_welch, Pxx_welch = signal.welch(x, fs=fs, nperseg=min(N, 128))
# Different AR orders
orders = [2, 4, 10, 20, 30, 50]
for ax, p in zip(axes.flat, orders):
r = r_full[center:center + p + 1]
# Levinson-Durbin
a, sigma2, refl = levinson_durbin(r, p)
# AR PSD
f_ar, Pxx_ar = ar_psd(a, sigma2, fs)
ax.semilogy(f_welch, Pxx_welch, 'b-', alpha=0.5, label='Welch')
ax.semilogy(f_ar, Pxx_ar, 'r-', linewidth=2, label=f'AR({p})')
ax.axvline(100, color='gray', linestyle='--', alpha=0.5)
ax.axvline(120, color='gray', linestyle='--', alpha=0.5)
ax.set_title(f'AR Order p = {p}')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD')
ax.set_xlim(0, 300)
ax.legend(fontsize=8)
ax.grid(True, which='both', alpha=0.3)
plt.suptitle(f'AR Spectral Estimation (N={N} samples)', fontsize=14)
plt.tight_layout()
plt.savefig('ar_spectral.png', dpi=150)
plt.close()
ar_spectral_estimation()
9.6 Burg's Method¶
Burg's method estimates AR parameters directly from the data (without first computing the autocorrelation), yielding better resolution:
def burg_method(x, order):
"""
Burg's method for AR parameter estimation.
Parameters:
x: input signal
order: AR model order
Returns:
a: AR coefficients
sigma2: prediction error power
"""
N = len(x)
# Initialize forward and backward prediction errors
ef = x.copy()
eb = x.copy()
a = np.zeros(order)
sigma2 = np.mean(np.abs(x)**2)
for m in range(order):
# Compute reflection coefficient
ef_m = ef[m+1:]
eb_m = eb[m:-1]
num = -2 * np.sum(ef_m * eb_m)
den = np.sum(ef_m**2) + np.sum(eb_m**2)
k = num / den if den > 0 else 0
# Update AR coefficients
a_new = np.zeros(order)
a_new[m] = k
for j in range(m):
a_new[j] = a[j] + k * a[m - 1 - j]
a = a_new
# Update prediction errors
ef_new = ef[m+1:] + k * eb[m:-1]
eb_new = eb[m:-1] + k * ef[m+1:]
ef[m+1:len(ef_new)+m+1] = ef_new
eb[m:len(eb_new)+m] = eb_new
# Update variance
sigma2 = sigma2 * (1 - k**2)
return a, sigma2
# Compare Burg vs Yule-Walker for a short data record
np.random.seed(42)
N = 64 # Very short!
fs = 1000
t = np.arange(N) / fs
# Two close sinusoids
x = np.sin(2*np.pi*100*t) + 0.8*np.sin(2*np.pi*115*t) + 0.3*np.random.randn(N)
p = 20 # AR order
# Yule-Walker
r_full = np.correlate(x, x, mode='full') / N
center = N - 1
r = r_full[center:center + p + 1]
a_yw, sigma2_yw, _ = levinson_durbin(r, p)
f_yw, Pxx_yw = ar_psd(a_yw, sigma2_yw, fs)
# Burg
a_burg, sigma2_burg = burg_method(x, p)
f_burg, Pxx_burg = ar_psd(a_burg, sigma2_burg, fs)
# Periodogram
f_per, Pxx_per = signal.periodogram(x, fs=fs)
fig, ax = plt.subplots(figsize=(12, 6))
ax.semilogy(f_per, Pxx_per, 'b-', alpha=0.5, label='Periodogram')
ax.semilogy(f_yw, Pxx_yw, 'g-', linewidth=2, label=f'Yule-Walker AR({p})')
ax.semilogy(f_burg, Pxx_burg, 'r-', linewidth=2, label=f'Burg AR({p})')
ax.axvline(100, color='gray', linestyle='--', alpha=0.5)
ax.axvline(115, color='gray', linestyle='--', alpha=0.5)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD')
ax.set_title(f'Burg vs Yule-Walker (N={N}, p={p})')
ax.set_xlim(0, 250)
ax.legend()
ax.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('burg_method.png', dpi=150)
plt.close()
10. Parametric Methods: MA and ARMA¶
10.1 Moving Average (MA) Model¶
An MA process of order $q$:
$$x[n] = w[n] + \sum_{k=1}^{q} b_k w[n-k]$$
PSD:
$$S_x(e^{j\omega}) = \sigma_w^2 \left|1 + \sum_{k=1}^{q} b_k e^{-j\omega k}\right|^2$$
MA models are useful for processes with spectral nulls (zeros in the spectrum).
10.2 ARMA Model¶
An ARMA($p, q$) process combines AR and MA:
$$x[n] = -\sum_{k=1}^{p} a_k x[n-k] + w[n] + \sum_{k=1}^{q} b_k w[n-k]$$
PSD:
$$S_x(e^{j\omega}) = \sigma_w^2 \frac{|B(e^{j\omega})|^2}{|A(e^{j\omega})|^2}$$
where $B(z) = 1 + b_1 z^{-1} + \cdots + b_q z^{-q}$ and $A(z) = 1 + a_1 z^{-1} + \cdots + a_p z^{-p}$.
ARMA models can represent both spectral peaks (from AR part) and spectral nulls (from MA part) with fewer parameters than a pure AR or MA model.
10.3 ARMA Parameter Estimation¶
from scipy.optimize import minimize
def arma_psd_estimate(x, p, q, fs, nfft=1024):
"""
Estimate ARMA(p,q) PSD using a two-stage method.
Stage 1: Estimate AR model of high order to approximate the signal
Stage 2: Fit ARMA(p,q) model to the AR-estimated autocorrelation
Parameters:
x: input signal
p: AR order
q: MA order
fs: sampling frequency
nfft: FFT size for PSD computation
Returns:
f: frequency vector
Pxx: PSD estimate
"""
N = len(x)
# Stage 1: High-order AR estimate
high_order = max(p + q, 30)
r_full = np.correlate(x, x, mode='full') / N
center = N - 1
r = r_full[center:center + high_order + 1]
a_high, sigma2_high, _ = levinson_durbin(r, high_order)
# Stage 2: Fit ARMA(p,q) - simplified approach using scipy
# We minimize the prediction error
def arma_cost(params):
a_params = params[:p]
b_params = params[p:p+q]
sigma = params[-1]
# Compute theoretical autocorrelation from ARMA params
# and compare with estimated autocorrelation
try:
# Simple cost: prediction error on the data
# Forward filter: A(z) * x[n] should be close to MA filtered noise
y = x.copy()
for k in range(p):
y[k+1:] += a_params[k] * x[:N-k-1]
return np.sum(y**2) / N
except Exception:
return 1e10
# Initial guess from AR model
x0 = np.zeros(p + q + 1)
if p > 0:
r_init = r_full[center:center + p + 1]
a_init, sigma_init, _ = levinson_durbin(r_init, p)
x0[:p] = a_init
x0[-1] = 1.0
# For simplicity, fall back to AR PSD as approximation
r_ar = r_full[center:center + p + 1]
a_ar, sigma2_ar, _ = levinson_durbin(r_ar, max(p, 1))
f = np.linspace(0, fs/2, nfft//2 + 1)
omega = 2 * np.pi * f / fs
A_freq = np.ones(len(f), dtype=complex)
for k in range(len(a_ar)):
A_freq += a_ar[k] * np.exp(-1j * (k+1) * omega)
Pxx = sigma2_ar / (np.abs(A_freq)**2 * fs)
return f, Pxx
# Demonstrate: signal that needs ARMA model
np.random.seed(42)
N = 1024
fs = 1000
# True ARMA(2,2) process
b_true = np.array([1.0, 0.5, 0.3]) # MA coefficients
a_true = np.array([1.0, -1.5, 0.85]) # AR coefficients
w = np.random.randn(N)
x_arma = signal.lfilter(b_true, a_true, w)
# Compare AR and ARMA estimates
fig, ax = plt.subplots(figsize=(12, 6))
# Welch (reference)
f_welch, Pxx_welch = signal.welch(x_arma, fs=fs, nperseg=256)
ax.semilogy(f_welch, Pxx_welch, 'b-', alpha=0.5, linewidth=1, label='Welch')
# AR estimates of different orders
for p in [2, 5, 10, 20]:
r = np.correlate(x_arma, x_arma, mode='full')[N-1:N-1+p+1] / N
a_est, sig2, _ = levinson_durbin(r, p)
f_ar, Pxx_ar = ar_psd(a_est, sig2, fs)
ax.semilogy(f_ar, Pxx_ar, linewidth=1.5, label=f'AR({p})')
# True ARMA PSD
f_true = np.linspace(0, fs/2, 512)
omega_t = 2 * np.pi * f_true / fs
B_freq = np.polyval(b_true[::-1], np.exp(-1j * omega_t))
A_freq = np.polyval(a_true[::-1], np.exp(-1j * omega_t))
Pxx_true = np.abs(B_freq)**2 / (np.abs(A_freq)**2 * fs)
ax.semilogy(f_true, Pxx_true, 'k--', linewidth=2, label='True ARMA(2,2)')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD')
ax.set_title('AR Approximation of ARMA Process')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('arma_comparison.png', dpi=150)
plt.close()
11. Model Order Selection¶
11.1 The Problem¶
Choosing the correct model order $p$ is critical: - Too low: underfitting, spectral peaks missed, biased estimate - Too high: overfitting, spurious peaks, high variance
11.2 Information Criteria¶
Akaike Information Criterion (AIC):
$$\text{AIC}(p) = N \ln(\hat{\sigma}_p^2) + 2p$$
Bayesian Information Criterion (BIC) (also called Schwarz criterion):
$$\text{BIC}(p) = N \ln(\hat{\sigma}_p^2) + p \ln(N)$$
where $\hat{\sigma}_p^2$ is the estimated prediction error variance for order $p$.
- AIC tends to overestimate the order
- BIC tends to give the correct order (consistent estimator)
- Both penalize complexity, but BIC penalizes more heavily for large $N$
11.3 Final Prediction Error (FPE)¶
$$\text{FPE}(p) = \hat{\sigma}_p^2 \cdot \frac{N + p + 1}{N - p - 1}$$
11.4 Implementation¶
def model_order_selection(x, max_order=50):
"""
Select AR model order using AIC, BIC, and FPE.
Parameters:
x: input signal
max_order: maximum order to test
Returns:
best_orders: dictionary of best orders by each criterion
criteria: dictionary of criterion values
"""
N = len(x)
# Autocorrelation
r_full = np.correlate(x, x, mode='full') / N
center = N - 1
orders = np.arange(1, max_order + 1)
aic_values = np.zeros(max_order)
bic_values = np.zeros(max_order)
fpe_values = np.zeros(max_order)
for i, p in enumerate(orders):
r = r_full[center:center + p + 1]
try:
a, sigma2, _ = levinson_durbin(r, p)
sigma2 = max(sigma2, 1e-15) # Prevent log of zero
except Exception:
sigma2 = 1e-15
aic_values[i] = N * np.log(sigma2) + 2 * p
bic_values[i] = N * np.log(sigma2) + p * np.log(N)
fpe_values[i] = sigma2 * (N + p + 1) / max(N - p - 1, 1)
best_aic = orders[np.argmin(aic_values)]
best_bic = orders[np.argmin(bic_values)]
best_fpe = orders[np.argmin(fpe_values)]
return {
'AIC': best_aic,
'BIC': best_bic,
'FPE': best_fpe,
}, {
'AIC': aic_values,
'BIC': bic_values,
'FPE': fpe_values,
'orders': orders,
}
# Example: known AR(4) process
np.random.seed(42)
N = 1024
fs = 1000
# True AR(4)
a_true = [1.0, -2.7607, 3.8106, -2.6535, 0.9238] # 4 poles
w = np.random.randn(N)
x_ar4 = signal.lfilter([1], a_true, w)
best_orders, criteria = model_order_selection(x_ar4, max_order=30)
print(f"True order: 4")
print(f"AIC selects: p = {best_orders['AIC']}")
print(f"BIC selects: p = {best_orders['BIC']}")
print(f"FPE selects: p = {best_orders['FPE']}")
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for ax, (name, values) in zip(axes, [('AIC', criteria['AIC']),
('BIC', criteria['BIC']),
('FPE', criteria['FPE'])]):
ax.plot(criteria['orders'], values, 'b-o', markersize=3, linewidth=1.5)
best_p = best_orders[name]
ax.axvline(best_p, color='r', linestyle='--', label=f'Best: p={best_p}')
ax.axvline(4, color='g', linestyle=':', label='True: p=4')
ax.set_xlabel('Model Order p')
ax.set_ylabel(name)
ax.set_title(f'{name} (Best: p={best_p})')
ax.legend()
ax.grid(True, alpha=0.3)
plt.suptitle('Model Order Selection for AR(4) Process', fontsize=14)
plt.tight_layout()
plt.savefig('model_order_selection.png', dpi=150)
plt.close()
12. Cross-Spectral Analysis and Coherence¶
12.1 Cross-Power Spectral Density¶
The cross-PSD between two signals $x[n]$ and $y[n]$ is:
$$S_{xy}(e^{j\omega}) = \sum_{k=-\infty}^{\infty} r_{xy}[k] e^{-j\omega k}$$
where $r_{xy}[k] = E\{x[n] y^*[n-k]\}$ is the cross-correlation.
The cross-PSD is generally complex-valued:
$$S_{xy}(e^{j\omega}) = |S_{xy}(e^{j\omega})| e^{j\phi_{xy}(\omega)}$$
12.2 Coherence¶
The magnitude-squared coherence (MSC) measures the linear relationship between two signals at each frequency:
$$C_{xy}(f) = \frac{|S_{xy}(f)|^2}{S_{xx}(f) \cdot S_{yy}(f)}$$
Properties: - $0 \leq C_{xy}(f) \leq 1$ - $C_{xy}(f) = 1$: perfectly linearly related at frequency $f$ - $C_{xy}(f) = 0$: no linear relationship at frequency $f$ - Analogous to $R^2$ (coefficient of determination) at each frequency
12.3 Phase Spectrum¶
The cross-spectral phase reveals the phase difference between corresponding frequency components:
$$\phi_{xy}(f) = \angle S_{xy}(f) = \arctan\frac{\text{Im}\{S_{xy}(f)\}}{\text{Re}\{S_{xy}(f)\}}$$
12.4 Implementation¶
def cross_spectral_analysis():
"""Demonstrate cross-spectral analysis and coherence."""
np.random.seed(42)
N = 8192
fs = 1000
# System: y is a filtered and delayed version of x, plus independent noise
t = np.arange(N) / fs
# Input: broadband noise + tone
x = np.random.randn(N) + 2*np.sin(2*np.pi*100*t)
# System: bandpass filter + delay
sos = signal.butter(4, [80, 200], btype='bandpass', fs=fs, output='sos')
delay_samples = 20
y_filtered = signal.sosfilt(sos, x)
y = np.zeros(N)
y[delay_samples:] = y_filtered[:N - delay_samples]
y += 0.5 * np.random.randn(N) # Add independent noise
# Compute cross-spectral quantities
f, Pxx = signal.welch(x, fs=fs, nperseg=512)
f, Pyy = signal.welch(y, fs=fs, nperseg=512)
f, Pxy = signal.csd(x, y, fs=fs, nperseg=512)
# Coherence
f_coh, Cxy = signal.coherence(x, y, fs=fs, nperseg=512)
# Phase
phase_xy = np.angle(Pxy)
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
# Auto-PSDs
axes[0, 0].semilogy(f, Pxx, 'b-', label='Sxx')
axes[0, 0].semilogy(f, Pyy, 'r-', label='Syy')
axes[0, 0].set_title('Auto Power Spectral Densities')
axes[0, 0].set_xlabel('Frequency (Hz)')
axes[0, 0].set_ylabel('PSD')
axes[0, 0].legend()
axes[0, 0].grid(True, which='both', alpha=0.3)
# Cross-PSD magnitude
axes[0, 1].semilogy(f, np.abs(Pxy), 'g-', linewidth=1.5)
axes[0, 1].set_title('Cross-PSD |Sxy|')
axes[0, 1].set_xlabel('Frequency (Hz)')
axes[0, 1].set_ylabel('|Sxy|')
axes[0, 1].grid(True, which='both', alpha=0.3)
# Coherence
axes[0, 2].plot(f_coh, Cxy, 'm-', linewidth=1.5)
axes[0, 2].axhline(1.0, color='gray', linestyle='--', alpha=0.5)
axes[0, 2].axhspan(80/500, 200/500, alpha=0.1, color='green',
transform=axes[0, 2].get_xaxis_transform())
axes[0, 2].set_title('Coherence Cxy(f)')
axes[0, 2].set_xlabel('Frequency (Hz)')
axes[0, 2].set_ylabel('Coherence')
axes[0, 2].set_ylim(0, 1.1)
axes[0, 2].grid(True, alpha=0.3)
# Phase spectrum (only where coherence is significant)
significant = Cxy > 0.3
f_sig = f_coh[significant]
phase_sig = phase_xy[significant]
axes[1, 0].plot(f_sig, phase_sig * 180 / np.pi, 'k.', markersize=3)
axes[1, 0].set_title('Cross-Spectral Phase (where C > 0.3)')
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('Phase (degrees)')
axes[1, 0].grid(True, alpha=0.3)
# Transfer function estimate
H_est = Pxy / (Pxx + 1e-15)
w_sys, H_sys = signal.sosfreqz(sos, worN=len(f), fs=fs)
axes[1, 1].plot(f, 20*np.log10(np.abs(H_est) + 1e-15), 'b-', alpha=0.7,
label='Estimated |H(f)|')
axes[1, 1].plot(w_sys, 20*np.log10(np.abs(H_sys) + 1e-15), 'r--',
linewidth=2, label='True |H(f)|')
axes[1, 1].set_title('Transfer Function Estimate')
axes[1, 1].set_xlabel('Frequency (Hz)')
axes[1, 1].set_ylabel('Magnitude (dB)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
# Time-domain
axes[1, 2].plot(t[:500]*1000, x[:500], 'b-', alpha=0.5, label='x')
axes[1, 2].plot(t[:500]*1000, y[:500], 'r-', alpha=0.5, label='y')
axes[1, 2].set_title('Time Domain')
axes[1, 2].set_xlabel('Time (ms)')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)
plt.suptitle('Cross-Spectral Analysis and Coherence', fontsize=14)
plt.tight_layout()
plt.savefig('cross_spectral.png', dpi=150)
plt.close()
cross_spectral_analysis()
13. Spectrogram¶
13.1 Time-Frequency Analysis¶
The spectrogram shows how the frequency content of a signal changes over time. It is computed as a sequence of short-time Fourier transforms (STFT):
$$S(t, f) = \left| \text{STFT}\{x[n]\}(t, f) \right|^2 = \left| \sum_{m=-\infty}^{\infty} x[m] w[m - t] e^{-j2\pi f m} \right|^2$$
The spectrogram is essentially a time-varying PSD obtained by sliding a window along the signal.
13.2 Parameters¶
- Window length (nperseg): determines frequency resolution ($\Delta f = f_s / \text{nperseg}$)
- Overlap: determines time resolution (more overlap = finer time steps)
- Window type: controls spectral leakage (same as in Welch's method)
13.3 Time-Frequency Resolution¶
The uncertainty principle limits simultaneous time and frequency resolution:
$$\Delta t \cdot \Delta f \geq \frac{1}{4\pi}$$
- Long window: good frequency resolution, poor time resolution
- Short window: good time resolution, poor frequency resolution
13.4 Implementation¶
def spectrogram_analysis():
"""Demonstrate spectrogram with different parameters."""
fs = 8000
duration = 2.0
t = np.arange(0, duration, 1/fs)
N = len(t)
# Create a time-varying signal
# Chirp (frequency sweep from 100 to 3000 Hz)
x_chirp = signal.chirp(t, f0=100, t1=duration, f1=3000, method='linear')
# Add some tone bursts
x_bursts = np.zeros(N)
burst_start = int(0.5 * fs)
burst_end = int(1.0 * fs)
x_bursts[burst_start:burst_end] = np.sin(2*np.pi*1500*t[burst_start:burst_end])
burst_start2 = int(1.2 * fs)
burst_end2 = int(1.5 * fs)
x_bursts[burst_start2:burst_end2] = 0.5*np.sin(2*np.pi*2500*t[burst_start2:burst_end2])
x = x_chirp + x_bursts + 0.1*np.random.randn(N)
fig, axes = plt.subplots(3, 2, figsize=(16, 14))
# Time domain
axes[0, 0].plot(t, x, 'b-', linewidth=0.5)
axes[0, 0].set_xlabel('Time (s)')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_title('Signal (Chirp + Tone Bursts + Noise)')
axes[0, 0].grid(True, alpha=0.3)
# Spectrum (overall)
f_welch, Pxx = signal.welch(x, fs=fs, nperseg=1024)
axes[0, 1].semilogy(f_welch, Pxx, 'b-', linewidth=1.5)
axes[0, 1].set_xlabel('Frequency (Hz)')
axes[0, 1].set_ylabel('PSD')
axes[0, 1].set_title('Overall PSD (Welch)')
axes[0, 1].grid(True, which='both', alpha=0.3)
# Spectrograms with different window sizes
for idx, (nperseg, title_prefix) in enumerate([(64, 'Short'), (256, 'Medium'),
(1024, 'Long')]):
row = 1 + idx // 2
col = idx % 2
if idx == 2:
row, col = 2, 0
f_spec, t_spec, Sxx = signal.spectrogram(x, fs=fs, nperseg=nperseg,
noverlap=nperseg*3//4,
window='hann')
im = axes[row, col].pcolormesh(t_spec, f_spec,
10*np.log10(Sxx + 1e-15),
shading='gouraud', cmap='viridis',
vmin=-60, vmax=0)
axes[row, col].set_xlabel('Time (s)')
axes[row, col].set_ylabel('Frequency (Hz)')
axes[row, col].set_title(f'{title_prefix} Window (nperseg={nperseg}, '
f'ฮf={fs/nperseg:.0f} Hz, '
f'ฮt={nperseg/fs*1000:.1f} ms)')
plt.colorbar(im, ax=axes[row, col], label='PSD (dB)')
# matplotlib specgram for comparison
axes[2, 1].specgram(x, NFFT=256, Fs=fs, noverlap=192,
cmap='viridis', vmin=-60, vmax=0)
axes[2, 1].set_xlabel('Time (s)')
axes[2, 1].set_ylabel('Frequency (Hz)')
axes[2, 1].set_title('matplotlib specgram (NFFT=256)')
plt.suptitle('Spectrogram: Time-Frequency Analysis', fontsize=14)
plt.tight_layout()
plt.savefig('spectrogram_analysis.png', dpi=150)
plt.close()
spectrogram_analysis()
13.5 Interactive Spectrogram Example¶
def music_like_spectrogram():
"""Create a spectrogram of a music-like signal."""
fs = 44100
duration = 3.0
t = np.arange(0, duration, 1/fs)
# Simulate a simple melody
notes_hz = {
'C4': 261.63, 'D4': 293.66, 'E4': 329.63, 'F4': 349.23,
'G4': 392.00, 'A4': 440.00, 'B4': 493.88, 'C5': 523.25,
}
melody = ['C4', 'E4', 'G4', 'C5', 'G4', 'E4', 'C4', 'D4',
'F4', 'A4', 'C5', 'A4']
note_duration = duration / len(melody)
x = np.zeros(len(t))
for i, note in enumerate(melody):
start_idx = int(i * note_duration * fs)
end_idx = int((i + 1) * note_duration * fs)
if end_idx > len(t):
end_idx = len(t)
t_note = t[start_idx:end_idx] - t[start_idx]
f0 = notes_hz[note]
# Envelope (ADSR-like)
env_len = end_idx - start_idx
attack = int(0.05 * fs)
decay = int(0.1 * fs)
release = int(0.05 * fs)
envelope = np.ones(env_len)
envelope[:attack] = np.linspace(0, 1, attack)
envelope[attack:attack+decay] = np.linspace(1, 0.7, decay)
envelope[-release:] = np.linspace(0.7, 0, release)
# Add harmonics
for harmonic in [1, 2, 3, 4]:
x[start_idx:end_idx] += (envelope *
np.sin(2*np.pi*f0*harmonic*t_note) / harmonic)
x += 0.02 * np.random.randn(len(x)) # Add slight noise
# Spectrogram
fig, axes = plt.subplots(2, 1, figsize=(14, 8))
axes[0].plot(t, x, 'b-', linewidth=0.3)
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Music-like Signal (Time Domain)')
axes[0].grid(True, alpha=0.3)
f_spec, t_spec, Sxx = signal.spectrogram(x, fs=fs, nperseg=4096,
noverlap=3072, window='hann')
im = axes[1].pcolormesh(t_spec, f_spec, 10*np.log10(Sxx + 1e-15),
shading='gouraud', cmap='magma',
vmin=-80, vmax=-10)
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Frequency (Hz)')
axes[1].set_ylim(0, 3000)
axes[1].set_title('Spectrogram')
plt.colorbar(im, ax=axes[1], label='PSD (dB)')
plt.tight_layout()
plt.savefig('music_spectrogram.png', dpi=150)
plt.close()
music_like_spectrogram()
14. Python Implementation¶
14.1 Comprehensive Spectral Analysis Toolkit¶
class SpectralAnalyzer:
"""Complete spectral analysis toolkit."""
def __init__(self, x, fs):
self.x = np.array(x, dtype=float)
self.fs = fs
self.N = len(x)
def periodogram(self, nfft=None, window='boxcar'):
"""Compute periodogram."""
return signal.periodogram(self.x, fs=self.fs, nfft=nfft, window=window)
def welch(self, nperseg=256, noverlap=None, window='hann'):
"""Compute Welch's PSD estimate."""
if noverlap is None:
noverlap = nperseg // 2
return signal.welch(self.x, fs=self.fs, nperseg=nperseg,
noverlap=noverlap, window=window)
def ar_psd(self, order=None, method='burg', nfft=1024):
"""Compute AR model PSD."""
if order is None:
# Auto-select using BIC
best_orders, _ = model_order_selection(self.x, max_order=50)
order = best_orders['BIC']
print(f"Auto-selected AR order: p = {order} (BIC)")
if method == 'burg':
a, sigma2 = burg_method(self.x, order)
else:
r_full = np.correlate(self.x, self.x, mode='full') / self.N
center = self.N - 1
r = r_full[center:center + order + 1]
a, sigma2, _ = levinson_durbin(r, order)
f = np.linspace(0, self.fs/2, nfft//2 + 1)
omega = 2 * np.pi * f / self.fs
A_freq = np.ones(len(f), dtype=complex)
for k in range(len(a)):
A_freq += a[k] * np.exp(-1j * (k+1) * omega)
Pxx = sigma2 / (np.abs(A_freq)**2 * self.fs)
return f, Pxx
def spectrogram(self, nperseg=256, noverlap=None, window='hann'):
"""Compute spectrogram."""
if noverlap is None:
noverlap = nperseg * 3 // 4
return signal.spectrogram(self.x, fs=self.fs, nperseg=nperseg,
noverlap=noverlap, window=window)
def comprehensive_analysis(self, save_path='comprehensive_spectral.png'):
"""Perform and plot comprehensive spectral analysis."""
fig, axes = plt.subplots(3, 2, figsize=(16, 14))
# Time domain
t = np.arange(self.N) / self.fs
axes[0, 0].plot(t, self.x, 'b-', linewidth=0.5)
axes[0, 0].set_title('Time Domain')
axes[0, 0].set_xlabel('Time (s)')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].grid(True, alpha=0.3)
# Periodogram
f_per, Pxx_per = self.periodogram()
axes[0, 1].semilogy(f_per, Pxx_per, 'b-', alpha=0.5, linewidth=0.5)
axes[0, 1].set_title('Periodogram')
axes[0, 1].set_xlabel('Frequency (Hz)')
axes[0, 1].set_ylabel('PSD')
axes[0, 1].grid(True, which='both', alpha=0.3)
# Welch
f_w, Pxx_w = self.welch(nperseg=min(self.N // 4, 512))
axes[1, 0].semilogy(f_w, Pxx_w, 'g-', linewidth=1.5)
axes[1, 0].set_title("Welch's Method")
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('PSD')
axes[1, 0].grid(True, which='both', alpha=0.3)
# AR PSD
f_ar, Pxx_ar = self.ar_psd()
axes[1, 1].semilogy(f_ar, Pxx_ar, 'r-', linewidth=1.5,
label='AR (Burg)')
axes[1, 1].semilogy(f_w, Pxx_w, 'g--', alpha=0.5, label='Welch')
axes[1, 1].set_title('AR Spectral Estimate')
axes[1, 1].set_xlabel('Frequency (Hz)')
axes[1, 1].set_ylabel('PSD')
axes[1, 1].legend()
axes[1, 1].grid(True, which='both', alpha=0.3)
# Spectrogram
f_spec, t_spec, Sxx = self.spectrogram(nperseg=min(self.N // 8, 256))
im = axes[2, 0].pcolormesh(t_spec, f_spec,
10*np.log10(Sxx + 1e-15),
shading='gouraud', cmap='viridis')
axes[2, 0].set_title('Spectrogram')
axes[2, 0].set_xlabel('Time (s)')
axes[2, 0].set_ylabel('Frequency (Hz)')
plt.colorbar(im, ax=axes[2, 0], label='PSD (dB)')
# Autocorrelation
max_lag = min(500, self.N // 4)
r = np.correlate(self.x, self.x, mode='full') / self.N
center = self.N - 1
lags = np.arange(-max_lag, max_lag + 1)
axes[2, 1].plot(lags / self.fs * 1000,
r[center-max_lag:center+max_lag+1], 'b-', linewidth=1)
axes[2, 1].set_title('Autocorrelation')
axes[2, 1].set_xlabel('Lag (ms)')
axes[2, 1].set_ylabel('r_xx[k]')
axes[2, 1].grid(True, alpha=0.3)
plt.suptitle('Comprehensive Spectral Analysis', fontsize=14)
plt.tight_layout()
plt.savefig(save_path, dpi=150)
plt.close()
return fig
# Example usage
np.random.seed(42)
fs = 1000
N = 4096
t = np.arange(N) / fs
# Complex signal: tones + AR noise + chirp
x = (np.sin(2*np.pi*50*t) +
0.5*np.sin(2*np.pi*120*t) +
signal.chirp(t, 200, t[-1], 300) * 0.3 +
0.3*np.random.randn(N))
analyzer = SpectralAnalyzer(x, fs)
analyzer.comprehensive_analysis()
14.2 Real-World Example: Vibration Analysis¶
def vibration_analysis_demo():
"""Simulate vibration analysis of a rotating machine."""
fs = 10000 # 10 kHz sampling
duration = 5.0
t = np.arange(0, duration, 1/fs)
N = len(t)
# Fundamental rotation frequency
f_rot = 30 # Hz (1800 RPM)
# Normal vibration: fundamental + harmonics
x_normal = (0.5 * np.sin(2*np.pi*f_rot*t) +
0.2 * np.sin(2*np.pi*2*f_rot*t) +
0.1 * np.sin(2*np.pi*3*f_rot*t) +
0.05 * np.random.randn(N))
# Faulty vibration: bearing defect frequency
f_bearing = 4.7 * f_rot # Ball pass frequency (outer race)
x_faulty = (0.5 * np.sin(2*np.pi*f_rot*t) +
0.2 * np.sin(2*np.pi*2*f_rot*t) +
0.1 * np.sin(2*np.pi*3*f_rot*t) +
0.3 * np.sin(2*np.pi*f_bearing*t) + # Bearing defect
0.15 * np.sin(2*np.pi*2*f_bearing*t) +
0.05 * np.random.randn(N))
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# PSDs
f_n, Pxx_n = signal.welch(x_normal, fs=fs, nperseg=2048)
f_f, Pxx_f = signal.welch(x_faulty, fs=fs, nperseg=2048)
axes[0, 0].semilogy(f_n, Pxx_n, 'b-', linewidth=1.5)
axes[0, 0].set_title('Normal Machine PSD')
axes[0, 0].set_xlabel('Frequency (Hz)')
axes[0, 0].set_ylabel('PSD')
axes[0, 0].set_xlim(0, 500)
for h in range(1, 4):
axes[0, 0].axvline(h*f_rot, color='g', linestyle='--', alpha=0.5)
axes[0, 0].grid(True, which='both', alpha=0.3)
axes[0, 1].semilogy(f_f, Pxx_f, 'r-', linewidth=1.5)
axes[0, 1].set_title('Faulty Machine PSD (Bearing Defect)')
axes[0, 1].set_xlabel('Frequency (Hz)')
axes[0, 1].set_ylabel('PSD')
axes[0, 1].set_xlim(0, 500)
for h in range(1, 4):
axes[0, 1].axvline(h*f_rot, color='g', linestyle='--', alpha=0.5)
axes[0, 1].axvline(f_bearing, color='red', linestyle=':', alpha=0.7,
label=f'BPFO = {f_bearing:.0f} Hz')
axes[0, 1].axvline(2*f_bearing, color='red', linestyle=':', alpha=0.5)
axes[0, 1].legend()
axes[0, 1].grid(True, which='both', alpha=0.3)
# Overlay comparison
axes[1, 0].semilogy(f_n, Pxx_n, 'b-', linewidth=1.5, label='Normal')
axes[1, 0].semilogy(f_f, Pxx_f, 'r-', linewidth=1.5, label='Faulty')
axes[1, 0].set_title('Comparison')
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('PSD')
axes[1, 0].set_xlim(0, 500)
axes[1, 0].legend()
axes[1, 0].grid(True, which='both', alpha=0.3)
# Spectrogram (faulty - developing fault)
# Simulate fault developing over time
x_developing = np.zeros(N)
fault_amplitude = np.linspace(0, 0.4, N) # Gradually increasing
x_developing = (0.5 * np.sin(2*np.pi*f_rot*t) +
0.2 * np.sin(2*np.pi*2*f_rot*t) +
fault_amplitude * np.sin(2*np.pi*f_bearing*t) +
0.05 * np.random.randn(N))
f_spec, t_spec, Sxx = signal.spectrogram(x_developing, fs=fs,
nperseg=2048, noverlap=1536)
im = axes[1, 1].pcolormesh(t_spec, f_spec, 10*np.log10(Sxx + 1e-15),
shading='gouraud', cmap='inferno',
vmin=-60, vmax=-10)
axes[1, 1].set_title('Developing Fault (Spectrogram)')
axes[1, 1].set_xlabel('Time (s)')
axes[1, 1].set_ylabel('Frequency (Hz)')
axes[1, 1].set_ylim(0, 500)
plt.colorbar(im, ax=axes[1, 1], label='PSD (dB)')
plt.suptitle('Vibration Spectral Analysis: Machine Condition Monitoring', fontsize=14)
plt.tight_layout()
plt.savefig('vibration_analysis.png', dpi=150)
plt.close()
vibration_analysis_demo()
15. Exercises¶
Exercise 1: Periodogram vs Welch¶
Generate a signal consisting of three sinusoids at 100, 150, and 200 Hz with amplitudes 1.0, 0.5, and 0.1 respectively, in additive white Gaussian noise (SNR = 10 dB). Use $f_s = 1000$ Hz and $N = 4096$ samples.
(a) Compute and plot the periodogram. Can you identify all three tones?
(b) Apply Welch's method with nperseg = 256 and 512. Which setting resolves the tones better?
(c) How does the choice of window function (boxcar, hann, blackmanharris) affect the ability to detect the weakest tone (200 Hz)?
(d) Compute the variance of each PSD estimate by running 100 Monte Carlo trials. Plot the variance as a function of frequency for each method.
Exercise 2: AR Model Identification¶
Generate an AR(4) process with known coefficients: $x[n] + 1.5x[n-1] - 0.75x[n-2] + 0.2x[n-3] - 0.05x[n-4] = w[n]$
(a) Generate $N = 512$ samples and estimate the AR parameters using the Yule-Walker method for orders $p = 2, 4, 6, 8, 12$.
(b) For each order, compute and plot the estimated PSD alongside the true PSD.
(c) Use AIC and BIC to select the optimal order. Do they agree?
(d) Repeat with $N = 64$ samples. How does the reduced data length affect the results?
Exercise 3: Cross-Spectral Analysis¶
Two sensors measure vibrations at different locations on a structure. The system transfer function between them is a 4th-order bandpass filter centered at 150 Hz.
(a) Simulate the input (broadband random excitation) and output (filtered + independent noise).
(b) Estimate the coherence function. At which frequencies is the coherence highest?
(c) Estimate the transfer function magnitude and phase from the cross-spectral density.
(d) Compute the group delay from the phase spectrum. Does it match the known filter group delay?
Exercise 4: Spectrogram of Bat Echolocation¶
Create a synthetic bat echolocation signal consisting of: - A series of 5 chirps, each sweeping from 100 kHz to 50 kHz in 5 ms - Inter-pulse interval of 50 ms - Additive noise at -20 dB
(a) Plot the spectrogram using a Hanning window. Experiment with nperseg = 64, 256, 1024 to find the best compromise.
(b) Measure the instantaneous frequency from the spectrogram peaks. Compare with the known chirp rate.
(c) Add an echo (attenuated, delayed copy). Can you detect the echo in the spectrogram? What is the minimum detectable delay?
Exercise 5: Resolution Limits¶
Two sinusoids at frequencies $f_1$ and $f_2 = f_1 + \Delta f$ with equal amplitude in noise.
(a) For $f_1 = 100$ Hz, $f_s = 1000$ Hz, $N = 256$: what is the minimum $\Delta f$ for which the periodogram can resolve the two tones? Verify experimentally.
(b) Compare the resolution of the periodogram, Welch's method ($\text{nperseg} = 128$), and AR Burg's method ($p = 20$) for $\Delta f = 5$ Hz.
(c) How does SNR affect resolution? Test with SNR = 0, 10, 20, 40 dB and plot the minimum resolvable $\Delta f$ for each method.
Exercise 6: Real-World Application: EEG Analysis¶
Simulate a simplified EEG signal containing: - Alpha band activity (8-13 Hz) during "eyes closed" - Beta band activity (13-30 Hz) during "eyes open" - Transition from eyes-closed to eyes-open at $t = 5$ s
(a) Create a 10-second signal at $f_s = 256$ Hz that transitions between the two states.
(b) Compute the spectrogram and identify the transition point.
(c) Compute the band power (total PSD integrated over a frequency band) for alpha and beta bands as a function of time. Plot these as time series.
(d) Compute the alpha/beta power ratio. How clearly can you detect the state transition?
Exercise 7: Burg vs Welch for Short Data¶
Compare Burg's AR method with Welch's method for increasingly short data records.
(a) Generate a signal with two tones at 100 Hz and 108 Hz ($f_s = 1000$ Hz). Vary $N$ from 32 to 4096 (powers of 2).
(b) For each $N$, compute both Welch (nperseg = N/2) and Burg (order = 20) PSD estimates.
(c) Define a "resolution success" criterion: the two peaks are resolved if there is a notch of at least 3 dB between them. Determine the minimum $N$ for each method.
(d) Plot the resolution success rate (over 100 Monte Carlo trials) as a function of $N$ for both methods.
References¶
- Oppenheim, A. V., & Schafer, R. W. (2010). Discrete-Time Signal Processing (3rd ed.). Pearson. Chapter 10.
- Kay, S. M. (1988). Modern Spectral Estimation: Theory and Application. Prentice Hall. [The classic reference for parametric methods]
- Stoica, P., & Moses, R. L. (2005). Spectral Analysis of Signals. Prentice Hall. [Excellent modern treatment]
- Welch, P. D. (1967). "The Use of Fast Fourier Transform for the Estimation of Power Spectra." IEEE Trans. Audio Electroacoustics, 15(2), 70-73.
- Burg, J. P. (1975). "Maximum Entropy Spectral Analysis." PhD Dissertation, Stanford University.
- Marple, S. L. (1987). Digital Spectral Analysis with Applications. Prentice Hall.
- SciPy Documentation -- Spectral analysis functions: https://docs.scipy.org/doc/scipy/reference/signal.html#spectral-analysis
Navigation¶
- Previous: 11. Multirate Signal Processing
- Next: 13. Adaptive Filters
- Back to Overview