# 用Python代码演示傅里叶变换和拉普拉斯变换：从理论到实战 如果你曾经在信号处理、控制系统或者电路分析的课程中，被那些抽象的积分公式和变换表搞得晕头转向，那么这篇文章就是为你准备的。我们不再仅仅停留在公式推导和理论证明的层面，而是直接动手，用Python代码将这两种强大的数学工具“可视化”和“可操作化”。无论是分析一段音频信号的频率成分，还是判断一个电路系统的稳定性，你都将看到，NumPy和SciPy这两个库如何让这些看似高深的理论变得触手可及。本文面向的是有一定Python基础，并希望将数学工具应用于实际工程问题的编程爱好者和学生。我们将从最基础的信号生成开始，一步步深入到频谱泄露、系统稳定性分析等核心实战问题，并提供可直接在Jupyter Notebook中运行的代码片段。 ## 1. 环境准备与基础概念速览 在开始编码之前，我们需要确保手头有趁手的工具。对于科学计算和信号处理，Python的生态系统提供了无与伦比的便利。我们将主要依赖`numpy`进行数值计算和数组操作，依赖`scipy`中的`fft`和`signal`模块进行快速傅里叶变换和拉普拉斯变换相关的系统分析，依赖`matplotlib`进行结果的可视化。如果你使用Anaconda发行版，这些库通常已经预装。如果没有，可以通过`pip install numpy scipy matplotlib`一键安装。 傅里叶变换和拉普拉斯变换，本质上都是将信号从一个域（通常是时域）映射到另一个域（频域或复频域）的数学工具。简单来说： * **傅里叶变换 (Fourier Transform)**：它告诉我们一个信号是由哪些不同频率的正弦波组成的，以及这些正弦波的“强度”（幅度）和“起始位置”（相位）是怎样的。它处理的是**绝对可积**或**能量有限**的信号，特别擅长分析稳态的周期或非周期信号。 * **拉普拉斯变换 (Laplace Transform)**：可以看作是傅里叶变换的“升级版”。它通过引入一个衰减因子（复指数中的实部），将一些不满足绝对可积条件（比如指数增长）的信号也变得“可变换”。这使得它不仅能分析频率特性，还能**分析系统的瞬态响应和稳定性**，在控制理论和电路分析中不可或缺。 下面的表格快速对比了二者的核心区别： | 特性 | 傅里叶变换 (FT) | 拉普拉斯变换 (LT) | | :--- | :--- | :--- | | **自变量** | 实频率 `ω` (或 `f`) | 复频率 `s = σ + jω` | | **变换域** | 频域 (Frequency Domain) | 复频域 (s-Domain) | | **核心思想** | 信号分解为正弦波 | 信号分解为复指数函数 | | **适用信号** | 能量有限/绝对可积信号 | 更广泛，尤其是指数阶信号 | | **主要应用** | 频谱分析、滤波、调制解调 | 系统稳定性分析、微分方程求解、瞬态响应 | | **与系统的关系** | 频率响应 `H(jω)` | 系统函数/传递函数 `H(s)` | > 提示：在Python的`scipy.signal`模块中，我们通常不直接计算连续的拉普拉斯变换，而是处理其离散化的“亲戚”——Z变换，或者直接使用传递函数（`TransferFunction`）对象进行系统分析。对于连续信号的拉普拉斯变换，我们更多是进行符号计算（可用`sympy`库），但本文侧重于数值计算和工程应用。 ## 2. 实战傅里叶变换：从简单信号到音频频谱分析 让我们从一个最简单的例子开始：生成一个包含多个频率分量的合成信号，然后用傅里叶变换看看它的“内在成分”。 ### 2.1 生成并观察一个合成信号 假设我们有一个信号，它由50Hz、120Hz和一个高频的随机噪声组成。我们将用代码生成它并绘制其时域波形。 ```python import numpy as np import matplotlib.pyplot as plt from scipy.fft import fft, fftfreq # 设置采样参数 Fs = 1000 # 采样频率，1000 Hz T = 1.0 # 信号总时长，1秒 N = int(Fs * T) # 采样点总数 t = np.linspace(0.0, T, N, endpoint=False) # 时间向量 # 生成信号：两个正弦波 + 噪声 freq1, amp1 = 50.0, 0.7 freq2, amp2 = 120.0, 1.0 signal = amp1 * np.sin(2.0 * np.pi * freq1 * t) + \ amp2 * np.sin(2.0 * np.pi * freq2 * t) # 添加一些随机噪声 np.random.seed(42) # 固定随机种子以便复现 noise = 0.5 * np.random.randn(N) signal_noisy = signal + noise # 绘制时域信号 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6)) ax1.plot(t, signal, 'b-', alpha=0.7, label='纯净信号') ax1.set_xlabel('时间 [秒]') ax1.set_ylabel('幅度') ax1.set_title('纯净信号时域图 (50Hz 和 120Hz)') ax1.legend() ax1.grid(True) ax2.plot(t, signal_noisy, 'r-', alpha=0.7, label='含噪信号') ax2.set_xlabel('时间 [秒]') ax2.set_ylabel('幅度') ax2.set_title('含噪信号时域图') ax2.legend() ax2.grid(True) plt.tight_layout() plt.show() ``` 运行这段代码，你会看到两个时域图。在含噪信号的图中，50Hz和120Hz的正弦波已经被噪声严重污染，几乎无法从时域波形中直接分辨出具体的频率成分。这就是傅里叶变换大显身手的时候。 ### 2.2 执行FFT并分析频谱 接下来，我们对这个含噪信号进行快速傅里叶变换（FFT），看看在频域里是什么样子。 ```python # 执行FFT yf = fft(signal_noisy) # 结果是复数数组，包含幅度和相位信息 xf = fftfreq(N, 1/Fs)[:N//2] # 生成正频率部分对应的频率轴 # 计算双边谱和单边谱的幅度 amplitude_spectrum = 2.0/N * np.abs(yf[:N//2]) # 单边谱幅度 # 绘制频谱图 plt.figure(figsize=(10, 4)) plt.plot(xf, amplitude_spectrum, 'g-') plt.xlabel('频率 [Hz]') plt.ylabel('幅度') plt.title('信号的单边幅度频谱') plt.grid(True) plt.xlim([0, 200]) # 聚焦在我们关心的频率范围 plt.axvline(x=freq1, color='b', linestyle='--', alpha=0.5, label=f'{freq1} Hz') plt.axvline(x=freq2, color='r', linestyle='--', alpha=0.5, label=f'{freq2} Hz') plt.legend() plt.show() ``` 在生成的频谱图上，你应该能在50Hz和120Hz附近看到两个清晰的尖峰，这正是我们合成信号时加入的两个正弦波的频率。噪声的能量则广泛分布在所有频率上，表现为频谱底部的“基底”。FFT成功地将时域中混叠在一起的成分，在频域中清晰地分离开来。 ### 2.3 深入探讨：频谱泄露与窗函数 在实际应用中，我们几乎不可能处理无限长的信号。我们只能对信号的一段有限时长（即一个“时间窗”）进行FFT。这会导致一个经典问题：**频谱泄露**。 假设我们生成一个频率不是采样率整数倍的信号，例如52.3Hz。 ```python # 生成一个非整数倍频的信号 freq_leak = 52.3 signal_leak = np.sin(2.0 * np.pi * freq_leak * t) # 不加窗直接进行FFT yf_leak_nowin = fft(signal_leak) amp_leak_nowin = 2.0/N * np.abs(yf_leak_nowin[:N//2]) # 应用汉宁窗 (Hanning Window) 后再进行FFT window = np.hanning(N) signal_leak_windowed = signal_leak * window yf_leak_win = fft(signal_leak_windowed) amp_leak_win = 2.0/N * np.abs(yf_leak_win[:N//2]) # 绘制对比 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4)) ax1.plot(xf, amp_leak_nowin, 'b-') ax1.set_xlabel('频率 [Hz]') ax1.set_ylabel('幅度') ax1.set_title('无窗函数 - 频谱泄露严重') ax1.set_xlim([40, 65]) ax1.grid(True) ax2.plot(xf, amp_leak_win, 'r-') ax2.set_xlabel('频率 [Hz]') ax2.set_ylabel('幅度') ax2.set_title('应用汉宁窗后 - 泄露被抑制') ax2.set_xlim([40, 65]) ax2.grid(True) plt.tight_layout() plt.show() ``` 你会观察到，在不加窗的情况下，52.3Hz信号的频谱能量“泄露”到了相邻的频率点上，形成了一个主瓣较宽、旁瓣明显的谱峰。而应用了汉宁窗之后，频谱泄露被显著抑制，主瓣变得更集中，尽管代价是主瓣宽度略有增加，频率分辨率稍有下降。 > 注意：选择窗函数是一种权衡。矩形窗（即不加窗）频率分辨率最高但泄露最严重；汉宁窗、汉明窗能有效抑制旁瓣，减少泄露，但会加宽主瓣。在实际的音频或振动信号分析中，根据信号特性选择合适的窗函数是获得准确频谱的关键一步。 ## 3. 实战拉普拉斯变换：系统稳定性分析与电路仿真 拉普拉斯变换在工程中更常见的用法是分析线性时不变系统的特性，尤其是**稳定性**。一个系统的传递函数 `H(s)` 的极点（分母为零的点）在复平面上的位置，直接决定了系统的稳定性。 ### 3.1 定义系统并分析其极点 让我们用`scipy.signal`来定义一个二阶低通滤波器的传递函数，并分析其极点。 ```python import numpy as np import matplotlib.pyplot as plt from scipy import signal import control as ct # 控制库，用于更丰富的系统分析，需额外安装: pip install control # 定义一个二阶系统的传递函数 H(s) = ω_n^2 / (s^2 + 2ζω_n s + ω_n^2) # 这是一个标准的低通滤波器，常用于模拟弹簧质量阻尼系统或RLC电路 omega_n = 10.0 # 自然频率 (rad/s) zeta = 0.25 # 阻尼比 <1，代表欠阻尼系统 # 创建传递函数对象 (分子系数， 分母系数) num = [omega_n**2] # 分子多项式系数: ω_n^2 den = [1, 2*zeta*omega_n, omega_n**2] # 分母多项式系数: s^2 + 2ζω_n s + ω_n^2 sys_tf = signal.TransferFunction(num, den) print(f"传递函数: H(s) = {omega_n**2} / (s^2 + {2*zeta*omega_n}s + {omega_n**2})") print(f"系统极点: {sys_tf.poles}") # 计算并绘制系统的单位阶跃响应 t_step, y_step = signal.step(sys_tf) plt.figure(figsize=(10, 4)) plt.plot(t_step, y_step, 'b-', linewidth=2) plt.xlabel('时间 [秒]') plt.ylabel('输出') plt.title(f'系统阶跃响应 (ζ={zeta}, ω_n={omega_n})') plt.grid(True) plt.show() ``` 运行后，控制台会打印出系统的极点。对于这个欠阻尼系统（ζ<1），极点是一对实部为负的共轭复数。**极点的实部（σ）决定了系统的衰减速度，虚部（ω）决定了振荡频率**。更重要的是，**所有极点都位于复平面左半平面（实部为负）是系统稳定的充要条件**。 ### 3.2 绘制极点-零点图与频率响应 为了更直观地理解极点位置与系统行为的关系，我们可以绘制极点-零点图和伯德图。 ```python # 绘制极点-零点图 plt.figure(figsize=(5,5)) plt.axhline(y=0, color='k', linestyle='-', alpha=0.3) plt.axvline(x=0, color='k', linestyle='-', alpha=0.3) # 绘制极点 (用'x'表示) plt.plot(sys_tf.poles.real, sys_tf.poles.imag, 'rx', markersize=10, label='Poles') # 绘制零点 (用'o'表示)，此系统无有限零点 if sys_tf.zeros.size > 0: plt.plot(sys_tf.zeros.real, sys_tf.zeros.imag, 'bo', markersize=10, label='Zeros') plt.xlabel('实部 (σ)') plt.ylabel('虚部 (jω)') plt.title('系统极点-零点图') plt.legend() plt.grid(True) plt.axis('equal') # 标记稳定区域（左半平面） plt.fill_between([-20, 0], -20, 20, color='green', alpha=0.1, label='稳定区域') plt.xlim([-omega_n*1.5, omega_n*0.5]) plt.ylim([-omega_n, omega_n]) plt.legend() plt.show() # 绘制伯德图 (Bode Plot) - 幅频和相频特性 w, mag, phase = signal.bode(sys_tf) fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8)) ax1.semilogx(w, mag, 'b-', linewidth=2) # 幅频图，对数坐标 ax1.set_ylabel('幅度 [dB]') ax1.set_title('伯德图 - 幅频特性') ax1.grid(True, which='both', axis='both', linestyle='--', alpha=0.7) ax1.axvline(x=omega_n, color='r', linestyle='--', alpha=0.5, label=f'ω_n={omega_n} rad/s') ax1.legend() ax2.semilogx(w, phase, 'r-', linewidth=2) # 相频图，对数坐标 ax2.set_xlabel('频率 [rad/s]') ax2.set_ylabel('相位 [度]') ax2.set_title('伯德图 - 相频特性') ax2.grid(True, which='both', axis='both', linestyle='--', alpha=0.7) ax2.axvline(x=omega_n, color='r', linestyle='--', alpha=0.5) plt.tight_layout() plt.show() ``` 从极点-零点图可以清晰看到，极点（红叉）位于左半平面，因此系统是稳定的。伯德图则展示了系统对不同频率正弦信号的响应：在自然频率ω_n附近，系统出现谐振峰（幅度增大），相位发生快速变化。这正是欠阻尼二阶系统的典型特征。 ### 3.3 电路系统仿真示例：RLC串联电路 让我们用一个更具体的例子——RLC串联电路，来展示拉普拉斯变换在电路分析中的应用。电路方程可以通过基尔霍夫定律建立，并转化为s域的传递函数。 假设我们有一个电压源 `V_in(s)`，串联电阻R、电感L、电容C，输出为电容两端的电压 `V_out(s)`。其传递函数为： `H(s) = V_out(s) / V_in(s) = 1 / (LCs^2 + RCs + 1)` ```python # 定义电路参数 R = 100.0 # 电阻，欧姆 L = 0.1 # 电感，亨利 C = 1e-6 # 电容，法拉 # 计算传递函数系数 (标准二阶形式: ω_n^2 / (s^2 + 2ζω_n s + ω_n^2)) omega_n_circuit = 1 / np.sqrt(L * C) zeta_circuit = R / (2 * np.sqrt(L / C)) num_circuit = [omega_n_circuit**2] den_circuit = [1, 2*zeta_circuit*omega_n_circuit, omega_n_circuit**2] sys_circuit = signal.TransferFunction(num_circuit, den_circuit) print(f"RLC电路参数: R={R}Ω, L={L}H, C={C}F") print(f"计算得到: 自然频率 ω_n = {omega_n_circuit:.2f} rad/s, 阻尼比 ζ = {zeta_circuit:.3f}") print(f"电路系统极点: {sys_circuit.poles}") # 仿真电路对方波输入的响应 t_sim = np.linspace(0, 0.01, 5000) # 10ms仿真时间 # 生成一个低频方波作为输入 input_signal = 0.5 * (signal.square(2 * np.pi * 500 * t_sim) + 1) # 500Hz方波，幅度0~1V # 使用lsim进行时域仿真 t_out, output_signal, _ = signal.lsim(sys_circuit, input_signal, t_sim) # 绘制结果 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 6), sharex=True) ax1.plot(t_sim*1000, input_signal, 'b-', label='输入电压 (方波)') ax1.set_ylabel('电压 [V]') ax1.set_title('RLC串联电路仿真 (电容电压输出)') ax1.legend() ax1.grid(True) ax2.plot(t_out*1000, output_signal, 'r-', label='输出电压 (电容电压)') ax2.set_xlabel('时间 [毫秒]') ax2.set_ylabel('电压 [V]') ax2.legend() ax2.grid(True) plt.tight_layout() plt.show() ``` 运行这段代码，你会看到电路对方波输入的响应。输出波形不再是清晰的方波，而是出现了**振荡和过冲**，这正是因为电路参数（R, L, C）构成了一个欠阻尼的二阶系统。通过分析传递函数的极点，我们可以预测这种振荡行为：极点的虚部决定了振荡频率，实部的绝对值决定了振荡衰减的快慢。如果增大电阻R，阻尼比ζ会增加，振荡会减弱，最终变为过阻尼（无振荡）的响应。 ## 4. 综合案例：音频信号降噪与滤波器设计 最后，我们将傅里叶变换和基于拉普拉斯变换（s域）设计的滤波器结合起来，完成一个简单的音频信号降噪任务。 ### 4.1 加载并分析含噪音频信号 我们将使用SciPy生成一段模拟的含噪音频信号，其中包含一个我们感兴趣的纯音和宽带噪声。 ```python from scipy.io import wavfile # 生成模拟音频信号 Fs_audio = 44100 # 标准音频采样率 duration = 3.0 # 3秒 t_audio = np.linspace(0, duration, int(Fs_audio * duration), endpoint=False) # 生成目标信号：一个440Hz的A4标准音 f_target = 440.0 target_signal = 0.5 * np.sin(2 * np.pi * f_target * t_audio) # 生成噪声：低频哼声 + 高频嘶嘶声 hum_noise = 0.1 * np.sin(2 * np.pi * 60.0 * t_audio) # 60Hz电源哼声 hiss_noise = 0.05 * np.random.randn(len(t_audio)) # 白噪声 noise_signal = hum_noise + hiss_noise # 混合信号 mixed_signal = target_signal + noise_signal # 归一化到[-1, 1]范围，模拟WAV文件格式 mixed_signal_normalized = np.int16(mixed_signal / np.max(np.abs(mixed_signal)) * 32767) # 可以保存为WAV文件供试听（可选） # wavfile.write('noisy_audio.wav', Fs_audio, mixed_signal_normalized) # 绘制时域波形和频谱 fig, axes = plt.subplots(2, 2, figsize=(12, 8)) # 时域波形 (前100ms) ax = axes[0, 0] ax.plot(t_audio[:4400] * 1000, mixed_signal[:4400], 'b-', alpha=0.7) ax.set_xlabel('时间 [毫秒]') ax.set_ylabel('幅度') ax.set_title('含噪音频信号时域波形 (前100ms)') ax.grid(True) # 整体频谱 N_audio = len(mixed_signal) yf_audio = fft(mixed_signal) xf_audio = fftfreq(N_audio, 1/Fs_audio)[:N_audio//2] amp_audio = 2.0/N_audio * np.abs(yf_audio[:N_audio//2]) ax = axes[0, 1] ax.semilogx(xf_audio, 20*np.log10(amp_audio + 1e-10), 'g-') # 转换为dB ax.set_xlabel('频率 [Hz]') ax.set_ylabel('幅度 [dB]') ax.set_title('含噪音频信号全频谱') ax.set_xlim([20, 20000]) # 音频范围 ax.grid(True) ax.axvline(x=f_target, color='r', linestyle='--', alpha=0.5, label=f'{f_target} Hz') ax.axvline(x=60, color='orange', linestyle='--', alpha=0.5, label='60 Hz Hum') ax.legend() ``` ### 4.2 设计并应用数字滤波器 从频谱中我们可以看到，除了440Hz的目标信号，还有明显的60Hz低频哼声和全频带的白噪声。我们可以设计一个**带通滤波器**，只保留目标频率附近的一个窄带。 在数字信号处理中，我们通常先设计一个模拟滤波器（s域），然后通过双线性变换等方法将其转换为数字滤波器（z域）。`scipy.signal`提供了便捷的函数来完成这些步骤。 ```python # 设计一个巴特沃斯带通滤波器，中心在440Hz，带宽约100Hz lowcut = 400 # 通带下限频率 (Hz) highcut = 480 # 通带上限频率 (Hz) order = 4 # 滤波器阶数 # 设计数字滤波器 (直接使用scipy的IIR设计) sos = signal.butter(order, [lowcut, highcut], btype='band', fs=Fs_audio, output='sos') # sos是二阶节形式，数值上更稳定 # 应用滤波器 filtered_signal = signal.sosfilt(sos, mixed_signal) # 分析滤波后信号的频谱 yf_filtered = fft(filtered_signal) amp_filtered = 2.0/N_audio * np.abs(yf_filtered[:N_audio//2]) # 绘制滤波前后的频谱对比和时域对比 ax = axes[1, 0] ax.semilogx(xf_audio, 20*np.log10(amp_audio + 1e-10), 'gray', alpha=0.5, label='原始含噪') ax.semilogx(xf_audio, 20*np.log10(amp_filtered + 1e-10), 'b-', label='滤波后') ax.set_xlabel('频率 [Hz]') ax.set_ylabel('幅度 [dB]') ax.set_title('滤波前后频谱对比') ax.set_xlim([20, 2000]) # 聚焦在低频段 ax.grid(True) ax.legend() ax.axvspan(lowcut, highcut, color='green', alpha=0.1, label='通带') ax = axes[1, 1] ax.plot(t_audio[:4400] * 1000, mixed_signal[:4400], 'gray', alpha=0.5, label='原始含噪') ax.plot(t_audio[:4400] * 1000, filtered_signal[:4400], 'r-', label='滤波后') ax.set_xlabel('时间 [毫秒]') ax.set_ylabel('幅度') ax.set_title('滤波前后时域波形对比 (前100ms)') ax.grid(True) ax.legend() plt.tight_layout() plt.show() # 计算并打印信噪比改善的简单估计 def estimate_snr(signal, noise_band_low=300, noise_band_high=600, target_freq=440, bw=20): """一个简单的SNR估计函数（仅用于演示）""" # 计算目标频带能量 (假设目标在target_freq附近±bw Hz) idx_target = np.where((xf_audio > target_freq - bw) & (xf_audio < target_freq + bw))[0] power_target = np.sum(amp_audio[idx_target]**2) # 计算噪声频带能量 (选择目标频带外的一个区域) idx_noise = np.where((xf_audio > noise_band_low) & (xf_audio < noise_band_high) & ~((xf_audio > target_freq - bw) & (xf_audio < target_freq + bw)))[0] power_noise = np.sum(amp_audio[idx_noise]**2) if len(idx_noise) > 0 else 1e-10 return 10 * np.log10(power_target / power_noise) snr_before = estimate_snr(mixed_signal) snr_after = estimate_snr(filtered_signal) print(f"估计的信噪比改善: {snr_after - snr_before:.2f} dB") ``` 通过频谱对比图可以清晰地看到，带通滤波器几乎完全滤除了60Hz的哼声和高频白噪声，只保留了440Hz附近的信号。时域波形也变得干净、平滑了许多。这个简单的例子展示了如何将频域分析（傅里叶变换找到干扰频率）与系统设计（基于拉普拉斯/模拟滤波器理论设计滤波器）结合起来解决实际问题。 在实际项目中，你可能会遇到更复杂的噪声，如周期性脉冲噪声、非平稳噪声等，这时可能需要更高级的技术，如自适应滤波、小波变换或深度学习降噪。但理解并掌握FFT和滤波器设计这两大基石，将为学习这些高级方法打下坚实的基础。