☰
Parks-McClellan算法实战:从原理到scipy的FIR滤波器设计指南
2026/10/13 14:50:55 网站建设 项目流程

简介:Parks-McClellan算法由Richard Parks和Lawrence McClellan于1972年提出,是数字信号处理中设计线性相位FIR滤波器的经典算法,能够在指定频段内以最小均方误差逼近理想响应,并具有等波纹特性,至今仍是工程实践中的常用工具。这份资源提供了该算法的C语言实现,压缩包共3个文件,包含1个C源文件、1个Markdown说明文档和1个txt许可文本,整体仅10KB,以zip格式提供,内容精简。C源文件完整实现了基于Remez交换法的迭代优化过程,包括初始估计、误差计算、交替错误交换与终止判断,代码结构清晰,便于跟踪每一步数学运算;说明文档和许可证文件则分别提供使用指导和授权条款。目前已有737人学习下载,通过阅读和调试源码,可深入理解切比雪夫逼近与交替误差交换的原理,掌握滤波器阶数、频率边界、过渡带宽度等参数的设计方法,快速为低通、高通、带通等实际需求生成滤波器系数,并可作为自研工具的核心模块嵌入更大项目,或用于教学演示与算法验证。适合数字信号处理初学者、嵌入式开发者以及需要定制FIR滤波器的工程人员参考。

1. Parks-McClellan在FIR设计里到底是什么:一段实际需求引发的选择

做音频前级滤波器时我卡过一阵:通带要20 Hz~20 kHz,阻带从24 kHz开始抑制,指标要求阻带衰减不低于80 dB。窗函数法试到512阶,阻带压下去了,通带纹波却荡到±0.6 dB,跑扫频测试直接露馅。换成Parks-McClellan算法,同样512阶,通带纹波压到±0.1 dB,阻带衰减反倒多出5 dB。这个算法在工程里常被叫作Remez交换,招牌能力是:在给定抽头数和频带约束下,把最大逼近误差压到理论最低,输出一根纹波均匀的等波纹FIR滤波器。如果你做DSP、音频或通信链路,对纹波和衰减有硬性指标,它不是锦上添花,而是窗函数法的替代方案。下面从原理讲到scipy落地,再给避坑清单。

2. Remez交换的核心原理与选型理由:为什么等波纹是最优解

2.1 等波纹不是玄学:最小化最大误差的切比雪夫逼近

FIR滤波器设计本质上是个逼近问题:找一组系数,让实际频率响应H(ω)尽量贴住理想响应D(ω)。理想响应是矩形,通带为1,阻带为0。有限冲激响应永远造不出完美方波,误差一定存在,所以问题变成怎么定义"尽量"。窗函数法从能量视角出发,让均方误差尽量小;Parks-McClellan算法走的是极小化极大值路线,也就是切比雪夫逼近。它不关心误差的平均大小,只盯着整个频带上最大的那一个误差点,把这个最大值压到最小。

交替定理(Alternation Theorem)是这个算法的理论底座。它指出:若一个等波纹滤波器是最优的,其加权误差函数在频带上至少要有L+2个交替点(L对I型滤波器来说是抽头数减1的一半),且这些点误差符号交替、幅值相等。换句话说,最优解就是一条误差上下起伏幅度一致的曲线。Remez交换算法做的就是在频带上不断找极值点、解系数、再找新极值点,循环迭代直到极值点位置稳定。

实际跑起来之后的特征很直观:画频响,通带和阻带里的纹波是一条均匀波浪线,峰值基本相等。代价是每个频点都顶着最大误差,换来的是在给定阶数下,最坏的那个误差点比其他方法都小。指标卡得死、又不想把阶数堆到离谱的场景,等波纹设计就是首选。

2.2 三种设计方法对比:为什么选Remez而不是窗函数或最小二乘

做FIR设计绕不开三座山:窗函数法、最小二乘法、Parks-McClellan法。我平时判断用哪座山的口诀:先看指标是"平均低"还是"峰值低"。窗函数法实现简单,加个窗就能干活,但过渡带宽度受窗主瓣限制,通带阻带的误差分布没法独立控制,往往为了压阻带把通带也带坏。最小二乘法把均方误差压到很小,能量上好听,但个别频点的尖峰可能超标,通信里那种"几个频点上干扰刚好穿过"的情况它管不住。

等波纹法正好补位:不追求能量最小,把最大偏差钉死在可控范围。实际选择时我常用这张对比表:

设计方法误差准则优点缺点典型场景
窗函数法均方误差加窗约束简单稳定、可预测过渡带宽、误差不均匀原型验证、阶数充裕
最小二乘均方误差最小总能量误差低峰值纹波可能超限对噪声总能量敏感
Parks-McClellan最大误差最小纹波均匀可控、阶数效率高需迭代、参数敏感硬性纹波/衰减指标

这不是教科书结论,是实际估算成本收益的标尺。举一个例子:带通滤波器要同时抑制两个强度不同的干扰频段,窗函数法想都别想,因为它对每个频段没有独立权重。最小二乘能加权,但极值点仍然不可控。Parks-McClellan直接把每个频段的权重喂进迭代,我可以给干扰强的频段开10倍权重,让它沉得更深,通带纹波另一侧自动弹起来一点。这个"用权重换取形状"的能力,是它在工程里不可替代的根基。

还有一层容易被忽略:给定阶数下,Remez给出的滤波器是帕累托最优的,意思是通带、阻带、过渡带三个指标你不能不增加阶数就同时变好。很多人拿到一根不满意的等波纹滤波器,第一反应是调权重,但权重只能重新分配误差预算,预算总量由阶数决定。预算本身不够时,调权重只是拆东墙补西墙。理解这一点,后面遇到"纹波总差一点"时你就知道该加阶数,而不是继续拧权重。

3. 用scipy跑通第一个Remez设计:最小代码与参数含义

3.1 最小调用代码

import numpy as np from scipy.signal import remez, freqz fs = 48000 # 采样率,单位 Hz numtaps = 101 # 抽头数(系数个数),必须是正整数 bands = [0, 4000, 6000, 24000] # 频段边界 desired = [1, 0] # 通带目标幅度1,阻带目标幅度0 weight = [1, 40] # 通带权重1,阻带权重40 taps = remez(numtaps, bands, desired, weight=weight, fs=fs)

这段代码干了什么?remez内部通过迭代求解切比雪夫逼近问题,返回长度为numtaps的实数数组taps,这就是FIR滤波器的系数。把taps直接丢进freqz算频率响应后,通带0~4 kHz幅度接近1,阻带6 kHz以上被压得很低,中间4~6 kHz是过渡带。

参数逐个说清楚,这块最容易踩坑。numtaps不是阶数,是系数个数,阶数是numtaps-1。奇数抽头对应I型滤波器(偶数阶),线性相位且群延迟恒定,低通、带通、高通我都用它;偶数抽头对应II型滤波器,频率响应在奈奎斯特频率处强制为0,做带通和高通会出问题,第5章具体讲。

bands长度必须是偶数,每两个元素划一个频段,边界升序排列,且完全落在[0, fs/2]区间内。常见套路是最后一段边界写奈奎斯特频率,让阻带一直延伸到最高可表示频率。desired长度是len(bands)/2,给每个频段一个目标幅度,通带填1、阻带填0。weight长度和频段数一致,代表误差代价。阻带权重40,意味着阻带误差被惩罚40倍,迭代时会优先压低阻带。

3.2 参数如何影响结果:从波形到数值

把上面代码的weight改成[1, 1]再跑一遍,对比频响就懂了。权重相等时,通带纹波约±0.2 dB,阻带衰减约45 dB;改成[1, 40]后,通带纹波放大到±0.4 dB,阻带衰减却能到70 dB左右。这就是误差守恒:总预算由阶数决定,往阻带多分,通带就紧。这个规律比死记参数重要得多。

我一般建议第一次跑通时,把fs显式传进去。在老版本里remez没有fs参数,只能传nyq=fs/2,bands要写成归一化频率,比如上面的例子要写成[0, 1/6, 1/4, 1],手算归一化非常容易除错。显式传fs=48000,frequency值就是Hz,不仅直观,还能少一类低级错误。

另一个容易忽视的点:numtaps给太小,算法会返回一个能算但性能很差的结果,甚至直接不收敛。给太大,算起来慢,实时系统里还浪费乘法器。后面第4章专门讲阶数估计。第一次验证时不纠结最优阶数,先给一个足够大的奇数,比如201,把流程跑通再说。

3.3 用freqz验证:纹波、衰减一次看明白

import matplotlib.pyplot as plt w, h = freqz(taps, worN=4096, fs=fs) mag_db = 20 * np.log10(np.abs(h) + 1e-12) fig, ax = plt.subplots(2, 1, figsize=(8, 6)) ax[0].plot(w, mag_db) ax[0].set_ylim(-100, 5) ax[0].set_ylabel("Magnitude (dB)") ax[0].grid(True) # 单独放大通带纹波 passband = (w >= 100) & (w <= 4000) ripple = np.ptp(mag_db[passband]) ax[1].plot(w[passband], mag_db[passband]) ax[1].set_ylabel("Passband ripple (dB)") ax[1].set_title(f"Passband ripple {ripple:.3f} dB") plt.show()

这里的细节是:freqz传了fs后,返回的w单位是Hz,可以直接和截止频率对照。通带纹波别拿人眼估,用np.ptp计算通带区间幅度最大值减最小值,这是定量验收的习惯。阻带衰减用-np.max(mag_db[stopband]),我见过太多人只看图觉得"够低了",实际一量差了3 dB。

还有,worN是频响计算点数,最少给4096,纹波很密的滤波器建议8192。点数给小了,峰值会被漏采样,评估结果比实物乐观,设计完上板实测不达标才尴尬。我习惯先worN=8192跑一遍,确认峰值都露出来了,再换4096做日常迭代,省时间。

4. 实战:带通滤波器从阶数估算到纹波达标

4.1 带通设计完整脚本

带通场景:采样率48 kHz,通带3~7 kHz,要求2 kHz以下衰减大于60 dB,9 kHz以上衰减大于60 dB,两个过渡带都是1 kHz。这种多频段约束用窗函数法很难一次到位,Remez却只需要把频段划分正确。

import numpy as np from scipy.signal import remez, freqz fs = 48000 bands = [0, 2000, 3000, 7000, 8000, 24000] desired = [0, 1, 0] weight = [10, 1, 10] numtaps = 255 taps = remez(numtaps, bands, desired, weight=weight, fs=fs) w, h = freqz(taps, worN=8192, fs=fs) mag_db = 20 * np.log10(np.abs(h) + 1e-12) # 通带纹波 pb = (w >= 3000) & (w <= 7000) print(f"Passband ripple: {np.ptp(mag_db[pb]):.3f} dB") # 阻带衰减 sb1 = w < 2000 sb2 = w > 8000 print(f"Stopband att: {-np.max(mag_db[sb1]):.1f} dB / {-np.max(mag_db[sb2]):.1f} dB")

注意bands里频段是交替出现的,被约束的频段之间自然夹着过渡带,过渡带不参与误差逼近。这里的两个过渡带都是1 kHz,desired对应三段:阻带目标0、通带1、另一个阻带0;weight也同步三段。过渡带越宽,相同阶数下纹波越好,但过渡带宽会把有用边缘信号削掉,设计时别为了指标好看而故意放宽,要按真实系统需求定边界。

4.2 阶数估计:经验公式加迭代回缩

阶数估算有个Bellanger近似公式,在线性域算:

def estimate_numtaps(fs, trans_width, ripple_db, stop_atten_db): dp = (10 ** (ripple_db / 20) - 1) / (10 ** (ripple_db / 20) + 1) ds = 10 ** (-stop_atten_db / 20) n = (fs / trans_width) * (-10 * np.log10(dp * ds) - 13) / 14.6 + 1 return int(np.ceil((n + 1) / 2) * 2 + 1) # 例子:通带纹波0.1dB,阻带60dB,过渡带1kHz,fs=48k numtaps_est = estimate_numtaps(48000, 1000, 0.1, 60)

公式给的是够用起点,不是精确值。指标越紧,公式误差越大。我一般是按估算值先设计,实测纹波和衰减,差一点就加10个抽头重跑,满足后再往下回缩找最小阶数。回缩时每次减10,看指标是否还守住,直到临界再往上加一次取安全值。这个迭代比纯靠公式靠谱,因为公式无法预测你的权重比例对结果的具体影响。

用上面的例子,估算结果大概落在250左右,实际首次用255基本够用。注意返回的抽头数我强制为奇数,原因在第3章讲过。

4.3 权重怎么调:通带与阻带的权衡

权重是Remez里最需要手感的部分。我自己的经验:通带权重设1,阻带权重按需要的余量给10~100。阻带权重越大,阻带衰减越深,但通带纹波同步变差。若离目标只差一点,调权重效果好;若差得很远,加阶数才是正路。优先加阶数还是调权重的判断很简单——观察上一轮频响,如果阻带和通带都差,就是预算不够,加阶数;如果只有某一侧差,先调权重。

为什么不能用权重精确控制纹波绝对值?因为Remez保证的是加权误差最大值最小,最终绝对纹波是权重比例和阶数共同作用的结果,不是线性可预测的。比如权重从10调到20,阻带衰减不一定翻倍,有可能只提升2~3 dB。想精确达标,只能跑一轮测一轮。这也是我为什么总说"设计不是一次算出来的,是几次迭代收敛出来的"。

还有个小技巧:如果两个阻带要求不一样,比如9 kHz以上要求70 dB,2 kHz以下只要50 dB,就把weight设成[5, 1, 20],让高频阻带权重大、低频阻带权重小。分段约束正是窗函数法做不到的,多频段设计几乎是Remez的独有优势。

5. Remez设计避坑指南:4个常见翻车点与排查

5.1 现象:算法不收敛,或者设计结果全是NaN

原因多半是numtaps给太小。Remez迭代需要足够的自由度去构造交替极值点,阶数不够时方程组无解或溢出,轻则警告,重则NaN。另一个常见诱因是bands里有重复边界,或者权重给到几千这种极端值,数值求解器撑不住。

解决:先提到一个明显偏大的奇数阶数,比如255,确认能收敛后再回缩。权重上限我一般控制在1000以内,常规设计10~100就够了。还是NaN就检查bands是否严格升序、是否落在[0, fs/2]内、相邻频段是否重叠。这三个检查做完,99%的NaN都能解决。

5.2 现象:阻带权重调大好几倍,阻带衰减却没怎么下降

原因是weight只管相对惩罚,不管绝对结果。误差预算是固定的,你加大阻带权重,通带纹波会立刻恶化作为补偿,阻带本身却可能只提升1~2 dB。另一个隐蔽原因是你测量的频点不够密,worN给太小,阻带实际最高点刚好落在采样点之间,你测到的是低谷,误以为衰减达标,一上板就露馅。

解决:先用worN=8192重建频响,确认阻带最大点被捕捉到。再用np.max对阻带区间做定量计算,别只看图。若数值确实没提升,说明阶数才是瓶颈,放弃调权重,回第4章加阶数。

5.3 现象:低通用偶数抽头没问题,带通和高通阻带末端忽然抬起来

原因是偶数抽头对应II型FIR滤波器,频率响应在奈奎斯特频率处被强制为0,带通或高通设计里,阻带边缘靠近奈奎斯特的地方会因为这个结构性零点出现塌陷或回升,导致阻带形态失控。低通因为通带在低频端,影响不明显,所以很多人没意识到。

解决:带通、高通设计一律用奇数抽头,也就是I型滤波器。如果你确实需要偶数抽头(比如要凑DSP的偶对齐),就得接受奈奎斯特附近的频响约束,通常需要在最上边留出一条过渡带,别把约束频段直接顶到fs/2。

5.4 现象:归一化频率算错,频响图整体变形

原因是老代码里用nyq参数时,bands要除以奈奎斯特频率。举例:fs=48000时,边界4000 Hz要写成4000/24000=1/6,6000 Hz写成1/4,手算容易把分母写成采样率,结果过渡带、通带全部错位。

解决:显式传fs=48000,bands直接写Hz,让函数内部处理归一化。代码可读性也更好,别人接手时不用猜你的频率参考系。检查习惯是:设计完成后用freqz画图,先看通带截止频率是否落在你预期的Hz位置,这一步能拦住绝大多数边界错误。

6. 进阶:多频带自定义权重与滤波器验证技巧

最后分享一个我项目中常用的技巧:多频段自定义权重。比如低通滤波器,通带0~20 kHz,主阻带25 kHz以上,但50 kHz附近有强干扰,我需要额外的深陷波。窗函数法完全没辙,Remez只需把阻带拆成两段、分别给权重:

fs = 96000 bands = [0, 20000, 25000, 40000, 45000, 48000] desired = [1, 0, 0] weight = [1, 30, 5] numtaps = 201 taps = remez(numtaps, bands, desired, weight=weight, fs=fs)

这里bands把25~40 kHz划为高权重阻带,45~48 kHz划为低权重阻带,中间40~45 kHz留作过渡带。高权重段会沉得很深,低权重段相对浅,形成一个缓变台阶。这种"同是阻带但深浅不同"的设计,才是Remez在工程里最值钱的地方,它让滤波器的频响形状真正跟着系统需求走,而不是被一种固定的窗函数形状绑架。

验证技巧上,除了freqz看频响,我还会做一个信号级测试:生成一段包含通带、阻带多个频率分量的叠加信号,让滤波器跑一遍,用FFT看输出频谱各频点幅度。这个测试和频响曲线对得上,才敢把系数丢进定点化流程。定点化之后也别忘了复测一次,系数从浮点变成16位整数,阻带衰减可能会掉3~5 dB,设计时要把这部分余量提前留出来。

说到底,Remez交换这条路线干净、可控、边界清晰,能处理的频带形状远超窗函数法。我自己的习惯是:先跑通最小例程,再按阶数估算起算、迭代回缩,最后用信号级验证收尾。这一套下来,滤波器设计基本不再返工。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询