Python FFT/IFFT频域滤波与四类滤波器设计实战

1. 先把问题摆清楚滤波和频域变换为什么总被绑在一起信号处理这块很多人第一次接触滤波器是在课本上看到那一堆传递函数然后被拉普拉斯变换劝退。我当年也是这样。真正让我搞明白低通、高通、带通、带阻这四类滤波器各自定位的不是推导公式而是某次处理一段带工频干扰的传感器数据——我需要把某个固定频率的“刺”拔掉同时保住低频的趋势项。那次之后低通、高通、带通、带阻在我脑子里才算各就各位。而把这件事做得又快又直观的路径就是用 Python 里的FFT 和 IFFT做频域分析、用scipy.signal做时域滤波设计两条路互相印证。这篇内容适合三类人刚学完数字信号处理课程、想把公式落成代码的学生手里有一堆采样数据、需要去噪或提特征的工程实践者以及写 Python 数据分析脚本、临时需要“把这段信号滤一下”的开发者。我会把那四个滤波器的参数怎么定、FFT 出来的频率轴怎么对上、IFFT 还原时哪里最容易翻车全部拆开讲代码可以直接抄。我先把最核心的一句话放在前面滤波器本质上是在频域里做取舍FFT 只是把这个取舍过程可视化并加速了。你理解了“哪些频率该留、哪些该扔”剩下的就是参数和工具的事。而绝大多数人卡住的地方恰恰不是滤波器公式而是几个很具体的工程细节——归一化频率算错了、双边谱没取对称、filtfilt用在了不该用的地方、IFFT 之后发现信号整体位移了。这些我在后面都会一个一个拎出来。2. 从采样率出发FFT 到底算出了什么频率轴怎么对上2.1 频率分辨率是手算出来的不是软件猜的很多人调用np.fft.fft(x)之后直接画图横轴随便标个 0 到 N然后看着一堆峰不知道对应多少赫兹。问题出在没有建立“采样率—点数—频率分辨率”这条链。假设采样率fs 1000 Hz采样点数N 1000那么这段数据在时间上覆盖 1 秒频率分辨率df fs / N 1 Hz。也就是说FFT 输出的第 k 个点对应的频率是k * df。如果你的数据是 2.5 秒、采样率 500 Hz那么N 1250df 500 / 1250 0.4 Hz。想提高分辨率只有两条路拉长采样时间或者降低采样率后者会牺牲高频信息通常不可取。这里有个反直觉的点——补零不能提高真实分辨率。你给 1000 个点补到 4096 个点再做 FFT频率轴看起来更密了但那只是频谱插值两个相邻频率成分依然分不开因为真实的时间窗长度没变。我见过不少人以为补零能“看得更清楚”结果在密集谱线场景下判断失误。import numpy as np fs 1000.0 N 1000 x np.random.randn(N) X np.fft.fft(x) freqs np.fft.fftfreq(N, d1.0 / fs) df fs / N print(f频率分辨率: {df} Hz) print(f第一个正频率点: {freqs[1]} Hz) print(f奈奎斯特频率: {fs / 2} Hz)fftfreq返回的数组前半段是 0 到接近fs/2的正频率后半段是负频率。奈奎斯特频率是能表示的最高频率等于fs/2超过这个频率的成分会被折叠混叠到低频这是采样定理的硬约束任何后处理都救不回来。所以采集端的抗混叠滤波永远排在第一位这一点我在多个项目里反复验证过。2.2 单边谱、双边谱和 fftshift 的三个坑np.fft.fft输出的是双边谱正负频率都有。对于实信号负频率部分是正频率的共轭镜像信息是冗余的。做幅频分析时我们通常只看单边谱并且要注意幅值换算单边谱的幅值要乘以 2直流分量和奈奎斯特点除外再除以 N。忘记除以 N画出来的“幅值”会随点数变化完全没法横向比较忘记乘 2非直流成分的幅值会少一半。X np.fft.fft(x) half N // 2 amp np.abs(X[:half]) / N amp[1:] * 2 # 除直流外乘2 freqs_half np.fft.fftfreq(N, 1/fs)[:half] # 或者用 rfft 一步到位 X_r np.fft.rfft(x) freqs_r np.fft.rfftfreq(N, 1/fs) amp_r np.abs(X_r) / N amp_r[1:-1] * 2np.fft.fftshift的作用是把零频搬到数组中间画图时更符合直觉。但要注意fftshift 之后频率轴的顺序也变了如果你后续要在这个数组上做掩码操作必须用fftshift后的fftfreq配合否则掩码会错位。这个错误极其隐蔽——图看起来很正常滤波结果却是乱的。我的习惯是分析画图用fftshift频域处理一律用未 shift 的原始顺序或者干脆用rfft系列省掉正负频率的心智负担。2.3 IFFT 还原前必须确认的三件事用 FFT 处理完频谱、准备ifft回时域时有三件事必须先确认。第一掩码是否保持了共轭对称性。如果你直接把负频率部分清零、只留正频率IFFT 出来的结果会是复数实部幅值只有原来的一半虚部还带着残余。第二长度是否一致。频域改了长度比如补零ifft的时候要用n参数指定正确的点数否则时域时长会变。第三是否做了必要的缩放。某些库的正逆变换归一化约定不同来回变换后幅值会差一个 N 倍务必用简单正弦做一次往返测试验证。# 正确的做法用 rfft / irfft天然处理实信号 X_r np.fft.rfft(x) f np.fft.rfftfreq(N, 1/fs) mask (f 100.0) # 低通掩码保留 100 Hz 以下 y np.fft.irfft(X_r * mask, nN) print(还原后是否为实数:, np.isrealobj(y))用rfft/irfft这一对等于把共轭对称性交给库去管我在这类频域滤波脚本里基本只用它。反过来如果你出于某种原因必须用fft/ifft那掩码一定要写成“以零频为中心对称”的形式比如mask np.abs(f) 100这样正负频率同时被保留还原出来才是实信号。这两种写法我都见过大量误用尤其是从网上抄来的代码片段往往只对一个特定长度成立换个数据就崩。3. 四类滤波器的参数怎么定从需求到 scipy 代码3.1 先明确“留什么、扔什么”再谈工具选型低通、高通、带通、带阻的定义其实一句话就能说清低通保留低频、高通保留高频、带通只保留某一段、带阻只挖掉某一段。真正的难点在于把模糊的工程需求翻译成具体的数值。比如“去掉高频噪声”高频是多高我的经验是先去画一次原始信号的频谱看噪声和有用信号的分界大概在哪个频率把这个分界作为截止频率的起点然后留出足够宽的过渡带。工具选型上Python 生态里做 IIR 滤波器设计首选scipy.signal的butter、cheby1、ellip做 FIR 用firwin或remez。巴特沃斯butter幅频响应最平滑通带内没有纹波代价是过渡带较宽切比雪夫cheby1过渡带更陡但通带会有纹波椭圆滤波器ellip最陡通带阻带都有纹波。我一般默认用巴特沃斯除非对过渡带宽度有硬性要求。原因很简单平滑的通带意味着信号波形失真小而多数去噪场景并不在意那点多出来的过渡带。滤波器类型通带特性过渡带陡度相位非线性典型适用场景巴特沃斯 butter最平坦无纹波较缓中等通用去噪、趋势提取切比雪夫I cheby1有等波纹较陡较强需要窄过渡带的场景椭圆 ellip通带阻带均有纹波最陡最强强干扰抑制、资源受限FIR firwin可做到线性相位取决于阶数可线性要求波形不失真的场景阶数是另一个关键参数。巴特沃斯的幅频响应在截止频率之后大约每阶每十倍频程衰减 20 dB。举例来说4 阶低通截止 100 Hz那么在 200 Hz 处的衰减大约是20 * 4 * log10(200/100) ≈ 24 dB也就是幅值降到约 6%。如果你需要 200 Hz 处衰减 60 dB那就得把阶数提到 10 阶左右但阶数越高数值稳定性和相位失真越成问题。一般工程上 2 到 6 阶就够用不够就考虑换更陡的滤波器类型而不是无脑堆阶数。3.2 低通和高通归一化频率是唯一的坑scipy.signal.butter最难记的一点是截止频率Wn必须是归一化到奈奎斯特频率的值范围 0 到 1。也就是说Wn fc / (fs / 2)。这个设计让同一套代码适配任意采样率但也是新手翻车最多的地方——直接传fc 100进去如果fs 1000归一化后应该是 0.2传 100 会直接报错或得到荒谬的结果。import numpy as np from scipy import signal fs 1000.0 fc_low 100.0 # 低通保留 100 Hz 以下 wn fc_low / (fs / 2) b_lp, a_lp signal.butter(N4, Wnwn, btypelow, analogFalse) # 高通保留 20 Hz 以上常用于去基线漂移 fc_high 20.0 wn_hp fc_high / (fs / 2) b_hp, a_hp signal.butter(N4, Wnwn_hp, btypehigh, analogFalse) # 查看实际频响验证设计是否符合预期 w, h signal.freqz(b_lp, a_lp, worN2048, fsfs)signal.freqz这一步千万别省。设计完滤波器立刻画一次幅频响应看清楚截止频率位置对不对、过渡带多宽、阻带衰减够不够。我见过有人把btype写错、把采样率记错全靠这一步检查抓出来。另外freqz传了fs参数之后返回的w直接就是赫兹不用再自己换算这个小细节能省不少事。高通的截止频率选择有个常见误区不要贴着有用信号的下边界设。比如有用信号最低到 5 Hz你把高通截止设成 5 Hz过渡带会把 5 到 10 Hz 的成分也削掉一部分。稳妥的做法是往下留一半截止设在 2 到 3 Hz。低频段的相对误差本来就大留余量比精确卡边重要得多。3.3 带通和带阻两个截止频率的顺序与写法差异带通和带阻都需要传两个截止频率区别在于btype和数组顺序。butter对Wn数组的要求是从小到大排列无论带通还是带阻都一样。带通的语义是“保留 low 到 high 之间”带阻是“挖掉 low 到 high 之间”。写反顺序会直接抛异常这个还算好排查真正麻烦的是把带通的语义理解反了结果滤出来的信号完全对不上。f1, f2 30.0, 80.0 wn_bp [f1 / (fs / 2), f2 / (fs / 2)] b_bp, a_bp signal.butter(4, wn_bp, btypebandpass) # 保留 30~80 Hz b_bs, a_bs signal.butter(4, wn_bp, btypebandstop) # 挖掉 30~80 Hz带阻滤波器最典型的应用就是工频干扰抑制。假设采样率 1000 Hz要挖掉 50 Hz 附近的干扰把带阻设成 48 到 52 Hz 就够了吗不一定。实际干扰往往有一定带宽而且电网频率本身会有微小漂移。我的经验是带宽设到 4 到 6 Hz也就是 47 到 53 Hz 左右同时阶数不要太高2 到 4 阶即可避免阻带边缘产生明显的振铃。如果干扰极其稳定、带宽很窄也可以考虑更窄的陷波器signal.iirnotch它专门针对单一频率点设计参数更少、更聚焦。提示带通和带阻的过渡带宽度直接受到两个截止频率间距的限制。间距很窄时IIR 滤波器的阶数需求会急剧上升数值稳定性也会变差。遇到这种情况优先考虑频域掩码法或者改用 FIR。3.4 零相位滤波 filtfilt好用的前提是你理解它的代价设计出b, a之后如果直接用signal.lfilter做因果滤波会引入相位延迟——波形整体向后平移而且不同频率的延迟量不同。对于“我只是想看波形”或者“特征提取”的场景这个延迟会污染结果。解决办法是signal.filtfilt它做前向加反向两次滤波相位正好抵消实现零相位。y_causal signal.lfilter(b_lp, a_lp, x) # 有相位延迟 y_zero signal.filtfilt(b_lp, a_lp, x) # 零相位幅频响应是平方代价有两个。第一幅频响应变成了原来的平方也就是说等效阶数翻倍、过渡带更陡、阻带衰减加倍。这通常是好事但意味着你按 4 阶设计出来的滤波器实际效果接近 8 阶如果不清楚这一点可能会误判参数。第二filtfilt在数据两端会有边界效应它通过内部填充来缓解但对很短的数据或者两端有突变的信号边缘仍然会失真。我的处理习惯是对数据两端各掐掉半个滤波器的瞬态长度再分析或者在做特征提取时避开首尾各 5% 的数据。还有一个容易忽略的点filtfilt要求输入长度至少是滤波器阶数的 3 倍左右太短的数据会报错。批量处理大量短片段时这个限制会突然冒出来提前加个长度判断能省很多调试时间。4. 频域滤波实操掩码、对称、IFFT 一条龙4.1 频域掩码的三种写法及适用边界用 FFT 做滤波核心就一步把频谱乘上一个 0/1 掩码或者平滑的权重函数。掩码的写法直接决定了结果对不对。第一种是理想矩形掩码通带内全 1、阻带内全 0对应理想滤波器过渡带为零时域会出现明显的吉布斯振铃。第二种是平滑过渡掩码比如升余弦窗牺牲一点过渡带换取时域波形干净。第三种是直接对幅值做缩放保留相位不动适合做增益调整而非硬性切除。N len(x) X np.fft.rfft(x) f np.fft.rfftfreq(N, d1/fs) # 写法一理想低通掩码 mask_ideal (f 100.0).astype(float) # 写法二升余弦过渡带90~110 Hz 之间平滑过渡 fc, tw 100.0, 10.0 mask_smooth np.ones_like(f) upper f (fc tw/2) lower f (fc - tw/2) mask_smooth[upper] 0.0 trans (~upper) (~lower) mask_smooth[trans] 0.5 * (1 np.cos(np.pi * (f[trans] - (fc - tw/2)) / tw)) y_ideal np.fft.irfft(X * mask_ideal, nN) y_smooth np.fft.irfft(X * mask_smooth, nN)理想掩码的振铃在什么场景下会出问题当信号里有明显的阶跃或脉冲成分时振铃会被放大成肉眼可见的波动看起来像是“滤出了假信号”。这时候换平滑掩码或者干脆回到时域用巴特沃斯效果都好得多。我个人的判断标准是做频谱分析、看能量分布用理想掩码没问题要把滤波后的波形拿去给人看或者做后续检测优先平滑掩码或时域 IIR。4.2 实信号频谱对称性与 irfft 的实数化处理前面提过实信号的频谱共轭对称所以用rfft只存一半是最省事的。但如果你出于某些原因用了完整fft掩码就必须写成np.abs(f) fc这种对称形式。只对正频率做掩码ifft的结果会是复数取.real之后幅值减半而且会残留一个随频率变化的相位误差波形看起来“像是对的但幅度不对”极其误导。# 错误示范只对正频率保留 X_full np.fft.fft(x) f_full np.fft.fftfreq(N, 1/fs) mask_wrong (f_full 0) (f_full 100) y_wrong np.fft.ifft(X_full * mask_wrong).real # 幅值约为正确值的一半 # 正确示范对称保留 mask_right np.abs(f_full) 100 y_right np.fft.ifft(X_full * mask_right).real即使掩码正确ifft的结果类型仍然是复数虚部是浮点误差残留量级通常在 1e-16 左右。直接.real取实部即可不需要额外处理。但如果虚部量级明显偏大比如 1e-3那说明掩码对称性或者长度处理出了问题得回头查不能简单取实部糊弄过去。4.3 一段含噪信号的完整处理流程把上面的东西串起来走一遍完整流程造一段包含低频有用信号、高频噪声和单频干扰的合成数据分别用频域掩码和时域 IIR 两条路处理对比结果。这个流程我建议每个刚上手的人都自己跑一遍比看十篇文章都管用。import numpy as np from scipy import signal import matplotlib.pyplot as plt fs, T 1000.0, 2.0 t np.arange(0, T, 1/fs) N len(t) # 合成信号5 Hz 有用成分 50 Hz 干扰 宽带噪声 useful 1.0 * np.sin(2*np.pi*5*t) interf 0.5 * np.sin(2*np.pi*50*t) noise 0.3 * np.random.randn(N) x useful interf noise # 路线一时域带阻去 50 Hz 低通去高频噪声 b_bs, a_bs signal.butter(4, [47/(fs/2), 53/(fs/2)], btypebandstop) b_lp, a_lp signal.butter(4, 30/(fs/2), btypelow) y_time signal.filtfilt(b_lp, a_lp, signal.filtfilt(b_bs, a_bs, x)) # 路线二频域掩码直接保留 0~30 Hz X np.fft.rfft(x) f np.fft.rfftfreq(N, 1/fs) mask (f 30.0).astype(float) y_freq np.fft.irfft(X * mask, nN) # 对比与干净有用信号的相关系数 corr_time np.corrcoef(y_time, useful)[0, 1] corr_freq np.corrcoef(y_freq, useful)[0, 1] print(f时域方案相关系数: {corr_time:.4f}) print(f频域方案相关系数: {corr_freq:.4f})实测下来两种方案的相关系数通常在 0.99 以上差别主要在中段和边缘时域方案边缘更干净但有轻微相位残留频域方案整体幅值更准但两端会有小幅波动环形卷积效应。频域滤波本质上是循环卷积它默认信号首尾相接所以当首尾幅值差异大时接缝处会出现跳变。解决办法是滤波前做端点延拓或者对结果做首尾裁剪。这个细节在任何教程里都很少提但实际做长序列处理时一定会撞上。5. 常见问题与排查速查表5.1 现象、原因、处理对照现象最可能的原因处理方式butter报错 Wn 超出范围把赫兹直接传进去了除以fs/2做归一化滤波后整体幅值减半掩码只保留了正频率改用np.abs(f)对称掩码或rfft波形整体向后平移用了lfilter而非filtfilt换filtfilt或补偿群延迟滤波后信号变“毛刺”理想掩码导致振铃改用平滑过渡掩码或降阶频谱峰位置和预期不符频率轴没用fftfreq生成用fftfreq(N, 1/fs)重建边缘出现异常波动循环卷积接缝效应端点延拓后滤波或裁剪首尾补零后频谱仍分不开两个峰补零不增加真实分辨率延长采集时间或降低采样率高通后基线反而更抖截止频率太靠近信号低频端截止频率下移留出过渡带这张表里的每一条我都至少踩过一次。最想强调的是第二条——幅值减半。这个错误太隐蔽了因为波形形状、频率成分看起来完全正确只是整体小了一半如果你没有和原始信号做幅值对比根本发现不了。养成习惯滤波后立刻算一下通带内某个已知频率成分的幅值和滤波前对比一秒钟就能验证。5.2 几个我实际踩过的坑第一个坑是采样率和时间轴不一致。有次数据文件里记录的时间戳和实际采样点数对不上我按点数算出来的fs偏差了 2%导致所有截止频率都偏了 2%。低频段 2% 影响不大但带阻位置偏了之后工频干扰根本没压住我一度以为是滤波器阶数不够白折腾了半天。后来加了校验abs(len(t) / (t[-1] - t[0]) - fs) 0.01不一致就报警。第二个坑是对非平稳信号用固定截止频率。一段信号前半段噪声强、后半段噪声弱用同一个滤波器处理后半段会被过度平滑丢掉真实细节。这种场景更适合自适应方法或者分段处理每段单独估计噪声水平再定参数。固定参数的滤波器适合“干扰特性稳定”的前提这一点在动手前要想清楚。第三个坑是忽略数据类型。从文件里读进来的数据有时候是整数类型int16 之类直接做 FFT 或者滤波中间运算会悄悄截断结果出现台阶状的失真。凡是做数值处理第一步统一x np.asarray(x, dtypenp.float64)这一行能挡掉一大类莫名其妙的精度问题。最后分享一个我常用的验证套路滤波前后各画一张频谱图把它们叠在一起看。通带内的谱线应该几乎重合如果用了filtfilt幅频是平方关系会有轻微差异阻带内的谱线应该明显下压。这比看时域波形直观得多也更容易发现参数设错的问题。我现在的习惯是任何滤波脚本里都内置这个对比图跑一次就心里有数。