简介：这是一份用 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的线性相位FIR（N为奇数，Type I对称），其幅度响应 A(ω) 可以写成 cos(ω) 的 M 次三角多项式，其中 M = (N−1)/2。于是设计问题归结为：找一个 M 次多项式，去加权逼近给定的目标函数 D(ω)，使最大加权误差最小。
这一节的数学结论是交替定理：若加权误差 E(ω) 在频带上有至少 M+2 个极值点，并且这些极值点处误差的符号正负交替、幅度都等于同一个值 δ，那么这个多项式就是唯一的最佳逼近，δ 就是能达到的最小最大误差。
这个“M+2个交替符号的极值点”就是最优解的指纹。直观理解：一个有 M 个自由参数的逼近，误差曲线如果只有 M+1 个极值，你还能调整参数把最高的峰压下去、导致另一个峰升高，直到峰与峰拉平。只有交替峰数达到 M+2，系统才“卡死”——任何微调都会让某个已经等于理论下限的峰更高，最优性被锁定了。注意极值点必须同时覆盖通带和阻带，两个频带的边界（过渡带边缘）天然是候选极值点。
3.2 Remez交换的四个步骤：网格、插值、误差、换点
交替定理给了判据，但怎么找到这组极值点？Remez交换算法是一个不断“猜测→验证→修正”的迭代过程：
- 初始化：在频率网格上等间隔取 M+2 个点作为初始极值点集合；
- 插值：假设这组极值点处误差都等于 ±δ（符号交替），用拉格朗日插值反解出当前的多项式系数和 δ 值；
- 求误差：在非常密集的全网格上计算加权误差 E(ω)，找出所有局部极值点；
- 交换：从这些局部极值点中选出绝对值最大且符合符号交替要求的 M+2 个点，替换掉旧极值点；
- 重复 2～4，直到连续两轮的极值点集合不再变化、δ 的波动小于阈值。
为什么这个循环能收敛？关键在于第 2 步解出的 δ 一定小于上一轮的真实最大误差——因为上一轮的极值点并非真正的等波动点，插值得到的“假定δ”不是全局最大，被低估了。第 3 步在全网格上扫描，必然能找到超过这个 δ 的新极值点，把它们换进来，下一轮的 δ 就会增大。这个单调过程最终收敛到理论最小最大误差。实际程序里通常几轮迭代 δ 就稳定了，Remez交换的收敛速度在工程上完全够用。
4. 核心代码拆解：C语言实现里的网格、插值与极值点搜索
一个典型的P-M算法C实现，主循环看起来并不长，难在几个细节函数。下面这段是主循环的核心骨架，注释标出了每步的机制。
/* remez 主迭代 * L = (N-1)/2, N为抽头数 * ext[] 存极值点在网格上的下标，大小 L+2 */ 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承担的是交替定理里“假定所有极值点误差相等”的角色：它把当前 M+2 个极值点看作已经处于等波动状态，用插值倒推多项式在某个参考点处的值和 δ。这里用重心形式的拉格朗日插值（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 的收敛轨迹和最终极值点符号序列，确认交替定理的条件确实满足了，再按指标去收紧阶数。这样踩过几次上面说的坑之后，我再也没让一个看起来“不收敛”的设计白白浪费过时间。这个验证习惯帮我少走了很多弯路，希望帮到你。
本文还有配套的精品资源，点击获取