1#!/usr/bin/env python3
2"""
3Fourier Series
4
5Demonstrates Fourier series analysis and synthesis:
6- Fourier coefficients for square wave, sawtooth, and triangular wave
7- Signal reconstruction using partial sums (Gibbs phenomenon)
8- Magnitude and phase spectra of the coefficients
9- Parseval's theorem: energy in time domain equals energy in frequency domain
10- Convergence behaviour with increasing number of harmonics
11
12The Fourier series of a periodic signal x(t) with period T is:
13 x(t) = sum_{n=-inf}^{inf} c_n * exp(j * 2*pi*n*t / T)
14
15The complex Fourier coefficients are:
16 c_n = (1/T) * integral_{-T/2}^{T/2} x(t) * exp(-j * 2*pi*n*t / T) dt
17
18Dependencies: numpy, matplotlib, scipy
19"""
20
21import numpy as np
22import matplotlib.pyplot as plt
23from scipy.integrate import quad
24
25
26# ---------------------------------------------------------------------------
27# Fourier coefficient computation
28# ---------------------------------------------------------------------------
29
30def fourier_coefficients_numerical(x_func, T: float, N_terms: int) -> np.ndarray:
31 """
32 Numerically compute complex Fourier coefficients c_n for n = -N..N.
33
34 c_n = (1/T) * integral_{0}^{T} x(t) * exp(-j * 2*pi*n*t / T) dt
35
36 Parameters
37 ----------
38 x_func : callable, periodic signal x(t)
39 T : period
40 N_terms : number of harmonics on each side (output has 2*N_terms+1 coeffs)
41
42 Returns
43 -------
44 c : complex array, shape (2*N_terms+1,), indices correspond to n = -N..N
45 """
46 c = np.zeros(2 * N_terms + 1, dtype=complex)
47 for idx, n in enumerate(range(-N_terms, N_terms + 1)):
48 def integrand_real(t, n=n):
49 return x_func(t) * np.cos(2 * np.pi * n * t / T)
50
51 def integrand_imag(t, n=n):
52 return -x_func(t) * np.sin(2 * np.pi * n * t / T)
53
54 real_part, _ = quad(integrand_real, 0, T)
55 imag_part, _ = quad(integrand_imag, 0, T)
56 c[idx] = (real_part + 1j * imag_part) / T
57 return c
58
59
60def reconstruct_signal(c: np.ndarray, T: float, t: np.ndarray) -> np.ndarray:
61 """
62 Reconstruct periodic signal from complex Fourier coefficients.
63
64 x(t) = sum_{n=-N}^{N} c_n * exp(j * 2*pi*n*t / T)
65
66 Parameters
67 ----------
68 c : complex Fourier coefficients, length 2*N+1
69 T : period
70 t : time array
71
72 Returns
73 -------
74 x_reconstructed : real-valued signal (imaginary part should be ~0)
75 """
76 N = (len(c) - 1) // 2
77 x_rec = np.zeros(len(t), dtype=complex)
78 for idx, n in enumerate(range(-N, N + 1)):
79 x_rec += c[idx] * np.exp(1j * 2 * np.pi * n * t / T)
80 return x_rec.real
81
82
83# ---------------------------------------------------------------------------
84# Signal definitions
85# ---------------------------------------------------------------------------
86
87def square_wave(t: np.ndarray, T: float = 1.0, duty: float = 0.5) -> np.ndarray:
88 """Square wave: +1 for first half-period, -1 for second half-period."""
89 phase = (t % T) / T # normalized phase in [0, 1)
90 return np.where(phase < duty, 1.0, -1.0)
91
92
93def sawtooth_wave(t: np.ndarray, T: float = 1.0) -> np.ndarray:
94 """Sawtooth wave: linear ramp from -1 to +1 over one period."""
95 return 2 * ((t / T) % 1) - 1
96
97
98def triangular_wave(t: np.ndarray, T: float = 1.0) -> np.ndarray:
99 """Triangular wave: linearly rises then falls symmetrically."""
100 phase = (t / T) % 1
101 return np.where(phase < 0.5, 4 * phase - 1, -4 * phase + 3)
102
103
104# ---------------------------------------------------------------------------
105# Demonstrations
106# ---------------------------------------------------------------------------
107
108def demo_gibbs_phenomenon():
109 """
110 Show Gibbs phenomenon: partial Fourier sums of a square wave overshoot
111 near discontinuities by ~9% regardless of the number of terms.
112
113 Analytical Fourier series of the square wave (odd harmonics only):
114 x(t) = (4/pi) * sum_{k=0}^{inf} sin((2k+1)*2*pi*t/T) / (2k+1)
115 """
116 print("=" * 60)
117 print("GIBBS PHENOMENON (SQUARE WAVE RECONSTRUCTION)")
118 print("=" * 60)
119
120 T = 1.0
121 t = np.linspace(0, 2 * T, 2000)
122 x_true = square_wave(t, T)
123
124 N_list = [1, 3, 7, 19, 49]
125
126 fig, axes = plt.subplots(len(N_list), 1, figsize=(10, 12), sharex=True)
127 fig.suptitle("Gibbs Phenomenon: Square Wave Partial Sums", fontsize=13, fontweight='bold')
128
129 for ax, N in zip(axes, N_list):
130 # Reconstruct using only odd harmonics up to 2N-1
131 x_rec = np.zeros_like(t)
132 for k in range(N):
133 n = 2 * k + 1 # odd harmonic index
134 x_rec += (4 / np.pi) * np.sin(n * 2 * np.pi * t / T) / n
135
136 overshoot = (x_rec.max() - 1.0) * 100
137 ax.plot(t, x_true, 'lightgray', linewidth=2, label='True')
138 ax.plot(t, x_rec, 'b', linewidth=1, label=f'N={N} terms')
139 ax.set_ylabel("Amplitude")
140 ax.set_ylim(-1.5, 1.5)
141 ax.legend(loc='upper right', fontsize=8)
142 ax.grid(True, alpha=0.3)
143 ax.set_title(f"N = {N} harmonics, overshoot ā {overshoot:.1f}%", fontsize=9)
144 print(f" N={N:2d}: overshoot = {overshoot:.2f}%")
145
146 axes[-1].set_xlabel("t (s)")
147 plt.tight_layout()
148 plt.show()
149
150 print("\nGibbs phenomenon: overshoot converges to ~9% (Ļ/4 - 1 ā 8.9%)")
151 print("Overshoot does NOT decrease with more terms (only gets narrower)")
152
153
154def demo_spectra():
155 """
156 Plot magnitude and phase spectra for square, sawtooth, and triangular waves.
157
158 The frequency spectrum shows which harmonics are present and their amplitudes.
159 """
160 print("\n" + "=" * 60)
161 print("MAGNITUDE AND PHASE SPECTRA")
162 print("=" * 60)
163
164 T = 1.0
165 N_terms = 15 # compute up to 15th harmonic
166 f0 = 1.0 / T # fundamental frequency
167
168 signals = {
169 'Square Wave': lambda t: np.where((t % T) / T < 0.5, 1.0, -1.0),
170 'Sawtooth Wave': lambda t: 2 * ((t / T) % 1) - 1,
171 'Triangular Wave': lambda t: triangular_wave(np.atleast_1d(t), T).item()
172 if np.isscalar(t) else triangular_wave(t, T),
173 }
174
175 fig, axes = plt.subplots(3, 2, figsize=(12, 10))
176 fig.suptitle("Fourier Spectra of Periodic Signals", fontsize=13, fontweight='bold')
177
178 for row, (name, func) in enumerate(signals.items()):
179 c = fourier_coefficients_numerical(func, T, N_terms)
180 n_vals = np.arange(-N_terms, N_terms + 1)
181 freqs = n_vals * f0
182
183 magnitude = np.abs(c)
184 phase = np.angle(c)
185 # Zero out negligible phases (below numerical noise floor)
186 phase[magnitude < 1e-10] = 0.0
187
188 # One-sided for clarity: show n >= 0 only
189 pos_mask = n_vals >= 0
190 n_pos = n_vals[pos_mask]
191 mag_pos = 2 * magnitude[pos_mask] # double-sided to one-sided
192 mag_pos[n_pos == 0] /= 2 # DC term is not doubled
193 pha_pos = phase[pos_mask]
194
195 axes[row, 0].stem(n_pos, mag_pos, basefmt='k-',
196 linefmt='C0-', markerfmt='C0o')
197 axes[row, 0].set_title(f"{name}: Magnitude Spectrum")
198 axes[row, 0].set_xlabel("Harmonic number n")
199 axes[row, 0].set_ylabel("|c_n|")
200 axes[row, 0].grid(True, alpha=0.3)
201
202 axes[row, 1].stem(n_pos, pha_pos, basefmt='k-',
203 linefmt='C1-', markerfmt='C1o')
204 axes[row, 1].set_title(f"{name}: Phase Spectrum")
205 axes[row, 1].set_xlabel("Harmonic number n")
206 axes[row, 1].set_ylabel("Phase (rad)")
207 axes[row, 1].set_ylim(-np.pi - 0.3, np.pi + 0.3)
208 axes[row, 1].axhline(0, color='k', linewidth=0.5)
209 axes[row, 1].axhline( np.pi/2, color='gray', linestyle='--', alpha=0.4)
210 axes[row, 1].axhline(-np.pi/2, color='gray', linestyle='--', alpha=0.4)
211 axes[row, 1].grid(True, alpha=0.3)
212
213 # Print dominant harmonics
214 print(f"{name}:")
215 for idx_n, (n, m) in enumerate(zip(n_pos[:8], mag_pos[:8])):
216 if m > 0.01:
217 print(f" n={n:2d}: |c_n|={m:.4f}, phase={pha_pos[idx_n]:.3f} rad")
218
219 plt.tight_layout()
220 plt.show()
221
222
223def demo_parseval():
224 """
225 Verify Parseval's theorem numerically.
226
227 Parseval's theorem states that the total signal power equals the sum
228 of squared magnitudes of Fourier coefficients:
229
230 (1/T) * integral_T |x(t)|^2 dt = sum_{n=-inf}^{inf} |c_n|^2
231
232 This is a statement of energy conservation between time and frequency domains.
233 """
234 print("\n" + "=" * 60)
235 print("PARSEVAL'S THEOREM VERIFICATION")
236 print("=" * 60)
237
238 T = 1.0
239 dt = 0.0001
240 t = np.arange(0, T, dt)
241 N_terms = 50 # use many terms for good approximation
242
243 signals = {
244 'Square Wave': square_wave(t, T),
245 'Sawtooth Wave': sawtooth_wave(t, T),
246 'Triangular Wave': triangular_wave(t, T),
247 }
248 signal_funcs = {
249 'Square Wave': lambda t: np.where((t % T) / T < 0.5, 1.0, -1.0),
250 'Sawtooth Wave': lambda t: 2 * ((t / T) % 1) - 1,
251 'Triangular Wave': lambda t: triangular_wave(np.atleast_1d(t), T).item()
252 if np.isscalar(t) else triangular_wave(t, T),
253 }
254
255 for name in signals:
256 x = signals[name]
257 func = signal_funcs[name]
258
259 # Time-domain power: (1/T) * integral |x(t)|^2 dt
260 power_time = np.mean(x ** 2)
261
262 # Frequency-domain power: sum |c_n|^2
263 c = fourier_coefficients_numerical(func, T, N_terms)
264 power_freq = np.sum(np.abs(c) ** 2)
265
266 print(f"{name}:")
267 print(f" Time-domain power: {power_time:.6f}")
268 print(f" Frequency-domain power: {power_freq:.6f}")
269 print(f" Relative error: {abs(power_time - power_freq) / power_time * 100:.3f}%")
270
271
272def demo_convergence():
273 """
274 Show convergence rate of Fourier series for different waveforms.
275
276 Convergence rate depends on signal smoothness:
277 - Discontinuities (square, sawtooth): coefficients decay as 1/n -> slow
278 - Continuous but non-smooth (triangular): coefficients decay as 1/n^2 -> faster
279 - Smooth signals: coefficients decay exponentially -> fastest
280 """
281 print("\n" + "=" * 60)
282 print("FOURIER SERIES CONVERGENCE COMPARISON")
283 print("=" * 60)
284
285 T = 1.0
286 dt = 0.001
287 t = np.linspace(0, T, int(T / dt), endpoint=False)
288
289 x_square = square_wave(t, T)
290 x_saw = sawtooth_wave(t, T)
291 x_tri = triangular_wave(t, T)
292
293 signal_funcs = {
294 'Square (1/n decay)': lambda t: np.where((t % T) / T < 0.5, 1.0, -1.0),
295 'Sawtooth (1/n decay)': lambda t: 2 * ((t / T) % 1) - 1,
296 'Triangular (1/nĀ² decay)': lambda t: triangular_wave(np.atleast_1d(t), T).item()
297 if np.isscalar(t) else triangular_wave(t, T),
298 }
299
300 N_max = 30
301 N_list = np.arange(1, N_max + 1)
302
303 fig, axes = plt.subplots(1, 2, figsize=(12, 5))
304 fig.suptitle("Fourier Series Convergence", fontsize=13, fontweight='bold')
305
306 true_signals = [x_square, x_saw, x_tri]
307
308 for (name, func), x_true in zip(signal_funcs.items(), true_signals):
309 errors = []
310 c_all = fourier_coefficients_numerical(func, T, N_max)
311
312 for N in N_list:
313 # Use coefficients up to harmonic N
314 c_trunc = np.zeros_like(c_all)
315 center = N_max
316 c_trunc[center - N: center + N + 1] = c_all[center - N: center + N + 1]
317 x_rec = reconstruct_signal(c_trunc, T, t)
318 rmse = np.sqrt(np.mean((x_rec - x_true) ** 2))
319 errors.append(rmse)
320
321 axes[0].plot(N_list, errors, marker='o', markersize=3, label=name)
322 axes[1].semilogy(N_list, errors, marker='o', markersize=3, label=name)
323
324 print(f"{name}: RMSE at N=1: {errors[0]:.4f}, N=10: {errors[9]:.4f}, N=30: {errors[-1]:.4f}")
325
326 axes[0].set_title("Convergence (linear scale)")
327 axes[0].set_xlabel("Number of harmonics N")
328 axes[0].set_ylabel("RMSE")
329 axes[0].legend()
330 axes[0].grid(True, alpha=0.3)
331
332 axes[1].set_title("Convergence (log scale)")
333 axes[1].set_xlabel("Number of harmonics N")
334 axes[1].set_ylabel("RMSE (log)")
335 axes[1].legend()
336 axes[1].grid(True, alpha=0.3)
337
338 plt.tight_layout()
339 plt.show()
340
341 print("\nSmoother signals require fewer harmonics for the same accuracy.")
342 print("Triangular wave converges faster (1/nĀ²) than square/sawtooth (1/n).")
343
344
345if __name__ == "__main__":
346 demo_gibbs_phenomenon()
347 demo_spectra()
348 demo_parseval()
349 demo_convergence()
350
351 print("\n" + "=" * 60)
352 print("All demonstrations completed!")
353 print("=" * 60)