简介：一份基于Parks-McClellan算法的FIR滤波器设计工具，使用C语言实现，基于Remez交换算法，面向数字信号处理开发者和嵌入式/信号处理领域工程师，用于生成线性相位FIR滤波器系数，解决频率响应设计与误差优化问题。算法通过迭代交替交换误差极值点，使指定频段内的均方误差达到最小。压缩包共3个文件，包括C源代码、README说明文档与LGPL许可证文本，三种类型分别覆盖代码、使用说明与授权信息，整体约10KB，结构精简，便于阅读和二次开发。已有738人学习该资源。源码实现完整，可直接编译运行，也支持通过调整滤波器阶数、频带边界、过渡带宽度等参数生成不同规格的滤波器；README与注释有助于理解Parks-McClellan和Remez交换算法的核心步骤，既可用于课程学习，也可作为实际工程设计的参考实现。

1. Parks-McClellan与Remez交换：为什么最优FIR看起来是“等纹波”的

同样是用有限的抽头数实现低通，窗函数法给你的总是“某个固定形态的误差”：换窗只能换旁瓣衰减的斜率，加长度只能把所有波纹一起压小，你没法告诉算法“通带误差必须小于0.1dB、阻带必须低于-80dB”。Parks-McClellan（也就是常说的Remez交换算法）却把FIR设计变成一道加权切比雪夫逼近问题——它不承诺平均误差小，只承诺一件事：在给定阶数下，让通带和阻带的最大加权误差同时达到最小。最终误差曲线会在所有极值点上整齐地顶到同一高度，等纹波（equiripple）就是最优解的指纹。搞懂它的关键是认清两件事：为什么“等纹波”必然等价于“最优”，以及Remez迭代凭什么能在不求逆矩阵的情况下一步步逼近这个答案。

2. 加权切比雪夫逼近：把滤波器指标翻译成误差函数的数学结构

2.1 线性相位FIR：为什么设计问题只和实系数余弦多项式有关

设计FIR时，我们真正能触碰的是脉冲响应h[n]。一旦限定线性相位，h[n]就必须满足对称关系h[n]=h[M-1-n]，相位变成固定的纯延迟e^{-jω(M-1)/2}，剩下的问题只剩一个实函数的形状：幅度响应A(ω)。

以奇数长度M、类型I为例，A(ω)可以写成

A(ω) = a[0] + 2·Σ_{n=1}^{L} a[n]·cos(nω)

其中L=(M-1)/2。换句话说，设计一个线性相位FIR滤波器，等价于挑选L+1个实系数a[n]，让这个余弦多项式的图像在某段频率区间上去贴合一段理想的描迹。相位问题被完全剥离，设计空间从“M个任意复数”压缩成“L+1个实数”。这个化简是整个设计方法能数学化的前提——没有它，切比雪夫逼近在一组复指数函数上做会麻烦得多。

需要注意，偶数长度的类型II、以及奇对称的类型III/IV，只是把cos(nω)换成cos((n+1/2)ω)或sin(nω)，极值点的计数规则略有变化，但迭代机制和本章讲的完全同构。所以后面所有推导都默认类型I，结论可以平移。

2.2 加权误差函数：通带和阻带“各说各话”

工程指标不是一个“数值”，而是一段约束：通带内纹波别超过δp，阻带内衰减别小于As。要同时控制两个频段的松紧，必须允许它们的重要程度不同。于是引入两样东西：

• 理想响应D(ω)：通带内为1，阻带内为0，过渡带不定义；
• 非负权重W(ω)：通带取Wp，阻带取Ws，过渡带为0。

误差函数就定义为

E(ω) = W(ω)·(D(ω) - A(ω))

于是|E(ω)|同时编码了两件事：偏离理想响应的幅度，以及这个偏离在指标里的“代价”。如果设计目标是最小化最大|E(ω)|，那么通带内|D-A|被压缩到δ/Wp，阻带内被压缩到δ/Ws。把δ规范成1，Wp和Ws的比例就直接由纹波指标反推：

Wp/Ws = δs/δp

例如阻带要求-80dB（线性幅度约0.0001），通带允许±0.1dB（线性约0.0116），两者比值接近1:116，阻带权重就应当比通带重约116倍。这一小段推导在很多现成代码里被藏在“设计滤波器指标换算”的注释里，但它决定了整个优化问题的方向。

过渡带则是一个“灰色地带”：D和W都是0，误差恒为0，不参与博弈。Remez算法会在迭代中自动把极值点集中在通带和阻带两个区间内，过渡带附近只负责“陡降”，不承担逼近责任。

2.3 交替定理：等纹波不是巧合，是最优解的必然指纹

为什么等纹波就代表最优？核心是加权切比雪夫逼近的交替定理（alternation theorem）。

设A(ω)是上面这个L+1参数余弦多项式的集合。若A(ω)是在逼近区间I上对D(ω)的最优一致逼近，那么加权误差E(ω)在I上至少存在L+2个极值点ω_0<ω_1<...<ω_{L+1}，满足：

• |E(ω_k)| = max|E(ω)|，即这些点上的误差幅度都相等且达到全局最大；
• E(ω_k)的符号交替出现：一正一负一正一负……

这个定理给出了两个极其有用的东西。第一，它是一个验证器：如果一条误差曲线满足“L+2个等幅交替极值”，那么它就是最优的，不需要再找别的候选解。第二，它把连续问题离散化了：我们没必要在无穷多个频率点上比较，只需要找到那L+2个关键频率点。

直观上可以这样理解：一个L+1参数的余弦多项式，它的自由度刚好够“顶住”L+2个约束。误差曲线想在某个局部降下来一点，就必然在另一个极值点抬上去一点。只有当所有极值点都被顶到同一高度、谁也没法再退让时，曲线才达到帕累托边界。等纹波正是这种“处处顶满”的形态——任何一段波纹低下去，都意味着别的地方一定会有某段波纹高过头。

交替定理也解释了Remez交换算法的名字来源：算法全程在寻找极值点并“交换”极值点集合，本质上就是在逐轮逼近这个指纹。

3. Remez迭代：重心插值如何在不求逆矩阵的前提下逼近最优

3.1 极值点方程组：L+2个方程，L+2个未知数

假设我们已经猜到了一组极值点ω_0,...,ω_{L+1}。那么在最优点上，误差应当满足

A(ω_k) = D(ω_k) + (-1)^k·δ / W(ω_k)

这里的符号交错来自交替定理，δ是要解出的最大加权误差。左边是一个L+1系数的余弦多项式在ω_k处的值，右边是已知的目标值+含δ的余项。

未知数一共L+2个：L+1个系数a[n]加上δ。方程恰好也是L+2个。理论上直接解这个线性方程组就能得到A(ω)和δ。但工程上不会这么干：矩阵是稠密的，每轮迭代极值点都变，迭代几十轮就要解几十次稠密方程组，L上百时完全是浪费。

更本质的问题是，解方程不是目的，在密网格上评估误差、找出新的极值点才是。直接求解系数再逐点算A(ω)，相当于绕远路。Remez算法的经典实现利用了一个结构：切比雪夫多项式逼近中的拉格朗日插值可以写成重心形式，于是“任意频率点上的A(ω)值”可以只用极值点和极值点上的函数值表示出来。

3.2 重心插值：极值点之间的值不靠解矩阵

令x=cosω，把频率变量换成x。类型I的A(ω)变成一个x的L次代数多项式P(x)。现在有L+2个插值点x_k=cosω_k，对应值

y_k = D_k + (-1)^k·δ/W_k

对任意x，P(x)的拉格朗日插值可以写成重心形式：

P(x) = (Σ w_k·y_k/(x-x_k)) / (Σ w_k/(x-x_k))

其中重心权重

w_k = 1 / Π_{j≠k}(x_k - x_j)

w_k在极值点集合固定后只算一次，O(L^2)代价。之后在任意网格点上求P(x)，只需一次O(L)的累加，几千个网格点加起来也只是O(L·Ngrid)，远比反复解稠密矩阵划算。

δ本身也可以从方程组里分离出来，不需要显式解系数：

δ = Σ w_k·D_k / Σ w_k·(-1)^k/W_k

这个公式在很多教科书实现里直接写在一个叫remez_calc_delta的函数里。注意分母上交替的(-1)^k除以权重——这正是交替符号在代数上留下的痕迹。

3.3 交换极值点：局部极值如何变成全局最优

有了P(x)后，在预置的密集频率网格上逐点计算加权误差E(ω)=W(ω)(D(ω)-P(ω))，然后做一个简单的局部搜索：凡是误差比左右两个邻居都大（或都小）的网格点，就是候选极值点。把这些候选点按频率升序排好，替换掉上一轮的极值点集合，下一轮迭代开始。

这个“替换”就是exchange的本意。每一轮实际上在做一件事：上一轮的极值点集合只是某个局部博弈的均衡，新的搜索结果把真正的“刺头”——误差冒头的地方——换进来，把已经服帖的点踢出去。随着迭代推进，极值点逐步移向真正的最优位置，所有局部峰值被拉平到同一高度。

理想情况下，δ在整个迭代过程中单调上升并收敛：初始估计的极值点位置不够准，解出的δ低于真实最优值；每轮把更大的误差峰换进来，δ就往上拱一截，直到所有峰顶齐平。如果实现正确且网格足够密，你最终会看到δ稳定在一个值上，不再变化。

3.4 收敛判据：看δ稳定，还是看极值点集合稳定

最常见的收敛条件是相邻两轮迭代的δ相对变化小于某个阈值，比如1e-10。但这里有一个容易忽视的细节：δ只反映“当前极值点集合下的解”，如果搜索阶段漏掉了真正的极端点，δ照样可以稳定，只是稳定在一个偏小的值上。所以更可靠的判据是同时比较新旧极值点集合是否完全一致。集合不再变化，意味着没有新的刺头被找到，博弈真正终结。

我一般把两条都保留：δ相对变化小于1e-12且极值点集合完全一致，才判定收敛。如果只信δ，很容易在网格太稀或搜索逻辑有bug时收获一个“看起来收敛了但实际不达标”的设计。

4. 源码拆解：remez-exchange里的网格、极值与关键循环

4.1 数据布局：六个数组撑起整个算法

一个典型的C实现，核心数据布局并不复杂：

网格初始化时要特别注意：通带区间[0, ωp]和阻带区间[ωs, π]上都要铺点，过渡带(ωp, ωs)直接跳过。若把过渡带也放进去并赋des=0，算法会让逼近多项式尝试拟合一个在过渡带内并不存在的目标，结果往往是δ瞬间放大、迭代漂移。

网格密度是个经验参数，通常取滤波器抽头数的8～16倍。过渡带越窄，误差曲线在边界附近越陡，越需要更密的网格去捕捉极值点。

4.2 迭代骨架：初始化、解δ、插值、搜索、交换

下面是用C风格写的核心循环骨架，省略了具体的频带划分，但保留了关键路径：

#define L ((M - 1) / 2) /* 余弦多项式阶数 */ #define NP (L + 2) /* 极值点个数 */ #define NG (16 * M) /* 网格密度经验值 */ static double grid[NG], des[NG], wt[NG]; static int ext[NP]; /* 极值点在 grid 里的下标 */ static double xk[NP], yk[NP], ww[NP]; static double remez_calc_delta(const double *xk, const double *des, const double *wt, const double *ww) { double num = 0.0, den = 0.0, sign = 1.0; for (int i = 0; i < NP; ++i) { num += ww[i] * des[i]; den += ww[i] * sign / wt[i]; sign = -sign; /* 交替符号 */ } return num / den; /* 最大加权误差 δ */ } static double bary_value(double x) { double num = 0.0, den = 0.0; for (int i = 0; i < NP; ++i) { double d = x - xk[i]; if (fabs(d) < 1e-14) /* 网格点恰好就是极值点 */ return yk[i]; double t = ww[i] / d; num += t * yk[i]; den += t; } return num / den; } int remez_loop(void) { double delta = 0.0, prev_delta = 0.0; int ext_old[NP], iter; init_grid(); /* 铺网格，赋 des[] 和 wt[] */ initial_ext(); /* 在通带/阻带内均匀选 NP 个极值点 */ for (iter = 0; iter < 100; ++iter) { for (int i = 0; i < NP; ++i) { xk[i] = cos(grid[ext[i]]); yk[i] = des[ext[i]]; } bary_weights(xk, ww, NP); /* 重算重心权重 */ delta = remez_calc_delta(xk, yk, wt, ww); for (int i = 0; i < NP; ++i) yk[i] += ((i & 1) ? -delta : delta) / wt[ext[i]]; /* 在每个极值点之间用重心插值评估误差并搜索新极值 */ for (int j = 0; j < NG; ++j) if (in_band(j)) err[j] = wt[j] * (des[j] - bary_value(cos(grid[j]))); memcpy(ext_old, ext, sizeof(ext)); search_extrema(err, ext); /* 局部极值搜索 + 交换 */ if (fabs(delta - prev_delta) < 1e-12 && same_set(ext, ext_old)) break; /* δ 稳定且集合不变 */ prev_delta = delta; } return 0; }

这段代码里最值得琢磨的是bary_value的数值保护：当采样点恰好重合极值点，d=0会让除法爆炸，所以要先做一次接近判断直接返回yk[i]。这在窄过渡带场景几乎是必然触发的路径，没有这行保护，迭代后期会出现莫名其妙的分母零错误。

δ的计算也藏了一个细节：公式里wt取的是极值点处的权重。通带和阻带的权重本身差距可能上百倍，所以分母会包含符号交替的大数相消。后面避坑章节会专门回到这个问题。

4.3 search_extrema的细节：带内找峰与去重

搜索函数是Remez里最容易写错的部分。标准的做法是：只在通带与阻带区间内（in_band为真）扫描，对每个网格点比较它与左右邻居的误差值，若同时大于两侧或同时小于两侧，就记为一个局部极值点。

static void search_extrema(const double *err, int *ext) { int n = 0; for (int j = 1; j < NG - 1; ++j) { if (!in_band(j)) continue; if ((err[j] > err[j-1] && err[j] >= err[j+1]) || (err[j] < err[j-1] && err[j] <= err[j+1])) ext[n++] = j; /* 记录局部极值 */ } /* 若局部极值数量超出 NP，保留绝对误差最大的 NP 个并保持交替符号 */ prune_to_np(ext, n); }

这里有个实现取舍：搜索出的局部极值数量往往大于NP，而交替定理要求极值点恰好NP个。常见的处理是保留幅度最大的NP个，并按符号交替筛选。如果搜索得到的极值数量不足NP，则需要从已有集合中按误差下降的优先级补点。这个“补点”和“裁剪”逻辑决定算法能否在窄过渡带下稳定收敛。

4.4 从极值点反推h[n]：余弦多项式到脉冲响应

迭代收敛后，我们有了NP个极值点上的A值。接下来要把多项式还原成滤波器系数。

常见做法不是在时域硬解方程，而是利用余弦多项式的频域结构：在M个等分频率点ω_m = πm/M上对A(ω)采样，再做一次M点逆傅里叶变换。因为A(ω)是实偶余弦和，逆变换结果直接是实数且对称，居中后就是h[n]。

具体到类型I，先由插值得到a[n]，再映射：

• h[L] = a[0]
• h[L-n] = a[n]/2，n=1,...,L
• h[L+n] = a[n]/2

这一步看似简单，但实际实现里如果用均匀频率采样，需要足够多的采样点避免混叠。经验值是采样点数至少取2L以上，否则A(ω)的高频余弦分量会被折回低频，设计出的滤波器通带边缘会出现意外的纹波。这个坑在调试时常常伪装成“算法结果和理论设计曲线对不上”，其实只是最后一步采样点数不够。

5. 避坑指南：极值簇、鬼峰与网格细度的排查清单

5.1 鬼峰：网格太稀，最大误差藏在两个采样点之间

现象：算法报告δ=0.001，把h[n]拿去做FFT实测，通带边缘附近却出现0.002的凸起，纹波图形状基本等纹波、唯独某处冒出一个尖。

原因：搜索函数只在离散网格点上找极值。如果真峰值恰好落在两个网格点之间，局部比较法会把它漏掉。这个点没进极值集合，就没有参与博弈，误差自然超标。过渡带越窄，误差曲线在边界附近斜率越大，峰值越容易“藏”进网格缝隙。

解决：把网格密度翻倍重跑，对比δ。若δ明显变大，说明上一轮根本没收敛到真最优。正是这个自检逻辑，应该固化成一套习惯流程。

5.2 极值点聚堆：窄过渡带下的数值抖动与NaN

现象：设计窄过渡带低通（比如过渡带宽度只有0.01π）时，迭代跑到二三十轮，δ开始上下抖动，甚至直接变成NaN。

原因：最佳等纹波解在过渡带两侧会布置多个距离极近的极值点。极值点在x=cosω域上挨得太近，重心权重公式里的乘积(x_i - x_j)趋近于零，权重w_k变得巨大且符号剧烈翻转。double精度下，重心插值的分子分母都成了大数相消，有效位数被耗尽。

解决：在搜索后合并靠得太近的极值点，保留误差更大的那个；同时把数值累加从double换成long double，或者对重心项的累加做Kahan求和。还有一种工程派做法：过渡带边界附近额外加密网格，让极值点不会在同一网格索引上反复拉锯。

5.3 权重比悬殊：δ公式的分母为何变得不稳定

现象：阻带要求-100dB、通带只要求0.05dB纹波时，权重比能到几百比一。算法跑出来的滤波器通带纹波达标，阻带却总是差几个dB。

原因：权重悬殊让δ公式里分母的交替项量级差距巨大，符号交替的大数相减消耗大量有效数字。误差曲线的动态范围也很大，搜索阶段比较“谁更大”时，小误差一侧的相对误差被放大。

解决：先将权重归一化，令最小权重为1；再把des和wt同时缩放，使得先验估计的δ接近1。归一化之后，数值条件会明显变好。这是多数商业设计工具内部都会偷偷做的一步，原始C库往往不做，需要自己补上。

5.4 收敛判据太宽松：δ不动了，但误差峰值还在涨

现象：δ相对变化已经小于1e-9，按常规判据应该收敛。但画出误差图，maxerr比δ高出0.5%，且位置在随机变动。

原因：δ只描述“当前极值点集合”的解。如果搜索阶段漏极值，δ会稳定在偏小的值上，真正的峰值根本没进入方程组。判据选错了统计量，只盯δ而不盯极值点集合。

解决：额外比较新旧极值点集合是否完全一致，集合不动才算真正收敛。同时检查maxerr与δ的比值，若maxerr比δ高出一个容忍量（比如0.1%），直接判定网格密度不足或搜索有bug。这个比值是比δ本身更诚实的体检指标。

5.5 初始极值点落进过渡带：第一轮迭代就被拉偏

现象：首次迭代δ大得离谱（超过0.5），后面虽然慢慢回落，但总迭代次数明显偏多，结果也偶尔卡在次优解上。

原因：初始化时如果简单地在全频段[0, π]均匀取点，极值点会有一部分落在过渡带内。过渡带没有理想的D值，方程组强行让多项式去拟合这些“无目标”的位置，把整个解空间拽向错误方向。

解决：初始化只从通带和阻带两个闭区间内取点，通带阻带各占一半，按符号交替顺序铺开。搜索时也要强制跳过过渡带，不让它在后续迭代混入极值集合。这是实现Remez的第一条纪律，很多新手翻车都翻在这。

6. 设计取舍：要等纹波最优，还是要最小能量误差

等纹波最优的本质是无穷范数意义下的极小化，它保证整个逼近区间内的最坏误差最小。这个性质对通信频谱模板、抗混叠滤波器这类“验收只看最坏点”的场景极其友好：只要设计报告里说最大纹波为δ，就真的没有任何频点超限。

但代价也很明确。把误差在所有极值点拉平，意味着算法不会“怜悯”那些本来可以更松的地方。宽阻带里远离过渡带的频段其实不需要那么低的衰减，等纹波设计却在那里耗费同样的误差配额，导致总误差能量（L2范数）偏高。如果你的性能指标是信号带内的SNR或整体失真能量，最小二乘类设计反而更划算。这就是为什么工程上经常两种工具交替着用：先Remez计算模板边界是否合规，再用最小二乘做精细的带内均衡。

从计算角度看，Remez交换算法在今天依然值得在手写C代码里用，是因为它把一个半无限约束优化问题压缩成了“局部插值+极值搜索”，内存占用只有O(M)量级，没有依赖矩阵库。现代凸优化也能解决同样的问题，但会引入求解器依赖和更大的运行时开销。在嵌入式平台、实时重配滤波器参数的场景里，一个纯粹的C实现Remez通常是更稳妥的选择。

我个人的习惯是：每次修改设计指标，先跑一遍“网格翻倍自检”——网格密度翻倍重跑，对比δ和误差峰值的变动。动超过了0.1%，就说明上一轮设计根本没在真最优解上，参数表里的数字都不可信。这个习惯帮我挡掉了不少次返工，也让我对“算法报告的值”始终保留一点健康的怀疑。希望帮到你。

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