1"""
2Fourier Analysis - Series, Transforms, and Spectral Analysis
3
4This script demonstrates:
5- Fourier series coefficients computation
6- Fast Fourier Transform (FFT)
7- Spectral analysis of signals
8- Filtering in frequency domain
9- Parseval's theorem verification
10- Windowing techniques
11"""
12
13import numpy as np
14
15try:
16 import matplotlib.pyplot as plt
17 HAS_MATPLOTLIB = True
18except ImportError:
19 HAS_MATPLOTLIB = False
20 print("Warning: matplotlib not available, skipping visualizations")
21
22
23def fourier_series_coefficients(f, T, n_harmonics):
24 """
25 Compute Fourier series coefficients for periodic function f with period T
26 f(t) = a0/2 + sum(an*cos(nωt) + bn*sin(nωt))
27 """
28 omega = 2 * np.pi / T
29 N = 1000 # Number of integration points
30
31 t = np.linspace(0, T, N)
32 dt = T / N
33
34 # a0 (DC component)
35 a0 = (2 / T) * np.sum(f(t)) * dt
36
37 # an and bn coefficients
38 an = np.zeros(n_harmonics)
39 bn = np.zeros(n_harmonics)
40
41 for n in range(1, n_harmonics + 1):
42 an[n-1] = (2 / T) * np.sum(f(t) * np.cos(n * omega * t)) * dt
43 bn[n-1] = (2 / T) * np.sum(f(t) * np.sin(n * omega * t)) * dt
44
45 return a0, an, bn
46
47
48def reconstruct_from_fourier(t, T, a0, an, bn):
49 """Reconstruct signal from Fourier coefficients"""
50 omega = 2 * np.pi / T
51 signal = a0 / 2 * np.ones_like(t)
52
53 for n in range(len(an)):
54 signal += an[n] * np.cos((n + 1) * omega * t)
55 signal += bn[n] * np.sin((n + 1) * omega * t)
56
57 return signal
58
59
60def parseval_theorem_check(f, T, a0, an, bn):
61 """
62 Verify Parseval's theorem:
63 (1/T) ∫|f(t)|² dt = (a0²/4) + (1/2)∑(an² + bn²)
64 """
65 # Left side: time domain energy
66 N = 1000
67 t = np.linspace(0, T, N)
68 dt = T / N
69 time_energy = (1 / T) * np.sum(np.abs(f(t))**2) * dt
70
71 # Right side: frequency domain energy
72 freq_energy = (a0**2) / 4
73 for n in range(len(an)):
74 freq_energy += 0.5 * (an[n]**2 + bn[n]**2)
75
76 return time_energy, freq_energy
77
78
79def apply_window(signal, window_type='hanning'):
80 """Apply window function to signal"""
81 N = len(signal)
82
83 if window_type == 'hanning':
84 window = 0.5 * (1 - np.cos(2 * np.pi * np.arange(N) / N))
85 elif window_type == 'hamming':
86 window = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(N) / N)
87 elif window_type == 'blackman':
88 n = np.arange(N)
89 window = 0.42 - 0.5 * np.cos(2*np.pi*n/N) + 0.08 * np.cos(4*np.pi*n/N)
90 elif window_type == 'rectangular':
91 window = np.ones(N)
92 else:
93 raise ValueError(f"Unknown window type: {window_type}")
94
95 return signal * window, window
96
97
98def lowpass_filter(signal, cutoff_freq, sampling_freq):
99 """Apply ideal lowpass filter in frequency domain"""
100 N = len(signal)
101 fft_signal = np.fft.fft(signal)
102 freqs = np.fft.fftfreq(N, 1/sampling_freq)
103
104 # Create filter
105 filter_mask = np.abs(freqs) <= cutoff_freq
106 fft_filtered = fft_signal * filter_mask
107
108 # Inverse FFT
109 filtered_signal = np.fft.ifft(fft_filtered).real
110
111 return filtered_signal, freqs, filter_mask
112
113
114# ============================================================================
115# MAIN DEMONSTRATIONS
116# ============================================================================
117
118print("=" * 70)
119print("FOURIER ANALYSIS - SERIES, TRANSFORMS, AND SPECTRAL ANALYSIS")
120print("=" * 70)
121
122# Test 1: Fourier series for square wave
123print("\n1. FOURIER SERIES - SQUARE WAVE")
124print("-" * 70)
125T = 2 * np.pi
126square_wave = lambda t: np.where(np.mod(t, T) < T/2, 1, -1)
127
128a0, an, bn = fourier_series_coefficients(square_wave, T, n_harmonics=10)
129print(f"Period T = {T:.4f}")
130print(f"DC component a0 = {a0:.6f}")
131print(f"\nFirst 5 Fourier coefficients:")
132for n in range(5):
133 print(f" a{n+1} = {an[n]:8.6f}, b{n+1} = {bn[n]:8.6f}")
134
135# Analytical result for square wave: bn = 4/(nπ) for odd n
136print(f"\nAnalytical (odd n): b_n = 4/(nπ)")
137for n in [1, 3, 5]:
138 analytical = 4 / (n * np.pi)
139 print(f" b{n} analytical = {analytical:.6f}")
140
141# Test 2: Fourier series for sawtooth wave
142print("\n2. FOURIER SERIES - SAWTOOTH WAVE")
143print("-" * 70)
144T = 2 * np.pi
145sawtooth = lambda t: (np.mod(t, T) / T) * 2 - 1
146
147a0, an, bn = fourier_series_coefficients(sawtooth, T, n_harmonics=10)
148print(f"DC component a0 = {a0:.6f}")
149print(f"\nFirst 5 Fourier coefficients:")
150for n in range(5):
151 print(f" a{n+1} = {an[n]:8.6f}, b{n+1} = {bn[n]:8.6f}")
152
153# Analytical: bn = -2/(nπ) for all n
154print(f"\nAnalytical: b_n = -2/(nπ)")
155for n in [1, 2, 3]:
156 analytical = -2 / (n * np.pi)
157 print(f" b{n} analytical = {analytical:.6f}")
158
159# Test 3: FFT of composite signal
160print("\n3. FAST FOURIER TRANSFORM (FFT)")
161print("-" * 70)
162# Create composite signal: 5Hz + 15Hz + 25Hz
163fs = 1000 # Sampling frequency
164t = np.linspace(0, 1, fs, endpoint=False)
165signal = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*15*t) + 0.3*np.sin(2*np.pi*25*t)
166
167# Add noise
168np.random.seed(42)
169signal_noisy = signal + 0.1 * np.random.randn(len(signal))
170
171# Compute FFT
172fft_result = np.fft.fft(signal_noisy)
173freqs = np.fft.fftfreq(len(signal_noisy), 1/fs)
174magnitude = np.abs(fft_result) / len(signal_noisy)
175
176# Find peaks
177positive_freqs = freqs[:len(freqs)//2]
178positive_magnitude = magnitude[:len(magnitude)//2]
179peak_indices = np.where(positive_magnitude > 0.1)[0]
180peak_freqs = positive_freqs[peak_indices]
181
182print(f"Sampling frequency: {fs} Hz")
183print(f"Signal components: 5Hz, 15Hz, 25Hz")
184print(f"\nDetected peaks in FFT:")
185for freq in sorted(peak_freqs):
186 if freq > 0:
187 idx = np.argmin(np.abs(positive_freqs - freq))
188 mag = positive_magnitude[idx]
189 print(f" {freq:.1f} Hz (magnitude: {mag:.3f})")
190
191# Test 4: Parseval's theorem
192print("\n4. PARSEVAL'S THEOREM")
193print("-" * 70)
194T = 2 * np.pi
195square_wave = lambda t: np.where(np.mod(t, T) < T/2, 1, -1)
196a0, an, bn = fourier_series_coefficients(square_wave, T, n_harmonics=20)
197
198time_energy, freq_energy = parseval_theorem_check(square_wave, T, a0, an, bn)
199print(f"Energy in time domain: {time_energy:.6f}")
200print(f"Energy in frequency domain: {freq_energy:.6f}")
201print(f"Relative difference: {abs(time_energy - freq_energy)/time_energy * 100:.2f}%")
202
203# Test 5: Windowing
204print("\n5. WINDOWING EFFECTS")
205print("-" * 70)
206# Signal: single frequency
207fs = 1000
208t = np.linspace(0, 1, fs, endpoint=False)
209signal = np.sin(2*np.pi*10*t)
210
211window_types = ['rectangular', 'hanning', 'hamming', 'blackman']
212print(f"Signal: 10 Hz sine wave")
213print(f"\nSpectral leakage (sidelobe levels):")
214
215for window_type in window_types:
216 windowed_signal, window = apply_window(signal, window_type)
217 fft_windowed = np.fft.fft(windowed_signal)
218 magnitude = np.abs(fft_windowed) / len(windowed_signal)
219
220 # Find main lobe and first sidelobe
221 main_lobe = np.max(magnitude)
222 # Look away from main peak
223 main_peak_idx = np.argmax(magnitude[:len(magnitude)//2])
224 sidelobe_region = magnitude[main_peak_idx+20:len(magnitude)//2]
225 if len(sidelobe_region) > 0:
226 sidelobe = np.max(sidelobe_region)
227 sidelobe_db = 20 * np.log10(sidelobe / main_lobe)
228 print(f" {window_type:12s}: {sidelobe_db:6.1f} dB")
229
230# Test 6: Filtering
231print("\n6. LOWPASS FILTERING")
232print("-" * 70)
233# Signal with multiple frequencies
234fs = 1000
235t = np.linspace(0, 1, fs, endpoint=False)
236signal = (np.sin(2*np.pi*5*t) +
237 np.sin(2*np.pi*50*t) +
238 np.sin(2*np.pi*120*t))
239
240print(f"Original signal: 5Hz + 50Hz + 120Hz")
241
242# Apply lowpass filter at 30Hz
243cutoff = 30
244filtered_signal, freqs, filter_mask = lowpass_filter(signal, cutoff, fs)
245
246# Check frequency content
247fft_original = np.fft.fft(signal)
248fft_filtered = np.fft.fft(filtered_signal)
249
250print(f"Lowpass filter cutoff: {cutoff} Hz")
251print(f"\nFrequency content after filtering:")
252
253for test_freq in [5, 50, 120]:
254 idx = np.argmin(np.abs(freqs[:len(freqs)//2] - test_freq))
255 original_mag = np.abs(fft_original[idx])
256 filtered_mag = np.abs(fft_filtered[idx])
257 attenuation = 20 * np.log10(filtered_mag / original_mag) if filtered_mag > 1e-6 else -100
258 print(f" {test_freq:3d} Hz: {attenuation:6.1f} dB")
259
260# Test 7: Gibbs phenomenon
261print("\n7. GIBBS PHENOMENON")
262print("-" * 70)
263T = 2 * np.pi
264square_wave = lambda t: np.where(np.mod(t, T) < T/2, 1, -1)
265
266print("Square wave reconstruction with different harmonics:")
267for n_harm in [5, 10, 50]:
268 a0, an, bn = fourier_series_coefficients(square_wave, T, n_harmonics=n_harm)
269
270 # Reconstruct at discontinuity
271 t_disc = np.array([T/2 - 0.01, T/2, T/2 + 0.01])
272 reconstructed = reconstruct_from_fourier(t_disc, T, a0, an, bn)
273
274 overshoot = (reconstructed[0] - 1) / 1 * 100 # Before discontinuity
275 print(f" n={n_harm:3d} harmonics: overshoot ≈ {abs(overshoot):.1f}%")
276
277print(f"\nGibbs overshoot approaches ~9% as n→∞")
278
279# Visualization
280if HAS_MATPLOTLIB:
281 print("\n" + "=" * 70)
282 print("VISUALIZATIONS")
283 print("=" * 70)
284
285 fig = plt.figure(figsize=(15, 12))
286
287 # Plot 1: Square wave Fourier series
288 ax1 = plt.subplot(3, 3, 1)
289 T = 2 * np.pi
290 t_plot = np.linspace(0, 2*T, 1000)
291 square_wave = lambda t: np.where(np.mod(t, T) < T/2, 1, -1)
292
293 ax1.plot(t_plot, square_wave(t_plot), 'k-', linewidth=2, label='Original', alpha=0.3)
294
295 for n_harm in [1, 3, 10]:
296 a0, an, bn = fourier_series_coefficients(square_wave, T, n_harmonics=n_harm)
297 reconstructed = reconstruct_from_fourier(t_plot, T, a0, an, bn)
298 ax1.plot(t_plot, reconstructed, label=f'n={n_harm}', linewidth=1.5)
299
300 ax1.set_xlabel('t')
301 ax1.set_ylabel('f(t)')
302 ax1.set_title('Fourier Series: Square Wave')
303 ax1.legend()
304 ax1.grid(True, alpha=0.3)
305 ax1.set_ylim(-1.5, 1.5)
306
307 # Plot 2: Sawtooth wave
308 ax2 = plt.subplot(3, 3, 2)
309 sawtooth = lambda t: (np.mod(t, T) / T) * 2 - 1
310
311 ax2.plot(t_plot, sawtooth(t_plot), 'k-', linewidth=2, label='Original', alpha=0.3)
312
313 for n_harm in [1, 5, 20]:
314 a0, an, bn = fourier_series_coefficients(sawtooth, T, n_harmonics=n_harm)
315 reconstructed = reconstruct_from_fourier(t_plot, T, a0, an, bn)
316 ax2.plot(t_plot, reconstructed, label=f'n={n_harm}', linewidth=1.5)
317
318 ax2.set_xlabel('t')
319 ax2.set_ylabel('f(t)')
320 ax2.set_title('Fourier Series: Sawtooth Wave')
321 ax2.legend()
322 ax2.grid(True, alpha=0.3)
323
324 # Plot 3: FFT spectrum
325 ax3 = plt.subplot(3, 3, 3)
326 fs = 1000
327 t = np.linspace(0, 1, fs, endpoint=False)
328 signal = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*15*t) + 0.3*np.sin(2*np.pi*25*t)
329 signal_noisy = signal + 0.1 * np.random.randn(len(signal))
330
331 fft_result = np.fft.fft(signal_noisy)
332 freqs = np.fft.fftfreq(len(signal_noisy), 1/fs)
333 magnitude = np.abs(fft_result) / len(signal_noisy)
334
335 ax3.plot(freqs[:len(freqs)//2], magnitude[:len(magnitude)//2], 'b-', linewidth=1.5)
336 ax3.set_xlabel('Frequency (Hz)')
337 ax3.set_ylabel('Magnitude')
338 ax3.set_title('FFT Spectrum: 5Hz + 15Hz + 25Hz')
339 ax3.grid(True, alpha=0.3)
340 ax3.set_xlim(0, 50)
341
342 # Plot 4: Window functions
343 ax4 = plt.subplot(3, 3, 4)
344 N = 256
345 signal_unit = np.ones(N)
346
347 for window_type in ['rectangular', 'hanning', 'hamming', 'blackman']:
348 _, window = apply_window(signal_unit, window_type)
349 ax4.plot(window, label=window_type, linewidth=1.5)
350
351 ax4.set_xlabel('Sample')
352 ax4.set_ylabel('Amplitude')
353 ax4.set_title('Window Functions')
354 ax4.legend()
355 ax4.grid(True, alpha=0.3)
356
357 # Plot 5: Windowing effect on spectrum
358 ax5 = plt.subplot(3, 3, 5)
359 fs = 1000
360 t = np.linspace(0, 1, fs, endpoint=False)
361 signal = np.sin(2*np.pi*10*t)
362
363 for window_type in ['rectangular', 'hanning', 'blackman']:
364 windowed_signal, _ = apply_window(signal, window_type)
365 fft_windowed = np.fft.fft(windowed_signal)
366 freqs = np.fft.fftfreq(len(windowed_signal), 1/fs)
367 magnitude_db = 20 * np.log10(np.abs(fft_windowed) / len(windowed_signal) + 1e-10)
368
369 ax5.plot(freqs[:len(freqs)//2], magnitude_db[:len(magnitude_db)//2],
370 label=window_type, linewidth=1.5, alpha=0.7)
371
372 ax5.set_xlabel('Frequency (Hz)')
373 ax5.set_ylabel('Magnitude (dB)')
374 ax5.set_title('Spectral Leakage with Different Windows')
375 ax5.legend()
376 ax5.grid(True, alpha=0.3)
377 ax5.set_xlim(0, 50)
378 ax5.set_ylim(-80, 10)
379
380 # Plot 6: Filtering demonstration
381 ax6 = plt.subplot(3, 3, 6)
382 fs = 1000
383 t = np.linspace(0, 0.2, 200)
384 signal = (np.sin(2*np.pi*5*t) +
385 np.sin(2*np.pi*50*t) +
386 np.sin(2*np.pi*120*t))
387
388 filtered_signal, _, _ = lowpass_filter(signal, 30, fs)
389
390 ax6.plot(t, signal, 'b-', linewidth=1.5, alpha=0.5, label='Original')
391 ax6.plot(t, filtered_signal, 'r-', linewidth=2, label='Filtered (30Hz)')
392 ax6.set_xlabel('Time (s)')
393 ax6.set_ylabel('Amplitude')
394 ax6.set_title('Lowpass Filtering')
395 ax6.legend()
396 ax6.grid(True, alpha=0.3)
397
398 # Plot 7: Frequency spectrum before/after filter
399 ax7 = plt.subplot(3, 3, 7)
400 t_full = np.linspace(0, 1, fs, endpoint=False)
401 signal_full = (np.sin(2*np.pi*5*t_full) +
402 np.sin(2*np.pi*50*t_full) +
403 np.sin(2*np.pi*120*t_full))
404 filtered_full, freqs, _ = lowpass_filter(signal_full, 30, fs)
405
406 fft_orig = np.fft.fft(signal_full)
407 fft_filt = np.fft.fft(filtered_full)
408 mag_orig = np.abs(fft_orig) / len(signal_full)
409 mag_filt = np.abs(fft_filt) / len(filtered_full)
410
411 ax7.plot(freqs[:len(freqs)//2], mag_orig[:len(mag_orig)//2],
412 'b-', linewidth=1.5, alpha=0.5, label='Original')
413 ax7.plot(freqs[:len(freqs)//2], mag_filt[:len(mag_filt)//2],
414 'r-', linewidth=2, label='Filtered')
415 ax7.axvline(x=30, color='k', linestyle='--', alpha=0.5, label='Cutoff')
416 ax7.set_xlabel('Frequency (Hz)')
417 ax7.set_ylabel('Magnitude')
418 ax7.set_title('Frequency Domain: Before/After Filter')
419 ax7.legend()
420 ax7.grid(True, alpha=0.3)
421 ax7.set_xlim(0, 150)
422
423 # Plot 8: Gibbs phenomenon
424 ax8 = plt.subplot(3, 3, 8)
425 T = 2 * np.pi
426 t_gibbs = np.linspace(T/2 - 0.5, T/2 + 0.5, 500)
427 square_wave = lambda t: np.where(np.mod(t, T) < T/2, 1, -1)
428
429 ax8.plot(t_gibbs, square_wave(t_gibbs), 'k-', linewidth=3, label='Original', alpha=0.3)
430
431 for n_harm in [5, 10, 50]:
432 a0, an, bn = fourier_series_coefficients(square_wave, T, n_harmonics=n_harm)
433 reconstructed = reconstruct_from_fourier(t_gibbs, T, a0, an, bn)
434 ax8.plot(t_gibbs, reconstructed, label=f'n={n_harm}', linewidth=1.5)
435
436 ax8.axhline(y=1.09, color='r', linestyle='--', alpha=0.5, linewidth=1)
437 ax8.set_xlabel('t')
438 ax8.set_ylabel('f(t)')
439 ax8.set_title('Gibbs Phenomenon (9% overshoot)')
440 ax8.legend()
441 ax8.grid(True, alpha=0.3)
442 ax8.set_ylim(-1.3, 1.3)
443
444 # Plot 9: Parseval's theorem
445 ax9 = plt.subplot(3, 3, 9)
446 harmonics = range(1, 51)
447 time_energies = []
448 freq_energies = []
449
450 for n_harm in harmonics:
451 a0, an, bn = fourier_series_coefficients(square_wave, T, n_harmonics=n_harm)
452 time_e, freq_e = parseval_theorem_check(square_wave, T, a0, an, bn)
453 time_energies.append(time_e)
454 freq_energies.append(freq_e)
455
456 ax9.plot(harmonics, time_energies, 'b-', linewidth=2, label='Time domain')
457 ax9.plot(harmonics, freq_energies, 'r--', linewidth=2, label='Frequency domain')
458 ax9.set_xlabel('Number of harmonics')
459 ax9.set_ylabel('Energy')
460 ax9.set_title("Parseval's Theorem Verification")
461 ax9.legend()
462 ax9.grid(True, alpha=0.3)
463
464 plt.tight_layout()
465 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/06_fourier.png', dpi=150)
466 print("Saved visualization: 06_fourier.png")
467 plt.close()
468
469print("\n" + "=" * 70)
470print("DEMONSTRATION COMPLETE")
471print("=" * 70)