Python信号处理实战:从FFT频谱分析到滤波器设计

📅 发布时间:2026/8/31 7:06:41
Python信号处理实战:从FFT频谱分析到滤波器设计 很多人学完《信号与系统》和《数字信号处理》两门课最容易产生一种错觉公式背了不少变换算得挺熟练但面对一段真实的采集数据却不知道该用什么工具做频谱分析不知道该怎样设计一个能用的滤波器甚至不清楚为什么调制、采样这些话题总是绕不开。这个现象很普遍。问题不在于公式本身而在于课程很少把“时域信号、频域表示、线性系统响应、采样离散化”这条主线连成一个可操作的闭环。本文想做的事情很直接用 Python 和几段可复制的代码把频谱分析、滤波器设计、调制解调、采样定理这些核心概念全部落到真实波形上让你看到信号在频域里到底是什么样子系统的频率响应为什么能决定滤波结果工程上理解这些概念时应该抓住哪些要点。先给一个明确判断信号与系统这门课真正重要的不是记住拉普拉斯变换、Z 变换的每一个性质而是理解“线性时不变系统 傅里叶变换”这条骨架。只要骨架通了滤波、调制、采样都不是孤立的记忆点而是同一个频域视角下的不同应用。读完这篇文章你至少能独立跑通一个完整的信号处理小实验生成混合信号、用 FFT 看频谱、设计低通滤波器去除噪声、再用调幅与包络检波理解通信链路的基本过程。1. 这篇文章真正要解决的问题先想一个问题为什么工程师需要频谱分析真实世界里的信号几乎不会是一条干净的理想正弦波。语音里有基频和谐波振动信号里有转频和故障特征频率射频信号里载波和基带调制混在一起传感器数据里除了有效成分还有噪声和工频干扰。这些信号全部是一堆频率分量叠加在一起的结果。时域波形只能告诉你幅值怎么随时间变化却很难告诉你这个信号由哪些频率成分构成、每种成分的能量有多少。频谱分析的价值就是把时间轴上的叠加拆成频率轴上的分布。从这个角度看频谱分析解决的其实是“信号可读性”问题让模糊的复合波形变成清晰的频率线索。而频域分析之所以重要是因为它给“系统行为”提供了一个非常直观的描述方式——一个线性时不变系统对正弦信号的响应依然是同频率的正弦信号只是幅度和相位会发生改变。这句话是整个频域分析方法的基石。滤波器、调制解调、采样率选择本质上都在围绕这个规律转。这篇文章适合这几类读者正在复习《信号与系统》《数字信号处理》希望把概念和实际计算对应起来的人。做嵌入式、音频处理、传感器采集、通信相关开发需要处理采样波形和噪声问题的人。想用 Python 快速验证一个频域想法但不太确定 FFT 之后该怎么归一化、滤波器参数怎么定的人。你不必是数学高手。文中会给出最小可运行的实验代码你能理解“输入波形—频谱视图—系统处理—输出波形”这条链路就够了。2. 核心概念时域、频域、频谱与傅里叶变换2.1 时域和频域是同一个信号的两种视角时域描述的是信号瞬时值随时间的变化典型表示是 x(t) 或 x[n]。频域描述的是信号功率或幅度在频率轴上的分布典型表示是 X(f) 或 X(k)。同一个信号在这两个域里携带的信息是完全等价的区别只是人们选择从哪个角度观察。打个比方一段音乐可以看成时间轴上连续变化的声波也可以看成无数不同频率音符的叠加。时域适合观察事件发生的顺序、突变和包络频域适合判断“这个信号里到底有哪些频率成分”。后者在很多工程场景里更有决策价值。初学者最容易犯的思维错误是“时域和频域是两个不同信号”。实际上它们只是同一个信号在不同坐标系下的投影不存在谁更真实的问题。2.2 从傅里叶变换到 FFT傅里叶变换的核心思想是任何满足一定条件的信号都可以分解成不同频率、不同幅度、不同相位的正弦波之和。用公式表示为连续傅里叶变换CFTX(f) ∫ x(t) · e^(-j2πft) dt在计算机中信号是离散采样后的有限长序列所以需要使用离散傅里叶变换DFT。DFT 的计算复杂度是 O(N²)当 N 很大时并不现实。快速傅里叶变换FFT是 DFT 的高效算法复杂度降为 O(N log N)。这就是为什么工程中只要做频谱分析几乎都是调用 FFT。Python 里对应的就是 NumPy 的numpy.fft.fft和numpy.fft.fftfreq。前者计算频域复数值后者生成频率轴。2.3 线性时不变系统为什么是核心主线考察一个系统是否容易用频域分析关键看它是否满足两个性质线性和时不变。线性对输入信号进行加权叠加输出也能按同样的权重叠加。时不变输入延时一段时间输出只发生相同延时不会改变自身形态。当系统同时满足这两个条件时它就被称为线性时不变系统简称 LTI 系统。LTI 系统对复指数信号 e^(j2πft) 的输出只会改变该频率分量的幅度和相位不会产生新的频率分量。这个性质非常关键复杂的输入信号可以拆成正弦波叠加每个频率分量通过系统的行为可以用一个复数值描述整体输出就是所有分量响应的叠加。于是就有了“频率响应” H(f) 的概念。低通滤波器、高通滤波器、带通滤波器本质上都是在设计一个特定的 H(f)保留想要的频段压制不想要的频段。调制器、积分器、微分器等系统行为也都可以放到这个框架里理解。3. 环境准备Python 信号处理开发环境本文的示例基于 Python需要以下依赖Python 3.8 及以上版本。NumPy负责数组运算和 FFT。SciPy提供滤波器设计与信号处理函数。Matplotlib用于绘制时域波形和频谱图。建议在虚拟环境中安装避免污染系统环境python -m venv signal-env source signal-env/bin/activate在 Windows 下激活命令为signal-env\Scripts\activate安装依赖pip install numpy scipy matplotlib验证环境是否正常import numpy as np import scipy import matplotlib print(numpy:, np.__version__) print(scipy:, scipy.__version__) print(matplotlib:, matplotlib.__version__)如果你能看到三个版本号环境就准备好了。后续所有代码都建议在同一个脚本或 Jupyter Notebook 中运行方便对比波形和频谱。4. 采样定理连续信号离散化必须跨过的门槛真实世界里的信号是连续的计算机只能处理离散样本。把连续信号变成离散序列的过程叫采样它有一个无法回避的约束采样率必须大于信号最高频率的两倍。这个结论就是奈奎斯特采样定理。为什么是两倍因为离散采样无法区分频率为 f 的信号和频率为 f k·fs 的信号这种现象叫混叠。如果信号频率超过采样频率的一半即奈奎斯特频率 fs/2高频分量就会被折叠到低频区域原本不存在的频率出现在频谱里而且无法通过后续滤波消除。用一段代码直观感受混叠。假设真实信号是 300 Hz 的正弦波我们分别用 1000 Hz 和 400 Hz 采样import numpy as np import matplotlib.pyplot as plt f_signal 300 duration 0.02 fs_good 1000 fs_bad 400 t_good np.arange(0, duration, 1/fs_good) t_bad np.arange(0, duration, 1/fs_bad) x_good np.sin(2 * np.pi * f_signal * t_good) x_bad np.sin(2 * np.pi * f_signal * t_bad) plt.figure(figsize(10, 4)) plt.subplot(1, 2, 1) plt.plot(t_good, x_good, markero, markersize3) plt.title(ffs {fs_good} Hz) plt.xlabel(time/s) plt.subplot(1, 2, 2) plt.plot(t_bad, x_bad, markero, markersize4) plt.title(ffs {fs_bad} Hz) plt.xlabel(time/s) plt.tight_layout() plt.show()400 Hz 采样时奈奎斯特频率是 200 Hz300 Hz 的信号不满足采样定理。把样本点连起来看它看起来像一个 100 Hz 的信号这个 100 Hz 就是混叠产物计算公式是 |f - fs| |300 - 400| 100 Hz。工程上的教训很明确在 ADC 采样之前必须先确认信号带宽必要时加抗混叠滤波器保证进入采样器的信号最高频率不超过 fs/2。一旦混叠发生后级任何数字处理都无法还原真实频谱。5. 频域分析完整示例从时域波形到频谱图现在做一个最典型的实验生成一个由 50 Hz 和 120 Hz 正弦波叠加的信号再叠加随机噪声然后用 FFT 做频谱分析。import numpy as np import matplotlib.pyplot as plt fs 1000 T 1.0 N int(fs * T) t np.linspace(0, T, N, endpointFalse) f1 50 f2 120 x np.sin(2 * np.pi * f1 * t) 0.6 * np.sin(2 * np.pi * f2 * t) x x 0.3 * np.random.randn(N) X np.fft.fft(x) freqs np.fft.fftfreq(N, 1/fs) mag np.abs(X) * 2 / N half N // 2 plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(t[:300], x[:300]) plt.title(time domain) plt.xlabel(time/s) plt.ylabel(amplitude) plt.subplot(1, 2, 2) plt.plot(freqs[:half], mag[:half]) plt.title(frequency domain) plt.xlabel(frequency/Hz) plt.ylabel(magnitude) plt.xlim(0, 250) plt.tight_layout() plt.show() print(Top 5 peak frequencies:) idx np.argsort(mag[:half])[::-1][:5] for i in idx: print(f{freqs[i]:.1f} Hz : {mag[i]:.3f})关键点解释np.fft.fft(x)得到的是复数值序列每个复数包含幅度和相位信息。np.fft.fftfreq(N, 1/fs)生成与 FFT 结果对应的频率轴单位是 Hz。np.abs(X) * 2 / N是单边谱的幅度归一化。FFT 结果在正负频率上对称取正频部分时幅度需要除以 N 并乘以 2才能还原原始正弦波的幅度。前面先取了mag[:half]所以只画正频率部分。运行后频谱图会在 50 Hz 和 120 Hz 附近出现明显峰值其他频率上则是噪声底。你可以直观看到虽然时域波形被噪声污染得无法直接读出频率但频域里特征非常明显。这就是频谱分析在故障诊断、语音识别、振动监测中被广泛使用的原因。需要注意两个细节频谱泄漏当信号频率不是 FFT 分辨率的整数倍时能量会泄漏到相邻频点。分辨率为 fs/N如果 N 不够大两个相近频率可能无法分辨。窗函数加窗可以抑制频谱泄漏但会降低主瓣分辨能力。工程中在时域截断之前选择合适窗函数例如汉宁窗、布莱克曼窗是常见手段。6. 滤波器设计与应用6.1 滤波器本质是频率响应设计滤波器不是一种“魔法”它是对 LTI 系统频率响应 H(f) 的工程设计。低通滤波器让低于截止频率的分量通过衰减高于截止频率的分量高通、带通、带阻同理。从实现角度看数字滤波器分为两大类FIR 滤波器有限脉冲响应容易实现线性相位但需要较高阶数。IIR 滤波器无限脉冲响应效率高、阶数低但相位通常是非线性的。在 SciPy 中最常用的设计方法是scipy.signal.butter它能生成巴特沃斯低通滤波器的系数。巴特沃斯滤波器的特点是通带内最大平坦适合大多数工程场景。6.2 一个可复用的低通滤波器函数以下代码设计了一个 4 阶巴特沃斯低通滤波器并用filtfilt做零相位滤波from scipy.signal import butter, filtfilt def butter_lowpass_filter(data, cutoff, fs, order4): nyq 0.5 * fs normal_cutoff cutoff / nyq b, a butter(order, normal_cutoff, btypelow, analogFalse) y filtfilt(b, a, data) return y cutoff 80 x_filtered butter_lowpass_filter(x, cutoff, fs, order4)代码说明cutoff是截止频率单位是 Hz。normal_cutoff cutoff / nyq内部转换为归一化频率范围在 0 到 1 之间。butter返回分子系数 b 和分母系数 a。filtfilt对数据做正向和反向两次滤波从而抵消相位偏移输出信号没有净相位延迟。这是离线信号处理中非常实用的方法。把滤波后的时域波形和频谱画出来X_filtered np.fft.fft(x_filtered) mag_filtered np.abs(X_filtered) * 2 / N plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(t, x, alpha0.5, labelraw) plt.plot(t, x_filtered, linewidth2, labelfiltered) plt.title(time domain) plt.legend() plt.subplot(1, 2, 2) plt.plot(freqs[:half], mag[:half], alpha0.5, labelraw) plt.plot(freqs[:half], mag_filtered[:half], linewidth2, labelfiltered) plt.xlim(0, 250) plt.title(frequency domain) plt.legend() plt.tight_layout() plt.show()在这个例子中80 Hz 的截止频率可以保留 50 Hz 分量同时显著削弱 120 Hz 分量随机噪声也会被大幅压制。你可以修改 cutoff 数值观察频谱中分量幅度如何变化从而理解“设计滤波器”到底在调什么。6.3 滤波器的常见坑butter的阶数过高会导致数值不稳定尤其当截止频率很低或很高时。建议先试 4 到 8 阶。filtfilt会引入边缘效应信号两端可能出现轻微摆动长数据的边缘影响较小。实时系统中不能直接使用filtfilt因为它需要整段数据。实时场景应使用lfilter或面向流式数据的sosfilt。7. 调制与解调把信号搬到需要的频段7.1 为什么需要调制基带信号的频率通常较低例如语音主要分布在 300 Hz 到 3400 Hz 之间。如果所有设备都在同一频段直接传输天线尺寸、频谱划分、多路复用都会成为问题。调制的本质是用一个高频载波把基带信号搬到合适的频段使天线尺寸合理、不同业务占用不同频段、多路信号可以共享信道。7.2 用 Python 实现调幅与包络检波调幅AM是最直观的调制方式。它的数学形式是s(t) [1 m·x(t)] · cos(2π·fc·t)其中 x(t) 是基带信号m 是调幅系数取值在 0 到 1 之间fc 是载波频率。下面用代码完成一个完整的 AM 调制与包络检波过程from scipy.signal import butter, filtfilt fs 8000 T 0.1 t np.arange(int(fs * T)) / fs fm 100 fc 1000 m 0.8 message np.sin(2 * np.pi * fm * t) carrier np.cos(2 * np.pi * fc * t) am_signal (1 m * message) * carrier envelope np.abs(am_signal) b, a butter(4, 2 * fm / fs) demod filtfilt(b, a, envelope) plt.figure(figsize(12, 8)) plt.subplot(3, 1, 1) plt.plot(t[:400], message[:400]) plt.title(baseband signal) plt.xlabel(time/s) plt.subplot(3, 1, 2) plt.plot(t[:400], am_signal[:400]) plt.title(AM modulated signal) plt.xlabel(time/s) plt.subplot(3, 1, 3) plt.plot(t[:400], demod[:400]) plt.title(demodulated signal) plt.xlabel(time/s) plt.tight_layout() plt.show()流程解释产生 100 Hz 的基带正弦波。用 1000 Hz 载波调幅调幅系数 0.8。载波频率远高于基带频率才能在时域中看到清晰的包络。包络检波的第一步是取绝对值此时信号中保留了载波频率和基带包络需要一个低通滤波器去掉载波成分。第四行butter(4, 2 * fm / fs)的归一化截止频率设置为 200 Hz即 2·fm / fs目的是让基带 100 Hz 通过同时压制 1000 Hz 载波。如果观察 AM 信号的频谱你会看到原始基带频谱被搬移到载频 1000 Hz 两侧形成以 fc 为中心的上边带和下边带。这个观察可以更深刻地理解“调制把频谱搬移”这句话的含义。实际通信系统还涉及频分复用、双边带、单边带、正交调制等更复杂方案但背后的核心都是傅里叶变换的频率搬移特性。8. 常见问题与排查思路问题现象可能原因排查方式解决方案频谱峰值位置不对采样率或频率轴计算错误检查 fs 和 fftfreq 参数确认 fs 单位是 Hz频率轴使用 fftfreq(N, 1/fs)FFT 后出现对称的波形没有取单边谱画了正负频率全部检查绘图坐标范围或数组切片只取前 N//2 个频率点绘制单边谱信号幅度与理论幅值不符没有正确归一化幅度打印峰值和输入幅度使用 np.abs(X)*2/N 归一化直流分量不乘 2滤波后波形出现摆动边缘效应或阶数过高观察信号首尾部分增加数据长度降低阶数或使用 pad 处理边缘采样后出现不存在的低频采样率不满足奈奎斯特定理发生混叠检查信号的最高频率和 fs 关系提高采样率或添加抗混叠滤波器后再采样滤波器相位失真明显使用了 lfilter 且滤波器的相位非线性对比原波形和输出波形的时间偏移离线分析使用 filtfilt实时系统允许一定延迟时使用线性相位 FIR包络检波得到的信号毛刺多低通截止频率设置过高载波泄漏检查滤波后频谱是否仍包含 fc 分量降低截止频率但需高于基带最高频率建议排查顺序先确认采样率与最高频率关系再检查频率轴再检查幅度归一化最后检查滤波器参数。大多数问题都出在前面两步。9. 工程实践建议从实验到真实项目9.1 采样率设计要留余量奈奎斯特定理给出的是理论下限工程上通常会让采样率远高于信号最高频率常见做法是最高频率的 5 到 10 倍。这样可以降低抗混叠滤波器的设计难度也能提供更宽松的频谱空间。9.2 先看频谱再决定滤波器参数很多工程师一上来就写滤波器结果参数全靠猜。正确做法是先对原始信号做 FFT 或功率谱估计观察哪些频率是有效成分哪些是噪声和干扰再根据有效频带确定滤波器的类型和截止频率。滤波器设计必须基于对频谱的充分理解而不是经验值。9.3 区分离线分析和实时处理Filtfilt 只适合离线数据因为它需要整段数据进行双向滤波。实时系统应使用scipy.signal.lfilter或sosfilt并结合块式处理保证延迟可控。实时处理时滤波器阶数、缓冲区大小、计算耗时都需要纳入设计。9.4 记录你的信号处理参数在工程项目中采样率、滤波器类型、截止频率、阶数、窗函数这些参数必须像代码一样管理起来。建议写入配置文件并在输出结果中记录。否则两个月后回看脚本很难还原当时为什么选这些参数。9.5 警惕非线性相位IIR 滤波器阶数低、效率高但相位非线性会导致波形失真尤其对脉冲信号和通信信号影响较大。对于强调波形保真度的场景优先选择 FIR 滤波器或使用零相位滤波。对于只关心幅度谱的应用IIR 通常足够。10. 总结信号与系统的“最小认知闭环”把前面所有内容串起来可以发现一条完整的思路真实信号是时域波形傅里叶变换把它映射到频域采样定理决定了数字系统能否无失真表达连续信号频谱分析告诉你信号里有什么滤波器设计改变了频谱形状调制则利用频谱搬移让信号适应信道。每一块都不是孤立知识而是同一个频域视角的不同操作。建议你现在就动手做一个小实验生成一个包含两个相近频率的信号分别改变采样点数、采样率和滤波截止频率观察频谱和时域波形如何变化。这个实验做一遍比抄十遍公式更能建立直觉。下一步可以继续深入这些方向窗函数和频谱泄漏的定量分析。功率谱密度估计理解随机信号处理。FIR 与 IIR 滤波器的性能对比。数字下变频、正交解调、QPSK 等通信系统基础。如果你在做嵌入式采集、音频处理或传感器数据分析把这一套流程跑通后你会发现自己面对数据时不再迷茫而是能清晰地说出“先看频谱再定方案”。这篇内容可以先收藏备用遇到相关问题直接回来对照代码和思路。