★ FEATURED ARTICLE 简介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 Hz20 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是这个算法的理论底座。它指出若一个等波纹滤波器是最优的其加权误差函数在频带上至少要有L2个交替点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, weightweight, fsfs)这段代码干了什么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参数只能传nyqfs/2bands要写成归一化频率比如上面的例子要写成[0, 1/6, 1/4, 1]手算归一化非常容易除错。显式传fs48000frequency值就是Hz不仅直观还能少一类低级错误。另一个容易忽视的点numtaps给太小算法会返回一个能算但性能很差的结果甚至直接不收敛。给太大算起来慢实时系统里还浪费乘法器。后面第4章专门讲阶数估计。第一次验证时不纠结最优阶数先给一个足够大的奇数比如201把流程跑通再说。3.3 用freqz验证纹波、衰减一次看明白import matplotlib.pyplot as plt w, h freqz(taps, worN4096, fsfs) 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(fPassband ripple {ripple:.3f} dB) plt.show()这里的细节是freqz传了fs后返回的w单位是Hz可以直接和截止频率对照。通带纹波别拿人眼估用np.ptp计算通带区间幅度最大值减最小值这是定量验收的习惯。阻带衰减用-np.max(mag_db[stopband])我见过太多人只看图觉得够低了实际一量差了3 dB。还有worN是频响计算点数最少给4096纹波很密的滤波器建议8192。点数给小了峰值会被漏采样评估结果比实物乐观设计完上板实测不达标才尴尬。我习惯先worN8192跑一遍确认峰值都露出来了再换4096做日常迭代省时间。4. 实战带通滤波器从阶数估算到纹波达标4.1 带通设计完整脚本带通场景采样率48 kHz通带3~7 kHz要求2 kHz以下衰减大于60 dB9 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, weightweight, fsfs) w, h freqz(taps, worN8192, fsfs) mag_db 20 * np.log10(np.abs(h) 1e-12) # 通带纹波 pb (w 3000) (w 7000) print(fPassband ripple: {np.ptp(mag_db[pb]):.3f} dB) # 阻带衰减 sb1 w 2000 sb2 w 8000 print(fStopband att: {-np.max(mag_db[sb1]):.1f} dB / {-np.max(mag_db[sb2]):.1f} dB)注意bands里频段是交替出现的被约束的频段之间自然夹着过渡带过渡带不参与误差逼近。这里的两个过渡带都是1 kHzdesired对应三段阻带目标0、通带1、另一个阻带0weight也同步三段。过渡带越宽相同阶数下纹波越好但过渡带宽会把有用边缘信号削掉设计时别为了指标好看而故意放宽要按真实系统需求定边界。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过渡带1kHzfs48k 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 dB2 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给太小阻带实际最高点刚好落在采样点之间你测到的是低谷误以为衰减达标一上板就露馅。解决先用worN8192重建频响确认阻带最大点被捕捉到。再用np.max对阻带区间做定量计算别只看图。若数值确实没提升说明阶数才是瓶颈放弃调权重回第4章加阶数。5.3 现象低通用偶数抽头没问题带通和高通阻带末端忽然抬起来原因是偶数抽头对应II型FIR滤波器频率响应在奈奎斯特频率处被强制为0带通或高通设计里阻带边缘靠近奈奎斯特的地方会因为这个结构性零点出现塌陷或回升导致阻带形态失控。低通因为通带在低频端影响不明显所以很多人没意识到。解决带通、高通设计一律用奇数抽头也就是I型滤波器。如果你确实需要偶数抽头比如要凑DSP的偶对齐就得接受奈奎斯特附近的频响约束通常需要在最上边留出一条过渡带别把约束频段直接顶到fs/2。5.4 现象归一化频率算错频响图整体变形原因是老代码里用nyq参数时bands要除以奈奎斯特频率。举例fs48000时边界4000 Hz要写成4000/240001/66000 Hz写成1/4手算容易把分母写成采样率结果过渡带、通带全部错位。解决显式传fs48000bands直接写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, weightweight, fsfs)这里bands把25~40 kHz划为高权重阻带45~48 kHz划为低权重阻带中间40~45 kHz留作过渡带。高权重段会沉得很深低权重段相对浅形成一个缓变台阶。这种同是阻带但深浅不同的设计才是Remez在工程里最值钱的地方它让滤波器的频响形状真正跟着系统需求走而不是被一种固定的窗函数形状绑架。验证技巧上除了freqz看频响我还会做一个信号级测试生成一段包含通带、阻带多个频率分量的叠加信号让滤波器跑一遍用FFT看输出频谱各频点幅度。这个测试和频响曲线对得上才敢把系数丢进定点化流程。定点化之后也别忘了复测一次系数从浮点变成16位整数阻带衰减可能会掉3~5 dB设计时要把这部分余量提前留出来。说到底Remez交换这条路线干净、可控、边界清晰能处理的频带形状远超窗函数法。我自己的习惯是先跑通最小例程再按阶数估算起算、迭代回缩最后用信号级验证收尾。这一套下来滤波器设计基本不再返工。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
