简介这是一份用 C 语言实现的 Parks-McClellan 算法remez-exchange源码面向数字信号处理学习者、嵌入式开发者以及需要设计线性相位 FIR 滤波器的工程师。算法采用切比雪夫逼近和交替误差交换策略通过迭代优化滤波器系数在指定频段内获得等波纹equiripple频率响应适合教学实验、滤波器原型验证、通信系统仿真以及工程移植等场景。压缩包共 3 个文件包含 C 语言源文件、Markdown 说明文档和许可证文本整体约 10KB结构精简便于快速阅读和直接编译。已有 738 人下载学习说明该源码在 FIR 滤波器设计场景中具有一定参考价值。使用者可直接运行源码生成滤波器系数也可结合说明文档理解迭代优化的完整步骤并根据滤波器阶数、通带/阻带边界和权重等需求修改参数从而深入掌握 Parks-McClellan 算法的原理及 C 语言实现细节。1. Parks-McClellan算法为什么“误差全部相等”才是FIR设计的最优解理想低通滤波器在物理上不可实现——频域是矩形时域就是无限长的sinc。工程上所有FIR设计方法本质都在回答同一个问题用一个有限长冲激响应去逼近那个不可能得到的理想响应逼近得“最好”。窗函数法靠截断频率采样法靠插值而Parks-McClellan算法给出的答案完全不同它把设计问题转化成切比雪夫逼近用Remez交换迭代去解最后得到的滤波器误差曲线是一条上下均匀起伏的“等波动”带通带和阻带里每个峰值误差都恰好一样大。这正是这个算法最反直觉的地方最优解不是让误差“平均最小”而是让“最大误差”最小并且所有最大误差点严格相等、符号交替。要理解为什么这个条件就是最优以及一段用C写出来的remez程序凭什么能迭代收敛需要从头拆开交替定理和Remez交换的机制。2. 从窗函数到切比雪夫逼近要“平均误差小”还是“最大误差小”2.1 窗函数法输在哪旁瓣不是你能控制的窗函数法的思路是把理想冲激响应 h_d[n] 直接乘以一个有限长窗口 w[n]得到 h[n] h_d[n]·w[n]。时域相乘等于频域卷积理想矩形频谱和窗口频谱卷积的结果决定了通带边缘的过渡带宽度和阻带衰减。问题在于一旦你选了窗口类型过渡带宽度和阻带衰减的配对关系就被锁死了。矩形窗阻带衰减只有 -21dB汉宁窗是 -44dB海明窗是 -53dB布莱克曼窗是 -74dB。你没法说“我要 -60dB 衰减、过渡带 0.01π”——窗函数只给你固定配方的折中指标是切不出来的。更麻烦的是通带和阻带共享同一个卷积核通带波纹和阻带波纹不可分别调节这对很多实际系统是不够用的。换窗能改善旁瓣但改善的机制是把能量从旁瓣推远代价是主瓣变宽、过渡带变缓。这是时域截断这个动作本身带来的物理限制不是实现精度问题。所以窗函数法适合快速出原型不适合做有严格指标的滤波器。2.2 频率采样法输在哪只在采样点上精确频率采样法把理想响应在N个等间隔频率点上采样然后 IDFT 得到系数。它的优点是设计过程零迭代、极快缺点是只在采样点上严格等于目标值两个采样点之间的响应完全不受控。当理想响应在过渡带边缘有陡峭跳变时点间会产生很大的振荡误差峰值甚至比窗函数法还大。可以这么理解频率采样法是一个在N个点上精确、在其他频率放任自流的插值逼近。它没有把“整个频带上的误差”作为一个整体来优化所以谈不上最优。2.3 极小极大准则从“平均误差”换到“最大误差”既然前两种方法都缺乏全局误差控制正确的数学化方式是什么逼近论里有一个现成的框架设理想响应为 D(ω)设计出的滤波器频率响应为 A(ω)误差为 E(ω) W(ω)·(A(ω) − D(ω))其中 W(ω) 是用户给的权重函数——你对哪个频带更敏感就给哪个频带更大的权重。切比雪夫逼近也叫极小极大逼近的目标是minimize max |E(ω)|, ω ∈ 通带 ∪ 阻带也就是说最小化整个频带上的最大加权误差。注意这里不含过渡带——过渡带没有目标值误差在那里不参与优化这正是它存在的意义过渡带越窄代价是通带和阻带的波纹越大反之亦然。用户用过渡带宽度换取带内误差这是其他方法给不了的旋钮。而“最大误差最小”这件事远比听起来苛刻。它要求最优解在所有局部极值点上的误差幅度完全一样——否则某个大误差点可以被压低直到它和其他极值点一样高。最终误差曲线必然是等波动形状这正是Parks-McClellan算法产物的视觉签名。3. 交替定理最佳逼近的“指纹”以及Remez交换为什么能收敛3.1 交替定理究竟说了什么对于抽头数为N的线性相位FIRN为奇数Type I对称其幅度响应 A(ω) 可以写成 cos(ω) 的 M 次三角多项式其中 M (N−1)/2。于是设计问题归结为找一个 M 次多项式去加权逼近给定的目标函数 D(ω)使最大加权误差最小。这一节的数学结论是交替定理若加权误差 E(ω) 在频带上有至少 M2 个极值点并且这些极值点处误差的符号正负交替、幅度都等于同一个值 δ那么这个多项式就是唯一的最佳逼近δ 就是能达到的最小最大误差。这个“M2个交替符号的极值点”就是最优解的指纹。直观理解一个有 M 个自由参数的逼近误差曲线如果只有 M1 个极值你还能调整参数把最高的峰压下去、导致另一个峰升高直到峰与峰拉平。只有交替峰数达到 M2系统才“卡死”——任何微调都会让某个已经等于理论下限的峰更高最优性被锁定了。注意极值点必须同时覆盖通带和阻带两个频带的边界过渡带边缘天然是候选极值点。3.2 Remez交换的四个步骤网格、插值、误差、换点交替定理给了判据但怎么找到这组极值点Remez交换算法是一个不断“猜测→验证→修正”的迭代过程初始化在频率网格上等间隔取 M2 个点作为初始极值点集合插值假设这组极值点处误差都等于 ±δ符号交替用拉格朗日插值反解出当前的多项式系数和 δ 值求误差在非常密集的全网格上计算加权误差 E(ω)找出所有局部极值点交换从这些局部极值点中选出绝对值最大且符合符号交替要求的 M2 个点替换掉旧极值点重复 24直到连续两轮的极值点集合不再变化、δ 的波动小于阈值。为什么这个循环能收敛关键在于第 2 步解出的 δ 一定小于上一轮的真实最大误差——因为上一轮的极值点并非真正的等波动点插值得到的“假定δ”不是全局最大被低估了。第 3 步在全网格上扫描必然能找到超过这个 δ 的新极值点把它们换进来下一轮的 δ 就会增大。这个单调过程最终收敛到理论最小最大误差。实际程序里通常几轮迭代 δ 就稳定了Remez交换的收敛速度在工程上完全够用。4. 核心代码拆解C语言实现里的网格、插值与极值点搜索一个典型的P-M算法C实现主循环看起来并不长难在几个细节函数。下面这段是主循环的核心骨架注释标出了每步的机制。/* remez 主迭代 * L (N-1)/2, N为抽头数 * ext[] 存极值点在网格上的下标大小 L2 */ static int remez_iter(double grid[], double des[], double wt[], int gridsize, int L, int ext[], double *delta, int max_iter) { int iter, i, n_ext, converged 0; for (iter 0; iter max_iter; iter) { /* 1. 由当前极值点解切比雪夫方程组算出假定误差上界 delta * 经典实现不显式解矩阵改用重心拉格朗日插值一次求出 delta */ *delta calc_delta(ext, L, grid, des, wt); /* 2. 在全网格上计算误差函数 err[i] wt[i] * (H(i) - des[i]) */ calc_error_curve(grid, des, wt, err, gridsize, ext, L); /* 3. 扫描全网格找符号交替、幅度超过 |delta| 的候选极值点 */ n_ext find_extrema(ext, grid, err, gridsize, L, *delta); /* 4. 极值点集合不再变化 - 收敛 */ if (n_ext L 2 same_extrema(ext, prev_ext)) { converged 1; break; } /* 5. 交换用新极值点集合覆盖旧集合进入下一轮 */ memcpy(prev_ext, ext, (L 2) * sizeof(int)); } return converged ? iter : -1; }calc_delta承担的是交替定理里“假定所有极值点误差相等”的角色它把当前 M2 个极值点看作已经处于等波动状态用插值倒推多项式在某个参考点处的值和 δ。这里用重心形式的拉格朗日插值barycentric interpolation而不是直接解矩阵方程因为后者在极值点位置靠近时矩阵条件数会爆炸而重心插值数值上稳得多。find_extrema则是交替定理的“验证器”它只接受符号正负交替的点连续两个同号极值点会被丢弃因为它们不可能构成最优解。网格密度和分段方式也是机制的一部分实现把 0 到 0.5归一化频率单位cycles/sample分成各频带每个频带内部再细分网格点默认每频带取约 16 倍的点数过渡带边缘点会强制加入极值点候选。网格太疏会漏掉真正的局部极值点导致迭代停滞在错误解太密则纯增计算量。这是这个算法里“工程直觉”最集中的地方——数学上要求扫描所有频率工程上只能在足够密的离散网格上近似。5. 避坑Remez迭代翻车的几个典型现场1. 迭代不收敛delta 上下抖动现象程序跑满 max_iter 次误差曲线仍然不是等波动delta 在相邻轮次之间反复跳变。原因最常见的是网格太稀疏局部极值点被漏掉导致 find_extrema 返回的候选点集合不稳定也可能是初始极值点个数不对比如抽头数给错导致 L 与网格大小不匹配。解决把每频带的网格密度从 16 倍翻到 32 倍再跑同时检查调用端传入的抽头数和频带数组是否一致尤其是奇数抽头数Type I不要误传成偶数。2. 阻带衰减始终达不到指标误差曲线“一边高一边低”现象设计结果通带波纹很小但阻带衰减差得远或者反过来阻带很好通带超差。原因加权函数 W(ω) 没有正确表达你的偏好。P-M算法最小化的是加权误差不是绝对误差——你不给阻带加权重算法就默认阻带误差和通带误差“同罪”阻带目标值接近0同样的绝对误差在阻带看来占比反而小算法自然懒得照顾阻带。解决按“阻带衰减需要比通带波纹高多少dB”来设定权重比例。例如要通带波纹 0.1dB、阻带衰减 60dB两者误差比约 1:1000权重就要反过来设约 1000:1。3. 系数出现不合理的异常大值滤波器输出发飘现象h[n] 序列里某些抽头值有几十甚至上百的数量级明明想要的通带增益只是 1。原因极值点之间靠得太近尤其是把过渡带边缘两点同时列为强制极值点时插值矩阵接近奇异数值解被放大。解决检查过渡带边缘点是否被重复加入候选使用重心插值版本必要时提高网格密度让相邻极值点间距不小于一个网格步长。4. 在 fs/2 处响应不为零设计结果和预期偏差很大现象用偶数抽头数Type II设计高通或带阻滤波器Nyquist频率处频率响应不为零甚至很差。原因Type II 对称FIR在 ωπ 处幅度响应被结构强制为零数学上不可能在那一点给出非零目标——再迭代也突破不了这个结构性约束。解决设计高通、带阻这类在 Nyquist 处要求非零响应的滤波器时强制使用奇数抽头数Type I或者在指标里避免在π处定义非零目标。5. 过渡带边缘位置总差几个采样点现象测出来的 -3dB 点、阻带起始频率和设计指标对不上位置整体偏移。原因频带数组的定义方式问题。实现里通带和阻带的边界必须精确落在网格点上如果换算归一化频率时用了错误的采样率或者过渡带边缘没有对齐网格算法只能近似移动边界。解决确认频带数值是按“cycles/sample”还是“rad/sample”计算的设计前先打印网格落点确认每个频带边缘对应的网格下标是整数。6. 设计取舍阶数、波纹、过渡带以及“最优”并不等于“好用”P-M算法给出了数学意义上的最优但“最优”是相对你给的权重和阶数而言的。权重如何分配、阶数取多少这些选择里藏着真正的工程判断。首先阶数是最大的取舍变量。阶数越高逼近多项式的次数越高理论能达到的最小最大误差 δ 越小但代价是指标收益递减且数值稳定性变差——高次多项式的极值对系数微小扰动极度敏感C语言里 double 的舍入误差会开始主导设计结果。实际经验是阶数一到几百尽量用浮点设计加量化分析两步走。其次等波动结构本身是一种偏好。它让最坏误差最小但代价是误差在整个频带上均匀分布。如果你需要的是“通带极度平坦阻带只要够用”等波动设计反而浪费了阻带的裕量——这时候用最小二乘设计误差能量最小得到的通带更平只是没有严格的最坏误差保证。P-M算法的“最优”只在极小极大准则下成立换一个目标函数最优解就换人了。第三加权函数是设计者唯一表达偏好的地方但权重比和实际波纹比不是线性对应。你设阻带权重是通带的 100 倍实际阻带衰减和通带波纹的差并不会正好是 40dB因为过渡带宽度、极值点分布都会影响最终结果。正确做法是先跑一轮不带权重的设计量出各频带实际误差再反推权重比迭代两三次才能对齐指标。我自己现在的习惯是设计一个新滤波器时先用一个宽松阶数跑通P-M程序打印每轮 delta 的收敛轨迹和最终极值点符号序列确认交替定理的条件确实满足了再按指标去收紧阶数。这样踩过几次上面说的坑之后我再也没让一个看起来“不收敛”的设计白白浪费过时间。这个验证习惯帮我少走了很多弯路希望帮到你。本文还有配套的精品资源点击获取