1#!/usr/bin/env python3
2"""
3Sampling, Aliasing, and Signal Reconstruction
4==============================================
5
6The Nyquist-Shannon Sampling Theorem states:
7
8 A bandlimited signal with maximum frequency f_max can be perfectly
9 reconstructed from samples taken at a rate fs > 2 * f_max.
10
11 The minimum rate fs = 2 * f_max is called the Nyquist rate.
12
13When fs < 2 * f_max (under-sampling), spectral copies overlap and the
14signal cannot be recovered — this is called ALIASING.
15
16Topics Covered:
17 1. Sampling a continuous-time signal at multiple rates
18 2. Aliasing demonstration (fs < 2*fmax violates Nyquist)
19 3. Proper sampling without aliasing (fs > 2*fmax)
20 4. Sinc (ideal) interpolation for perfect reconstruction
21 5. Spectrum before and after sampling (showing spectral copies)
22 6. Anti-aliasing (low-pass) filter applied before sampling
23
24Author: Educational example for Signal Processing
25License: MIT
26"""
27
28import numpy as np
29import matplotlib.pyplot as plt
30from scipy.signal import butter, filtfilt, resample
31
32
33# =============================================================================
34# 1. Signal and spectrum utilities
35# =============================================================================
36
37def make_continuous_signal(t):
38 """
39 Create a test bandlimited signal composed of two sinusoids.
40
41 x(t) = sin(2π·3·t) + 0.5·sin(2π·7·t)
42
43 Maximum frequency: f_max = 7 Hz
44 Nyquist rate: fs_nyquist = 14 Hz
45 """
46 f1, f2 = 3.0, 7.0 # Hz
47 return np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
48
49
50def compute_spectrum(x, fs):
51 """
52 Compute one-sided magnitude spectrum using FFT.
53
54 Args:
55 x (ndarray): Signal samples
56 fs (float): Sampling frequency (Hz)
57
58 Returns:
59 tuple: (f, mag) — positive frequency axis and magnitude
60 """
61 N = len(x)
62 X = np.fft.rfft(x)
63 f = np.fft.rfftfreq(N, d=1.0 / fs)
64 mag = (2.0 / N) * np.abs(X) # two-sided → one-sided amplitude
65 mag[0] /= 2 # DC component is not doubled
66 return f, mag
67
68
69def sinc_interpolation(x_samples, t_samples, t_recon):
70 """
71 Reconstruct a continuous signal from samples using ideal sinc interpolation.
72
73 Whittaker-Shannon formula:
74 x_r(t) = sum_n x[n] * sinc( fs*(t - n/fs) )
75 = sum_n x[n] * sinc( (t - t_n) * fs )
76
77 This is exact for bandlimited signals sampled above the Nyquist rate.
78
79 Args:
80 x_samples (ndarray): Sampled signal values
81 t_samples (ndarray): Sample time instants
82 t_recon (ndarray): Time points for reconstruction
83
84 Returns:
85 ndarray: Reconstructed signal at t_recon
86 """
87 fs = 1.0 / (t_samples[1] - t_samples[0])
88 # Broadcast: (len(t_recon), len(t_samples))
89 T = t_recon[:, np.newaxis] - t_samples[np.newaxis, :]
90 # np.sinc(x) = sin(pi*x) / (pi*x)
91 sinc_matrix = np.sinc(T * fs)
92 return sinc_matrix @ x_samples
93
94
95def lowpass_filter(x, fs, cutoff_hz, order=6):
96 """
97 Apply a Butterworth low-pass filter (anti-aliasing filter).
98
99 Args:
100 x (ndarray): Input signal
101 fs (float): Sampling frequency (Hz)
102 cutoff_hz (float): Cutoff frequency (Hz)
103 order (int): Filter order
104
105 Returns:
106 ndarray: Filtered signal
107 """
108 nyq = fs / 2.0
109 normal_cutoff = cutoff_hz / nyq
110 b, a = butter(order, normal_cutoff, btype='low', analog=False)
111 return filtfilt(b, a, x)
112
113
114# =============================================================================
115# 2. Demonstration functions
116# =============================================================================
117
118def demo_aliasing():
119 """
120 Show what happens when we sample below the Nyquist rate.
121
122 The signal has f_max = 7 Hz, so Nyquist rate = 14 Hz.
123 We sample at fs = 10 Hz (< 14 Hz) and observe aliasing:
124 the 7 Hz component aliases to 10-7 = 3 Hz, adding to the real 3 Hz tone.
125 """
126 print("=" * 60)
127 print("Demo 1: Aliasing (Under-Sampling)")
128 print("=" * 60)
129
130 # "Continuous-time" reference: very high sampling rate
131 fs_ref = 1000.0
132 t_ref = np.arange(0, 2.0, 1.0 / fs_ref)
133 x_ref = make_continuous_signal(t_ref)
134
135 # Under-sampling: fs < 2 * f_max
136 fs_low = 10.0 # Hz (Nyquist = 14 Hz → aliased!)
137 t_low = np.arange(0, 2.0, 1.0 / fs_low)
138 x_low = make_continuous_signal(t_low)
139
140 # Nyquist-rate sampling
141 fs_nyq = 14.0
142 t_nyq = np.arange(0, 2.0, 1.0 / fs_nyq)
143 x_nyq = make_continuous_signal(t_nyq)
144
145 # Over-sampling (safe)
146 fs_high = 40.0
147 t_high = np.arange(0, 2.0, 1.0 / fs_high)
148 x_high = make_continuous_signal(t_high)
149
150 # Spectra
151 f_ref, S_ref = compute_spectrum(x_ref, fs_ref)
152 f_low, S_low = compute_spectrum(x_low, fs_low)
153 f_nyq, S_nyq = compute_spectrum(x_nyq, fs_nyq)
154 f_high, S_high = compute_spectrum(x_high, fs_high)
155
156 # Aliased frequency for 7 Hz component at fs=10:
157 # alias = fs - f = 10 - 7 = 3 Hz (coincides with 3 Hz component!)
158 print(f" f_max = 7 Hz, Nyquist rate = 14 Hz")
159 print(f" Sampling at fs={fs_low} Hz: 7 Hz aliases to {fs_low-7.0:.0f} Hz")
160 print(f" Sampling at fs={fs_nyq} Hz: exactly at Nyquist (edge case)")
161 print(f" Sampling at fs={fs_high} Hz: safely above Nyquist — no aliasing")
162
163 fig, axes = plt.subplots(2, 2, figsize=(13, 8))
164 fig.suptitle("Aliasing: Effect of Sampling Rate", fontsize=14, fontweight='bold')
165
166 configs = [
167 (f_ref, S_ref, fs_ref, "Reference (fs=1000 Hz)", "green"),
168 (f_low, S_low, fs_low, f"Under-sampled (fs={fs_low} Hz) — ALIASED", "red"),
169 (f_nyq, S_nyq, fs_nyq, f"Nyquist rate (fs={fs_nyq} Hz)", "orange"),
170 (f_high, S_high, fs_high, f"Over-sampled (fs={fs_high} Hz) — safe", "blue"),
171 ]
172
173 for ax, (fi, Si, fsi, title, color) in zip(axes.flat, configs):
174 ax.stem(fi, Si, linefmt=color, markerfmt=color+'o', basefmt='k-')
175 ax.axvline(7.0, color='gray', linestyle='--', alpha=0.5, label='f_max=7 Hz')
176 ax.axvline(3.0, color='purple', linestyle=':', alpha=0.5, label='f=3 Hz')
177 if fsi <= 14.0:
178 ax.axvline(fsi / 2, color='red', linestyle='--', alpha=0.4,
179 label=f'fs/2={fsi/2:.0f} Hz')
180 ax.set_xlim(0, min(fsi / 2 + 1, 25))
181 ax.set_title(title, fontsize=10)
182 ax.set_xlabel("Frequency (Hz)")
183 ax.set_ylabel("Amplitude")
184 ax.legend(fontsize=7)
185 ax.grid(True, alpha=0.3)
186
187 plt.tight_layout()
188 plt.savefig('/tmp/05_aliasing_spectra.png', dpi=150)
189 print(" Saved plot to /tmp/05_aliasing_spectra.png")
190 plt.show()
191
192
193def demo_reconstruction():
194 """
195 Demonstrate perfect reconstruction via sinc interpolation when fs > 2*fmax.
196
197 Compare:
198 - Under-sampled (aliased): reconstruction fails
199 - Properly sampled: reconstruction matches original perfectly
200 """
201 print("\n" + "=" * 60)
202 print("Demo 2: Signal Reconstruction via Sinc Interpolation")
203 print("=" * 60)
204
205 # Original "continuous" signal
206 t_ref = np.arange(0, 1.0, 0.001)
207 x_ref = make_continuous_signal(t_ref)
208
209 # Reconstruction time axis
210 t_recon = np.arange(0.05, 0.95, 0.002)
211
212 # --- Under-sampling (aliased) ---
213 fs_bad = 10.0
214 t_bad = np.arange(0, 1.0, 1.0 / fs_bad)
215 x_bad = make_continuous_signal(t_bad)
216 x_recon_bad = sinc_interpolation(x_bad, t_bad, t_recon)
217
218 # --- Proper sampling ---
219 fs_good = 30.0
220 t_good = np.arange(0, 1.0, 1.0 / fs_good)
221 x_good = make_continuous_signal(t_good)
222 x_recon_good = sinc_interpolation(x_good, t_good, t_recon)
223
224 # Reconstruction error
225 x_ref_recon_pts = make_continuous_signal(t_recon)
226 err_bad = np.sqrt(np.mean((x_recon_bad - x_ref_recon_pts) ** 2))
227 err_good = np.sqrt(np.mean((x_recon_good - x_ref_recon_pts) ** 2))
228
229 print(f" Under-sampled (fs={fs_bad} Hz): RMS reconstruction error = {err_bad:.4f}")
230 print(f" Properly sampled (fs={fs_good} Hz): RMS reconstruction error = {err_good:.6f}")
231
232 fig, axes = plt.subplots(1, 2, figsize=(13, 5))
233 fig.suptitle("Sinc Interpolation: Reconstruction Quality", fontsize=14, fontweight='bold')
234
235 for ax, (fs, t_s, x_s, x_r, err, title_extra, color) in zip(axes, [
236 (fs_bad, t_bad, x_bad, x_recon_bad, err_bad, "ALIASED", "red"),
237 (fs_good, t_good, x_good, x_recon_good, err_good, "Correct", "blue"),
238 ]):
239 ax.plot(t_ref, x_ref, 'k-', lw=1.5, label='Original', alpha=0.6)
240 ax.stem(t_s, x_s, linefmt=color, markerfmt=color+'o',
241 basefmt='none', label=f'Samples (fs={fs} Hz)')
242 ax.plot(t_recon, x_r, color, lw=2, linestyle='--', label='Reconstructed')
243 ax.set_xlim(0, 1)
244 ax.set_ylim(-2, 2)
245 ax.set_title(f"fs={fs} Hz — {title_extra}\nRMS error = {err:.4f}")
246 ax.set_xlabel("Time (s)")
247 ax.set_ylabel("Amplitude")
248 ax.legend(fontsize=8)
249 ax.grid(True, alpha=0.3)
250
251 plt.tight_layout()
252 plt.savefig('/tmp/05_reconstruction.png', dpi=150)
253 print(" Saved plot to /tmp/05_reconstruction.png")
254 plt.show()
255
256
257def demo_antialias_filter():
258 """
259 Show how a low-pass anti-aliasing filter prevents aliasing before sampling.
260
261 Workflow:
262 1. Original broadband signal (contains high-frequency components)
263 2. Apply anti-aliasing LPF with cutoff = fs/2
264 3. Sample the filtered signal — no aliasing occurs
265 """
266 print("\n" + "=" * 60)
267 print("Demo 3: Anti-Aliasing Filter")
268 print("=" * 60)
269
270 # Broadband signal: 3 Hz + 7 Hz + 15 Hz + 22 Hz
271 fs_ref = 2000.0
272 t_ref = np.arange(0, 1.0, 1.0 / fs_ref)
273 x_broad = (np.sin(2 * np.pi * 3 * t_ref) +
274 np.sin(2 * np.pi * 7 * t_ref) +
275 np.sin(2 * np.pi * 15 * t_ref) +
276 0.5 * np.sin(2 * np.pi * 22 * t_ref))
277
278 # Target sampling rate
279 fs_target = 20.0 # Hz — Nyquist = 10 Hz
280
281 # Without anti-aliasing: sample directly
282 t_sampled = np.arange(0, 1.0, 1.0 / fs_target)
283 x_no_aa = np.interp(t_sampled, t_ref, x_broad)
284
285 # With anti-aliasing: apply LPF first
286 cutoff = fs_target / 2.0 * 0.9 # slightly below Nyquist
287 x_filtered = lowpass_filter(x_broad, fs_ref, cutoff_hz=cutoff)
288 x_with_aa = np.interp(t_sampled, t_ref, x_filtered)
289
290 # Spectra
291 f_broad, S_broad = compute_spectrum(x_broad, fs_ref)
292 f_filtered, S_filtered = compute_spectrum(x_filtered, fs_ref)
293 f_no_aa, S_no_aa = compute_spectrum(x_no_aa, fs_target)
294 f_with_aa, S_with_aa = compute_spectrum(x_with_aa, fs_target)
295
296 print(f" Broadband signal: 3, 7, 15, 22 Hz components")
297 print(f" Target fs = {fs_target} Hz → LPF cutoff = {cutoff:.1f} Hz")
298 print(f" Without AA filter: components at 15 Hz alias to {fs_target-15:.0f} Hz")
299 print(f" components at 22 Hz alias to {fs_target-22+fs_target:.0f} Hz -> {abs(22-fs_target):.0f} Hz")
300
301 fig, axes = plt.subplots(2, 2, figsize=(13, 8))
302 fig.suptitle("Anti-Aliasing Filter Effect", fontsize=14, fontweight='bold')
303
304 # Original broadband spectrum
305 mask_broad = f_broad <= 30
306 axes[0, 0].plot(f_broad[mask_broad], S_broad[mask_broad], 'b')
307 axes[0, 0].axvline(fs_target / 2, color='red', linestyle='--', label=f'fs/2={fs_target/2} Hz')
308 axes[0, 0].set_title("Original Broadband Spectrum")
309 axes[0, 0].set_xlabel("Frequency (Hz)")
310 axes[0, 0].set_ylabel("Amplitude")
311 axes[0, 0].legend()
312 axes[0, 0].grid(True, alpha=0.3)
313
314 # After LPF
315 axes[0, 1].plot(f_filtered[mask_broad], S_filtered[mask_broad], 'g')
316 axes[0, 1].axvline(fs_target / 2, color='red', linestyle='--', label=f'LPF cutoff≈{cutoff:.0f} Hz')
317 axes[0, 1].set_title("After Anti-Aliasing LPF")
318 axes[0, 1].set_xlabel("Frequency (Hz)")
319 axes[0, 1].set_ylabel("Amplitude")
320 axes[0, 1].legend()
321 axes[0, 1].grid(True, alpha=0.3)
322
323 # Sampled without AA — aliased
324 axes[1, 0].stem(f_no_aa, S_no_aa, linefmt='r', markerfmt='ro', basefmt='k-')
325 axes[1, 0].set_title(f"Sampled at fs={fs_target} Hz (NO anti-aliasing) — ALIASED")
326 axes[1, 0].set_xlabel("Frequency (Hz)")
327 axes[1, 0].set_ylabel("Amplitude")
328 axes[1, 0].set_xlim(0, fs_target / 2 + 1)
329 axes[1, 0].grid(True, alpha=0.3)
330
331 # Sampled with AA — clean
332 axes[1, 1].stem(f_with_aa, S_with_aa, linefmt='g', markerfmt='go', basefmt='k-')
333 axes[1, 1].set_title(f"Sampled at fs={fs_target} Hz (WITH anti-aliasing) — Clean")
334 axes[1, 1].set_xlabel("Frequency (Hz)")
335 axes[1, 1].set_ylabel("Amplitude")
336 axes[1, 1].set_xlim(0, fs_target / 2 + 1)
337 axes[1, 1].grid(True, alpha=0.3)
338
339 plt.tight_layout()
340 plt.savefig('/tmp/05_antialias.png', dpi=150)
341 print(" Saved plot to /tmp/05_antialias.png")
342 plt.show()
343
344
345if __name__ == "__main__":
346 print("Sampling, Aliasing, and Signal Reconstruction")
347 print("=" * 60)
348 print("Test signal: x(t) = sin(2π·3·t) + 0.5·sin(2π·7·t)")
349 print(" f_max = 7 Hz → Nyquist rate = 14 Hz\n")
350
351 demo_aliasing()
352 demo_reconstruction()
353 demo_antialias_filter()
354
355 print("\n" + "=" * 60)
356 print("Key Takeaways:")
357 print(" - Nyquist theorem: fs > 2*f_max required for perfect reconstruction")
358 print(" - Aliasing folds high frequencies back: f_alias = |f - k*fs|")
359 print(" - Sinc interpolation gives exact reconstruction (bandlimited signals)")
360 print(" - Anti-aliasing LPF removes components above fs/2 before sampling")
361 print(" - Practical systems oversample (typically fs > 5*f_max) for safety")
362 print("=" * 60)