1#!/usr/bin/env python3
2"""
3Continuous Fourier Transform (CTFT) - Numerical Approximation
4=============================================================
5
6The Continuous-Time Fourier Transform (CTFT) of a signal x(t) is defined as:
7
8 X(f) = integral_{-inf}^{inf} x(t) * exp(-j*2*pi*f*t) dt
9
10We numerically approximate this integral using the trapezoidal rule over a
11finite time window, sampled at a high rate.
12
13Topics Covered:
14 1. CTFT of common signals: rectangular pulse, Gaussian, exponential decay
15 2. Fourier transform properties: linearity, time shift, frequency shift
16 3. Parseval's theorem: energy conservation in time and frequency domains
17 4. Duality: rect(t) <--> sinc(f) pair
18 5. Bandwidth and energy spectral density
19
20Author: Educational example for Signal Processing
21License: MIT
22"""
23
24import numpy as np
25import matplotlib.pyplot as plt
26from scipy.signal import gausspulse
27
28
29# =============================================================================
30# 1. Numerical CTFT computation
31# =============================================================================
32
33def ctft(x, t):
34 """
35 Numerically approximate the Continuous-Time Fourier Transform.
36
37 Uses the trapezoidal rule: X(f) ≈ sum( x(t) * exp(-j2πft) ) * dt
38
39 Args:
40 x (ndarray): Signal samples
41 t (ndarray): Uniformly spaced time array
42
43 Returns:
44 tuple: (f, X) frequency axis and complex spectrum
45 """
46 dt = t[1] - t[0]
47 N = len(t)
48 # Frequency axis: centered at 0, resolution = 1/(N*dt)
49 f = np.fft.fftfreq(N, d=dt)
50 # FFT approximates the CTFT integral when multiplied by dt
51 X = np.fft.fft(x) * dt
52 # Shift zero-frequency component to center for display
53 f = np.fft.fftshift(f)
54 X = np.fft.fftshift(X)
55 return f, X
56
57
58# =============================================================================
59# 2. Common signal definitions
60# =============================================================================
61
62def rectangular_pulse(t, width=1.0, center=0.0):
63 """
64 Rectangular pulse of given width centered at 'center'.
65
66 rect(t/T) = 1 if |t - center| <= T/2
67 0 otherwise
68
69 CTFT: T * sinc(f*T) where sinc(x) = sin(pi*x)/(pi*x)
70 """
71 return np.where(np.abs(t - center) <= width / 2, 1.0, 0.0)
72
73
74def gaussian_signal(t, sigma=0.5, center=0.0):
75 """
76 Gaussian pulse: x(t) = exp(-t^2 / (2*sigma^2))
77
78 The Fourier transform of a Gaussian is also a Gaussian:
79 CTFT: sigma*sqrt(2*pi) * exp(-2*pi^2*sigma^2*f^2)
80
81 Key property: time-bandwidth product is constant (sigma * sigma_f = 1/(2*pi))
82 """
83 return np.exp(-((t - center) ** 2) / (2 * sigma ** 2))
84
85
86def exponential_decay(t, alpha=2.0):
87 """
88 One-sided exponential decay: x(t) = exp(-alpha*t) * u(t)
89
90 CTFT: 1 / (alpha + j*2*pi*f)
91 Magnitude: 1 / sqrt(alpha^2 + (2*pi*f)^2)
92 """
93 return np.where(t >= 0, np.exp(-alpha * t), 0.0)
94
95
96# =============================================================================
97# 3. Main demonstration functions
98# =============================================================================
99
100def demo_common_signals():
101 """Compute and plot CTFT for three common signals."""
102 print("=" * 60)
103 print("Demo 1: CTFT of Common Signals")
104 print("=" * 60)
105
106 # Time axis: long window, high sampling rate for good approximation
107 dt = 0.005
108 t = np.arange(-10, 10, dt)
109
110 # --- Rectangular pulse (width = 1 s) ---
111 T_rect = 1.0
112 x_rect = rectangular_pulse(t, width=T_rect)
113 f_rect, X_rect = ctft(x_rect, t)
114 # Analytical: T * sinc(f*T)
115 with np.errstate(divide='ignore', invalid='ignore'):
116 X_rect_analytical = T_rect * np.sinc(f_rect * T_rect)
117
118 # --- Gaussian (sigma = 0.4 s) ---
119 sigma = 0.4
120 x_gauss = gaussian_signal(t, sigma=sigma)
121 f_gauss, X_gauss = ctft(x_gauss, t)
122 # Analytical: sigma*sqrt(2*pi)*exp(-2*pi^2*sigma^2*f^2)
123 X_gauss_analytical = sigma * np.sqrt(2 * np.pi) * np.exp(
124 -2 * np.pi ** 2 * sigma ** 2 * f_gauss ** 2
125 )
126
127 # --- Exponential decay (alpha = 3) ---
128 alpha = 3.0
129 x_exp = exponential_decay(t, alpha=alpha)
130 f_exp, X_exp = ctft(x_exp, t)
131 # Analytical magnitude: 1/sqrt(alpha^2 + (2*pi*f)^2)
132 X_exp_analytical = 1.0 / np.sqrt(alpha ** 2 + (2 * np.pi * f_exp) ** 2)
133
134 # Report bandwidths (3 dB point)
135 fmax = 5.0 # display limit
136 mask = np.abs(f_rect) <= fmax
137
138 print(f" Rect pulse (T={T_rect}s): first null at f = {1/T_rect:.2f} Hz")
139 print(f" Gaussian (sigma={sigma}s): spectral sigma = {1/(2*np.pi*sigma):.3f} Hz")
140 print(f" Exp decay (alpha={alpha}): half-power bandwidth = {alpha/(2*np.pi):.3f} Hz")
141
142 # Plot
143 fig, axes = plt.subplots(3, 2, figsize=(12, 10))
144 fig.suptitle("CTFT of Common Signals", fontsize=14, fontweight='bold')
145
146 signals = [
147 ("Rectangular Pulse", t, x_rect, f_rect, X_rect, X_rect_analytical, "b"),
148 ("Gaussian Pulse", t, x_gauss, f_gauss, X_gauss, X_gauss_analytical, "g"),
149 ("Exponential Decay", t, x_exp, f_exp, X_exp, X_exp_analytical, "r"),
150 ]
151
152 for i, (name, ti, xi, fi, Xi, Xi_anal, color) in enumerate(signals):
153 mask = np.abs(fi) <= fmax
154 # Time domain
155 axes[i, 0].plot(ti, xi, color=color)
156 axes[i, 0].set_xlim(-3, 5)
157 axes[i, 0].set_title(f"{name}: Time Domain")
158 axes[i, 0].set_xlabel("Time (s)")
159 axes[i, 0].set_ylabel("Amplitude")
160 axes[i, 0].grid(True, alpha=0.3)
161 # Frequency domain: magnitude
162 axes[i, 1].plot(fi[mask], np.abs(Xi[mask]), color=color, label='Numerical', lw=2)
163 axes[i, 1].plot(fi[mask], np.abs(Xi_anal[mask]), 'k--', label='Analytical', lw=1.5)
164 axes[i, 1].set_title(f"{name}: |X(f)|")
165 axes[i, 1].set_xlabel("Frequency (Hz)")
166 axes[i, 1].set_ylabel("|X(f)|")
167 axes[i, 1].legend(fontsize=8)
168 axes[i, 1].grid(True, alpha=0.3)
169
170 plt.tight_layout()
171 plt.savefig('/tmp/04_ctft_signals.png', dpi=150)
172 print(" Saved plot to /tmp/04_ctft_signals.png")
173 plt.show()
174
175
176def demo_fourier_properties():
177 """
178 Illustrate key Fourier transform properties:
179 - Linearity: a*x1(t) + b*x2(t) <--> a*X1(f) + b*X2(f)
180 - Time shift: x(t - t0) <--> X(f)*exp(-j2*pi*f*t0)
181 - Freq shift: x(t)*exp(j2*pi*f0*t)<--> X(f - f0) (modulation)
182 """
183 print("\n" + "=" * 60)
184 print("Demo 2: Fourier Transform Properties")
185 print("=" * 60)
186
187 dt = 0.005
188 t = np.arange(-10, 10, dt)
189 fmax = 6.0
190
191 # Base signal: Gaussian
192 sigma = 0.3
193 x1 = gaussian_signal(t, sigma=sigma, center=0.0)
194 x2 = rectangular_pulse(t, width=0.8, center=0.0)
195
196 f, X1 = ctft(x1, t)
197 _, X2 = ctft(x2, t)
198
199 # --- Linearity ---
200 a, b = 0.7, 0.5
201 x_linear = a * x1 + b * x2
202 _, X_linear = ctft(x_linear, t)
203 X_linear_theory = a * X1 + b * X2
204
205 # --- Time shift: shift x1 by t0 seconds ---
206 t0 = 1.5
207 x_shifted = gaussian_signal(t, sigma=sigma, center=t0)
208 _, X_shifted = ctft(x_shifted, t)
209 # Theory: X1(f) * exp(-j2*pi*f*t0)
210 X_shifted_theory = X1 * np.exp(-1j * 2 * np.pi * f * t0)
211
212 # --- Frequency shift (modulation): multiply x1 by cos(2*pi*f0*t) ---
213 f0 = 3.0
214 x_modulated = x1 * np.cos(2 * np.pi * f0 * t)
215 _, X_modulated = ctft(x_modulated, t)
216 # Theory: 0.5*(X1(f-f0) + X1(f+f0)), approximated by shifting X1
217 print(f" Linearity error: {np.max(np.abs(X_linear - X_linear_theory)):.4e}")
218 print(f" Time-shift error: {np.max(np.abs(X_shifted - X_shifted_theory)):.4e}")
219
220 mask = np.abs(f) <= fmax
221 fig, axes = plt.subplots(3, 2, figsize=(12, 9))
222 fig.suptitle("Fourier Transform Properties", fontsize=14, fontweight='bold')
223
224 # Linearity
225 axes[0, 0].plot(t, x_linear)
226 axes[0, 0].set_xlim(-3, 3)
227 axes[0, 0].set_title(f"Linearity: {a}·x₁(t) + {b}·x₂(t)")
228 axes[0, 0].set_xlabel("Time (s)")
229 axes[0, 0].grid(True, alpha=0.3)
230 axes[0, 1].plot(f[mask], np.abs(X_linear[mask]), label='Numerical', lw=2)
231 axes[0, 1].plot(f[mask], np.abs(X_linear_theory[mask]), 'r--', label='a·X₁+b·X₂', lw=1.5)
232 axes[0, 1].set_title("Spectrum: Linearity verified")
233 axes[0, 1].legend(fontsize=8)
234 axes[0, 1].grid(True, alpha=0.3)
235
236 # Time shift
237 axes[1, 0].plot(t, x1, label='x₁(t)', alpha=0.7)
238 axes[1, 0].plot(t, x_shifted, label=f'x₁(t-{t0})', alpha=0.7)
239 axes[1, 0].set_xlim(-3, 5)
240 axes[1, 0].set_title(f"Time Shift: t₀ = {t0} s")
241 axes[1, 0].legend()
242 axes[1, 0].grid(True, alpha=0.3)
243 # Phase of X_shifted should be linear: -2*pi*f*t0
244 phase_numerical = np.angle(X_shifted[mask])
245 phase_theory = np.angle(X_shifted_theory[mask])
246 axes[1, 1].plot(f[mask], phase_numerical, label='Numerical phase')
247 axes[1, 1].plot(f[mask], phase_theory, 'r--', label='-2πf·t₀')
248 axes[1, 1].set_title("Phase Shift from Time Delay")
249 axes[1, 1].set_xlabel("Frequency (Hz)")
250 axes[1, 1].set_ylabel("Phase (rad)")
251 axes[1, 1].legend(fontsize=8)
252 axes[1, 1].grid(True, alpha=0.3)
253
254 # Modulation
255 axes[2, 0].plot(t, x_modulated)
256 axes[2, 0].set_xlim(-3, 3)
257 axes[2, 0].set_title(f"Modulation: x₁(t)·cos(2π·{f0}·t)")
258 axes[2, 0].set_xlabel("Time (s)")
259 axes[2, 0].grid(True, alpha=0.3)
260 axes[2, 1].plot(f[mask], np.abs(X_modulated[mask]), label='Modulated', lw=2)
261 axes[2, 1].plot(f[mask], np.abs(X1[mask]) * 0.5, 'r--', label='0.5·|X₁(f)|', lw=1.5)
262 axes[2, 1].axvline(f0, color='gray', linestyle=':', alpha=0.6, label=f'f₀={f0} Hz')
263 axes[2, 1].axvline(-f0, color='gray', linestyle=':', alpha=0.6)
264 axes[2, 1].set_title("Spectrum shifts to ±f₀ (Freq. Shift)")
265 axes[2, 1].set_xlabel("Frequency (Hz)")
266 axes[2, 1].legend(fontsize=8)
267 axes[2, 1].grid(True, alpha=0.3)
268
269 plt.tight_layout()
270 plt.savefig('/tmp/04_fourier_properties.png', dpi=150)
271 print(" Saved plot to /tmp/04_fourier_properties.png")
272 plt.show()
273
274
275def demo_parseval_and_duality():
276 """
277 Demonstrate Parseval's theorem and the rect <--> sinc duality.
278
279 Parseval's theorem:
280 integral |x(t)|^2 dt = integral |X(f)|^2 df
281 (energy is conserved between time and frequency domains)
282
283 Duality:
284 If x(t) <--> X(f), then X(t) <--> x(-f)
285 For real symmetric signals: rect(t/T) <--> T*sinc(f*T)
286 and by duality: sinc(t*W) <--> (1/W)*rect(f/W)
287 """
288 print("\n" + "=" * 60)
289 print("Demo 3: Parseval's Theorem and Rect-Sinc Duality")
290 print("=" * 60)
291
292 dt = 0.002
293 t = np.arange(-15, 15, dt)
294
295 # --- Parseval's Theorem ---
296 T_rect = 1.0
297 x = rectangular_pulse(t, width=T_rect)
298 f, X = ctft(x, t)
299 df = f[1] - f[0]
300
301 energy_time = np.trapz(np.abs(x) ** 2, t)
302 energy_freq = np.trapz(np.abs(X) ** 2, f)
303
304 print(f" Parseval's Theorem (rect pulse, T={T_rect}):")
305 print(f" Time-domain energy: {energy_time:.6f}")
306 print(f" Freq-domain energy: {energy_freq:.6f}")
307 print(f" Relative error: {abs(energy_time - energy_freq)/energy_time:.2e}")
308
309 # --- Duality: sinc in time → rect in frequency ---
310 # x(t) = sinc(W*t) has Fourier transform (1/W)*rect(f/W)
311 W = 2.0 # bandwidth in Hz
312 x_sinc = np.sinc(W * t) # np.sinc(x) = sin(pi*x)/(pi*x)
313 f_sinc, X_sinc = ctft(x_sinc, t)
314
315 fmax = 6.0
316 mask = np.abs(f_sinc) <= fmax
317
318 # Plot
319 fig, axes = plt.subplots(2, 2, figsize=(12, 8))
320 fig.suptitle("Parseval's Theorem and Rect-Sinc Duality", fontsize=14, fontweight='bold')
321
322 # Energy spectral density |X(f)|^2
323 axes[0, 0].plot(t, x, 'b')
324 axes[0, 0].set_xlim(-3, 3)
325 axes[0, 0].set_title(f"rect(t/{T_rect}): Energy = {energy_time:.3f} J")
326 axes[0, 0].set_xlabel("Time (s)")
327 axes[0, 0].set_ylabel("Amplitude")
328 axes[0, 0].grid(True, alpha=0.3)
329
330 axes[0, 1].plot(f[mask], np.abs(X[mask]) ** 2, 'b')
331 axes[0, 1].fill_between(f[mask], np.abs(X[mask]) ** 2, alpha=0.3)
332 axes[0, 1].set_title(f"Energy Spectral Density |X(f)|² (∫ = {energy_freq:.3f} J)")
333 axes[0, 1].set_xlabel("Frequency (Hz)")
334 axes[0, 1].set_ylabel("|X(f)|²")
335 axes[0, 1].grid(True, alpha=0.3)
336
337 # Duality: sinc(Wt) <--> (1/W)*rect(f/W)
338 axes[1, 0].plot(t, x_sinc, 'g')
339 axes[1, 0].set_xlim(-4, 4)
340 axes[1, 0].set_title(f"sinc({W}t) in time domain")
341 axes[1, 0].set_xlabel("Time (s)")
342 axes[1, 0].set_ylabel("Amplitude")
343 axes[1, 0].grid(True, alpha=0.3)
344
345 axes[1, 1].plot(f_sinc[mask], np.abs(X_sinc[mask]), 'g', label='Numerical', lw=2)
346 # Analytical: (1/W)*rect(f/W)
347 X_sinc_analytical = (1.0 / W) * rectangular_pulse(f_sinc, width=W)
348 axes[1, 1].plot(f_sinc[mask], X_sinc_analytical[mask], 'r--', label=f'(1/{W})·rect(f/{W})', lw=1.5)
349 axes[1, 1].set_title("Duality: sinc(Wt) ↔ (1/W)·rect(f/W)")
350 axes[1, 1].set_xlabel("Frequency (Hz)")
351 axes[1, 1].set_ylabel("|X(f)|")
352 axes[1, 1].legend()
353 axes[1, 1].grid(True, alpha=0.3)
354
355 plt.tight_layout()
356 plt.savefig('/tmp/04_parseval_duality.png', dpi=150)
357 print(" Saved plot to /tmp/04_parseval_duality.png")
358 plt.show()
359
360
361if __name__ == "__main__":
362 print("Continuous Fourier Transform (CTFT) - Numerical Approximation")
363 print("=" * 60)
364 print("The CTFT integral is approximated numerically using FFT * dt.")
365 print("A fine time grid (high sampling rate) over a long window")
366 print("gives good agreement with analytical formulas.\n")
367
368 demo_common_signals()
369 demo_fourier_properties()
370 demo_parseval_and_duality()
371
372 print("\n" + "=" * 60)
373 print("Key Takeaways:")
374 print(" - CTFT ≈ FFT * dt (for uniformly sampled, windowed signals)")
375 print(" - Linearity: superposition holds in both domains")
376 print(" - Time shift introduces linear phase: exp(-j2πft₀)")
377 print(" - Modulation shifts spectrum to carrier frequency ±f₀")
378 print(" - Parseval: total energy is conserved across domains")
379 print(" - Duality: rect ↔ sinc (fundamental transform pair)")
380 print("=" * 60)