Exercises

ex-sp-ch07-01

Easy

Generate a signal $x(t) = 3\cos(2\pi \cdot 80\,t) + \sin(2\pi \cdot 200\,t)$ sampled at $f_s = 1000$ Hz for 1 second. Compute the FFT, plot the magnitude spectrum for positive frequencies, and verify Parseval's theorem: $\sum|x[n]|^2 = \frac{1}{N}\sum|X[k]|^2$ .

Show Hint

Use np.fft.fft and np.fft.fftfreq.

Parseval check: np.isclose(np.sum(np.abs(x)**2), np.sum(np.abs(X)**2)/N).

Solution

Implementation

import numpy as np

fs = 1000
t = np.arange(0, 1, 1/fs)
x = 3*np.cos(2*np.pi*80*t) + np.sin(2*np.pi*200*t)

X = np.fft.fft(x)
N = len(x)
freqs = np.fft.fftfreq(N, 1/fs)

# Parseval's theorem
energy_time = np.sum(np.abs(x)**2)
energy_freq = np.sum(np.abs(X)**2) / N
print(f"Time energy: {energy_time:.4f}")
print(f"Freq energy: {energy_freq:.4f}")
print(f"Match: {np.isclose(energy_time, energy_freq)}")

ex-sp-ch07-02

Easy

Construct the $16 \times 16$ DFT matrix $\mathbf{F}_{16}$ and verify: (a) $\mathbf{F}_{16}\mathbf{x}$ matches np.fft.fft(x) for a random vector, (b) $\frac{1}{\sqrt{16}}\mathbf{F}_{16}$ is unitary.

Show Hint

Use np.exp(-2j * np.pi * k * n / N) with outer product indexing.

Solution

Implementation

import numpy as np

N = 16
n = np.arange(N)
k = n.reshape(-1, 1)
F = np.exp(-2j * np.pi * k * n / N)

x = np.random.randn(N)
print(f"FFT match: {np.allclose(F @ x, np.fft.fft(x))}")
F_norm = F / np.sqrt(N)
print(f"Unitary: {np.allclose(F_norm.conj().T @ F_norm, np.eye(N))}")

ex-sp-ch07-03

Easy

Design a 5th-order Butterworth lowpass filter with cutoff 150 Hz for a signal sampled at 1000 Hz. Apply it to $x(t) = \sin(2\pi\cdot50\,t) + 0.5\sin(2\pi\cdot300\,t)$ and verify that the 300 Hz component is removed.

Show Hint

Use scipy.signal.butter with output="sos" and sosfilt.

Solution

Implementation

from scipy.signal import butter, sosfilt
import numpy as np

fs = 1000
sos = butter(5, 150, btype='low', fs=fs, output='sos')

t = np.arange(0, 1, 1/fs)
x = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*300*t)
y = sosfilt(sos, x)

# Check: FFT of filtered signal should have no 300 Hz peak
Y = np.fft.fft(y)
freqs = np.fft.fftfreq(len(y), 1/fs)
idx_300 = np.argmin(np.abs(freqs - 300))
print(f"300 Hz power ratio: {np.abs(Y[idx_300])/np.max(np.abs(Y)):.4f}")

ex-sp-ch07-04

Easy

Apply a Hann window to a 256-sample sinusoid at 100 Hz (sampled at 1024 Hz) and compare the FFT with and without the window. Measure the sidelobe level in dB for both cases.

Show Hint

Sidelobe level: find the maximum magnitude away from the main lobe.

Solution

Implementation

from scipy.signal import windows
import numpy as np

fs = 1024
N = 256
t = np.arange(N) / fs
x = np.sin(2*np.pi*100*t)

X_rect = np.fft.fft(x, n=4096)
X_hann = np.fft.fft(windows.hann(N) * x, n=4096)

mag_rect = 20*np.log10(np.abs(X_rect[:2048]) / np.max(np.abs(X_rect)))
mag_hann = 20*np.log10(np.abs(X_hann[:2048]) / np.max(np.abs(X_hann)))

print(f"Rect peak sidelobe: {np.sort(mag_rect)[-2]:.1f} dB")
print(f"Hann peak sidelobe: {np.sort(mag_hann)[-2]:.1f} dB")

ex-sp-ch07-05

Easy

Resample a 1-second sinusoid at 50 Hz from $f_s = 500$ Hz to $f_s = 800$ Hz using both resample and resample_poly. Compute the maximum error between the two results.

Show Hint

For resample_poly, simplify 800/500 = 8/5.

Solution

Implementation

from scipy.signal import resample, resample_poly
import numpy as np

fs_in = 500
t = np.arange(0, 1, 1/fs_in)
x = np.sin(2*np.pi*50*t)
n_out = int(len(x) * 800 / 500)

y_fft = resample(x, n_out)
y_poly = resample_poly(x, 8, 5)
min_len = min(len(y_fft), len(y_poly))
print(f"Max difference: {np.max(np.abs(y_fft[:min_len] - y_poly[:min_len])):.6f}")

ex-sp-ch07-06

Medium

Implement linear convolution via the FFT: zero-pad $x[n]$ (length $N$ ) and $h[n]$ (length $M$ ) to length $\ge N+M-1$ , multiply their FFTs, and take the inverse FFT. Verify the result matches np.convolve.

Show Hint

The FFT-based approach computes circular convolution, which equals linear convolution when both sequences are zero-padded.

Solution

Implementation

import numpy as np

x = np.random.randn(100)
h = np.random.randn(30)
L = len(x) + len(h) - 1

Y = np.fft.ifft(np.fft.fft(x, n=L) * np.fft.fft(h, n=L))
y_direct = np.convolve(x, h)
print(f"Match: {np.allclose(Y.real, y_direct)}")

ex-sp-ch07-07

Medium

Design an equiripple bandpass FIR filter using scipy.signal.remez with passband [1000, 2000] Hz, transition bands of 200 Hz on each side, and sampling rate 8000 Hz. Use 101 taps. Plot the magnitude response and verify the passband ripple and stopband attenuation.

Show Hint

Bands: [0, 800, 1000, 2000, 2200, 4000] with desired [0, 1, 0].

Solution

Implementation

from scipy.signal import remez, freqz
import numpy as np

fs = 8000
b = remez(101, bands=[0, 800, 1000, 2000, 2200, fs/2],
          desired=[0, 1, 0], fs=fs)
w, H = freqz(b, worN=4096, fs=fs)
mag_db = 20*np.log10(np.abs(H) + 1e-15)

passband = (w >= 1000) & (w <= 2000)
stopband = (w <= 800) | (w >= 2200)
print(f"Passband ripple: {np.ptp(mag_db[passband]):.2f} dB")
print(f"Stopband atten: {np.max(mag_db[stopband]):.1f} dB")

ex-sp-ch07-08

Medium

Compare the Welch PSD estimate with the periodogram for a signal containing three sinusoids (50, 120, 200 Hz) in white noise at SNR = 0 dB. Use segment lengths of 256, 512, and 1024 samples with 50% overlap. Which segment length gives the best trade-off between variance and resolution?

Show Hint

Generate at least 10 seconds of data for meaningful averaging.

Solution

Implementation

from scipy.signal import welch, periodogram
import numpy as np

fs = 1000
t = np.arange(0, 10, 1/fs)
x = (np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t)
     + np.sin(2*np.pi*200*t))
noise_power = np.var(x)
x += np.sqrt(noise_power) * np.random.randn(len(t))

for nperseg in [256, 512, 1024]:
    f, Pxx = welch(x, fs=fs, nperseg=nperseg, noverlap=nperseg//2)
    print(f"nperseg={nperseg}: df={f[1]:.2f} Hz, "
          f"peak SNR={10*np.log10(np.max(Pxx)/np.median(Pxx)):.1f} dB")

ex-sp-ch07-09

Medium

Implement a matched filter for a 13-chip Barker code $[+1, +1, +1, +1, +1, -1, -1, +1, +1, -1, +1, -1, +1]$ . Embed the code at a random position in a noise signal (SNR = $-5$ dB) and show that the matched filter detects the correct position.

Show Hint

The matched filter is correlate(received, barker, mode="full").

The Barker code has peak sidelobe ratio of $1/N$ .

Solution

Implementation

from scipy.signal import correlate
import numpy as np

barker13 = np.array([1,1,1,1,1,-1,-1,1,1,-1,1,-1,1], dtype=float)
N = 200
true_pos = 87
received = np.zeros(N)
received[true_pos:true_pos+13] = barker13
# SNR = -5 dB
noise_std = np.sqrt(13 * 10**(5/10))
received += noise_std * np.random.randn(N)

mf_out = correlate(received, barker13, mode='full')
detected_pos = np.argmax(np.abs(mf_out)) - 12
print(f"True: {true_pos}, Detected: {detected_pos}")
print(f"Match: {true_pos == detected_pos}")

ex-sp-ch07-10

Medium

Decompose a noisy signal using the DWT with the 'db4' wavelet at 5 levels. Apply soft thresholding with the universal threshold and reconstruct the denoised signal. Report the MSE improvement.

Show Hint

Use pywt.wavedec, pywt.threshold, and pywt.waverec.

Estimate $\\sigma$ from the MAD of the finest detail coefficients.

Solution

Implementation

import pywt
import numpy as np

np.random.seed(42)
N = 1024
t = np.linspace(0, 1, N)
x_clean = np.sin(4*np.pi*t) * np.exp(-3*t)
sigma = 0.2
x_noisy = x_clean + sigma * np.random.randn(N)

coeffs = pywt.wavedec(x_noisy, 'db4', level=5)
sigma_est = np.median(np.abs(coeffs[-1])) / 0.6745
threshold = sigma_est * np.sqrt(2*np.log(N))

coeffs_t = [coeffs[0]] + [pywt.threshold(c, threshold, 'soft')
                            for c in coeffs[1:]]
x_den = pywt.waverec(coeffs_t, 'db4')[:N]

mse_noisy = np.mean((x_clean - x_noisy)**2)
mse_den = np.mean((x_clean - x_den)**2)
print(f"MSE noisy: {mse_noisy:.4f}, denoised: {mse_den:.4f}")
print(f"Improvement: {10*np.log10(mse_noisy/mse_den):.1f} dB")

ex-sp-ch07-11

Medium

Compute the spectrogram of a signal composed of two chirps: one sweeping 100-400 Hz (0-1 s) and another sweeping 400-100 Hz (1-2 s). Use scipy.signal.spectrogram and experiment with different nperseg values (64, 256, 1024) to observe the time-frequency resolution trade-off.

Show Hint

Use scipy.signal.chirp for each segment and concatenate.

Solution

Implementation

from scipy.signal import chirp, spectrogram
import numpy as np

fs = 2000
t1 = np.arange(0, 1, 1/fs)
t2 = np.arange(0, 1, 1/fs)
x = np.concatenate([chirp(t1, 100, 1, 400), chirp(t2, 400, 1, 100)])

for nperseg in [64, 256, 1024]:
    f, t, Sxx = spectrogram(x, fs=fs, nperseg=nperseg, noverlap=nperseg-1)
    print(f"nperseg={nperseg}: {len(f)} freq bins, {len(t)} time bins")

ex-sp-ch07-12

Hard

Implement the overlap-add method for efficient convolution of a very long signal ( $N = 10^6$ ) with a short filter ( $M = 128$ ). Break the signal into blocks of length $L$ , zero-pad each block to $L+M-1$ , convolve via FFT, and accumulate the overlapping tails. Verify the result matches fftconvolve and benchmark both approaches.

Show Hint

Choose $L$ as a power of 2 for FFT efficiency.

Pre-compute the FFT of the zero-padded filter.

Solution

Implementation

from scipy.signal import fftconvolve
import numpy as np
import time

def overlap_add(x, h, L=1024):
    M = len(h)
    N_fft = L + M - 1
    H = np.fft.fft(h, N_fft)
    y = np.zeros(len(x) + M - 1)
    for i in range(0, len(x), L):
        block = x[i:i+L]
        Y_block = np.fft.ifft(np.fft.fft(block, N_fft) * H).real
        y[i:i+N_fft] += Y_block[:len(Y_block)]
    return y

x = np.random.randn(1_000_000)
h = np.random.randn(128)

t0 = time.perf_counter()
y_oa = overlap_add(x, h)
t_oa = time.perf_counter() - t0

t0 = time.perf_counter()
y_ref = fftconvolve(x, h)
t_ref = time.perf_counter() - t0

print(f"Match: {np.allclose(y_oa[:len(y_ref)], y_ref, atol=1e-10)}")
print(f"OA: {t_oa:.3f}s, fftconvolve: {t_ref:.3f}s")

ex-sp-ch07-13

Hard

Design a Chebyshev Type I bandpass filter (order 6, 1 dB passband ripple) for the band [500, 1500] Hz at $f_s = 8000$ Hz. Compare its magnitude response, group delay, and step response with an equivalent-order Butterworth bandpass filter.

Show Hint

Use scipy.signal.cheby1 with btype="band" and compare with butter.

Compute group delay with scipy.signal.group_delay.

Solution

Implementation

from scipy.signal import cheby1, butter, sosfreqz, sosfilt, group_delay
import numpy as np

fs = 8000
sos_cheb = cheby1(6, 1, [500, 1500], btype='band', fs=fs, output='sos')
sos_butt = butter(6, [500, 1500], btype='band', fs=fs, output='sos')

w_c, H_c = sosfreqz(sos_cheb, worN=4096, fs=fs)
w_b, H_b = sosfreqz(sos_butt, worN=4096, fs=fs)

print(f"Cheby passband ripple: {np.ptp(20*np.log10(np.abs(H_c[(w_c>=500)&(w_c<=1500)]))):.2f} dB")
print(f"Butter passband ripple: {np.ptp(20*np.log10(np.abs(H_b[(w_b>=500)&(w_b<=1500)]))):.2f} dB")

ex-sp-ch07-14

Hard

Implement the 2D FFT of a $256 \times 256$ grayscale image. Zero out all frequency components outside a circular region of radius $R$ in the frequency domain (ideal lowpass filter), then reconstruct with the inverse FFT. Show the effect for $R = 10, 30, 60, 120$ .

Show Hint

Create a circular mask using np.sqrt(kx**2 + ky**2) <= R.

Solution

Implementation

import numpy as np

# Create test image
x = np.linspace(-1, 1, 256)
X, Y = np.meshgrid(x, x)
img = np.sin(10*X) * np.cos(8*Y) + 0.5*np.exp(-(X**2+Y**2)/0.1)

F = np.fft.fft2(img)
F_shifted = np.fft.fftshift(F)

kx = np.arange(-128, 128)
KX, KY = np.meshgrid(kx, kx)
dist = np.sqrt(KX**2 + KY**2)

for R in [10, 30, 60, 120]:
    mask = (dist <= R).astype(float)
    F_filtered = np.fft.ifftshift(F_shifted * mask)
    img_filtered = np.fft.ifft2(F_filtered).real
    mse = np.mean((img - img_filtered)**2)
    print(f"R={R:3d}: MSE={mse:.6f}, kept={np.sum(mask)/mask.size*100:.1f}%")

ex-sp-ch07-15

Hard

Compare wavelet denoising with Wiener filtering for a piecewise-smooth signal in white Gaussian noise at SNR = 5 dB. Use the 'sym8' wavelet with soft thresholding and a Wiener filter estimated from the noisy signal's PSD. Report the MSE for both methods.

Show Hint

For Wiener filtering, estimate the signal PSD with Welch, then apply $H(f) = S_{xx}/(S_{xx} + S_{nn})$ .

Solution

Implementation

import pywt
import numpy as np
from scipy.signal import welch

np.random.seed(42)
N = 2048
t = np.linspace(0, 1, N)
x_clean = np.sign(np.sin(6*np.pi*t)) * (1 + 0.5*np.sin(20*np.pi*t))
snr_db = 5
sigma = np.std(x_clean) / 10**(snr_db/20)
x_noisy = x_clean + sigma * np.random.randn(N)

# Wavelet denoising
coeffs = pywt.wavedec(x_noisy, 'sym8', level=6)
sig_est = np.median(np.abs(coeffs[-1])) / 0.6745
thr = sig_est * np.sqrt(2*np.log(N))
coeffs_t = [coeffs[0]] + [pywt.threshold(c, thr, 'soft') for c in coeffs[1:]]
x_wav = pywt.waverec(coeffs_t, 'sym8')[:N]

# Wiener filter (frequency domain)
X_noisy = np.fft.fft(x_noisy)
f_w, Pxx = welch(x_noisy, nperseg=256)
Pnn = sigma**2  # known noise PSD (flat)
# Interpolate Welch PSD to FFT bins
from scipy.interpolate import interp1d
psd_interp = interp1d(f_w, Pxx, fill_value='extrapolate')
freqs = np.fft.fftfreq(N)
S_noisy = psd_interp(np.abs(freqs))
S_signal = np.maximum(S_noisy - Pnn, 0)
H_wiener = S_signal / (S_signal + Pnn)
x_wien = np.fft.ifft(X_noisy * H_wiener).real

mse_wav = np.mean((x_clean - x_wav)**2)
mse_wien = np.mean((x_clean - x_wien)**2)
print(f"Wavelet MSE: {mse_wav:.4f}, Wiener MSE: {mse_wien:.4f}")

ex-sp-ch07-16

Hard

Implement CWT-based time-frequency analysis of a signal with three components: a 50 Hz sinusoid (0-1 s), a 150 Hz sinusoid (0.5-1.5 s), and a transient Gaussian pulse at $t = 0.75$ s. Compare the CWT scalogram (Morlet wavelet) with the STFT spectrogram. Discuss which method better resolves the transient.

Show Hint

Use pywt.cwt with Morlet wavelet and scipy.signal.stft.

Solution

Implementation

import pywt
import numpy as np
from scipy.signal import stft

fs = 1000
t = np.arange(0, 2, 1/fs)
x = np.zeros(len(t))
x[t <= 1] += np.sin(2*np.pi*50*t[t <= 1])
x[(t >= 0.5) & (t <= 1.5)] += np.sin(2*np.pi*150*t[(t >= 0.5) & (t <= 1.5)])
x += 5 * np.exp(-((t - 0.75)*200)**2)

# CWT
scales = np.arange(1, 200)
coeffs, freqs = pywt.cwt(x, scales, 'morl', sampling_period=1/fs)

# STFT
f_stft, t_stft, Zxx = stft(x, fs=fs, nperseg=128, noverlap=120)

print(f"CWT shape: {coeffs.shape}, STFT shape: {Zxx.shape}")
print("CWT resolves the transient better due to adaptive resolution")

ex-sp-ch07-17

Challenge

Implement an OFDM transmitter and receiver in Python:

Generate random QPSK symbols for 64 subcarriers
Apply a 64-point IFFT to create the time-domain OFDM symbol
Add a cyclic prefix of length 16
Pass through a multipath channel $h = [1, 0, 0.5, 0, 0, 0.2]$
Add AWGN at SNR = 20 dB
Remove the cyclic prefix and apply the 64-point FFT
Equalize by dividing by the channel frequency response
Compute the bit error rate

Verify that the cyclic prefix eliminates inter-symbol interference.

Show Hint

QPSK symbols: $(\\pm 1 \\pm j)/\\sqrt{2}$ .

Channel frequency response: np.fft.fft(h, 64).

CP must be longer than the channel delay spread.

Solution

Implementation

import numpy as np

N = 64
cp_len = 16
n_symbols = 1000
snr_db = 20
h = np.array([1, 0, 0.5, 0, 0, 0.2])
H = np.fft.fft(h, N)

# Transmit
bits = np.random.randint(0, 4, (n_symbols, N))
qpsk_map = np.array([1+1j, 1-1j, -1+1j, -1-1j]) / np.sqrt(2)
X = qpsk_map[bits]

total_errors = 0
for i in range(n_symbols):
    # IFFT + CP
    x_time = np.fft.ifft(X[i])
    x_cp = np.concatenate([x_time[-cp_len:], x_time])

    # Channel + noise
    y = np.convolve(x_cp, h)[:len(x_cp)]
    noise_power = np.mean(np.abs(y)**2) / 10**(snr_db/10)
    y += np.sqrt(noise_power/2) * (np.random.randn(len(y))
                                    + 1j*np.random.randn(len(y)))

    # Remove CP + FFT + equalize
    y_no_cp = y[cp_len:]
    Y = np.fft.fft(y_no_cp)
    X_hat = Y / H

    # Detect
    det = np.argmin(np.abs(X_hat.reshape(-1,1) - qpsk_map.reshape(1,-1)), axis=1)
    total_errors += np.sum(det != bits[i])

ber = total_errors / (n_symbols * N)
print(f"BER: {ber:.6f}")

ex-sp-ch07-18

Challenge

Implement a complete spectral analysis pipeline for a non-stationary signal:

Generate a 10-second signal with three regimes:
- 0-3 s: two sinusoids at 100 Hz and 250 Hz
- 3-7 s: linear chirp from 100 Hz to 400 Hz
- 7-10 s: white noise filtered through a bandpass [150, 350] Hz
Compute the spectrogram, CWT scalogram, and Welch PSD for each regime
Implement an automatic regime detector that identifies transitions based on spectral feature changes
Compare the time-frequency resolution of the spectrogram vs CWT

Discuss which representation is best for each regime type.

Show Hint

Use spectral flatness (ratio of geometric to arithmetic mean) to distinguish tonal vs noise-like regimes.

Track spectral centroid over time to detect the chirp.

Solution

Implementation outline

import numpy as np
from scipy.signal import chirp, butter, sosfilt, spectrogram, welch
import pywt

fs = 2000
# Regime 1: tones
t1 = np.arange(0, 3, 1/fs)
x1 = np.sin(2*np.pi*100*t1) + np.sin(2*np.pi*250*t1)
# Regime 2: chirp
t2 = np.arange(0, 4, 1/fs)
x2 = chirp(t2, 100, 4, 400)
# Regime 3: filtered noise
t3 = np.arange(0, 3, 1/fs)
sos = butter(5, [150, 350], btype='band', fs=fs, output='sos')
x3 = sosfilt(sos, np.random.randn(len(t3)))
x = np.concatenate([x1, x2, x3])

# Spectrogram
f_s, t_s, Sxx = spectrogram(x, fs=fs, nperseg=256, noverlap=240)

# Spectral flatness per time frame
geometric_mean = np.exp(np.mean(np.log(Sxx + 1e-20), axis=0))
arithmetic_mean = np.mean(Sxx, axis=0)
flatness = geometric_mean / (arithmetic_mean + 1e-20)

# Regime detection
tonal_mask = flatness < 0.1
noise_mask = flatness > 0.3
print(f"Detected {np.sum(tonal_mask)} tonal frames, "
      f"{np.sum(noise_mask)} noise frames")

Chapter Summary References & Further Reading