从零实现ButterWorth滤波器：C/C++工程实践与核心算法解析

1. 项目概述：为什么我们需要亲手实现ButterWorth滤波器？

如果你正在处理音频信号、图像处理、传感器数据，或者任何需要从嘈杂背景中提取有效信息的领域，那么“滤波器”这个词对你来说一定不陌生。在众多滤波器类型中，巴特沃斯（ButterWorth）滤波器以其在通带内具有最大平坦的幅度响应特性而闻名，这意味着它在通带内几乎没有纹波，信号通过时失真极小。这个特性使得它在需要保持信号原始形状的应用中（如生物医学信号处理、音频均衡、通信系统）备受青睐。

市面上有很多现成的库，比如MATLAB的butter函数，Python的scipy.signal，它们用起来很方便。但当你需要在嵌入式系统、高性能实时处理（如音频DSP）、或者没有这些庞大运行时环境的C/C++项目中应用时，自己动手实现一个ButterWorth滤波器就成了必备技能。这不仅能让你彻底理解从模拟原型到数字实现的完整链路，更能让你对滤波器的频率响应、稳定性、计算效率有最直接的掌控。尤其是在资源受限的环境下，一个经过精心优化、没有外部依赖的C/C++实现，其价值远超一个黑盒函数调用。

这篇指南的目标，就是带你从零开始，理解巴特沃斯滤波器的设计原理，并用纯C/C++语言实现一个可配置阶数、可指定截止频率的通用IIR（无限脉冲响应）滤波器类。我们会避开复杂的数学推导，聚焦于工程实现中的关键步骤和那些容易踩坑的细节。无论你是正在学习数字信号处理的学生，还是需要在产品中集成滤波功能的工程师，这篇文章都将提供一条清晰的路径。

2. 滤波器设计核心思路与方案选型

在动手写代码之前，我们必须先理清设计一个数字巴特沃斯滤波器的完整流程。这个过程可以概括为“四步走”：确定规格、设计模拟原型、进行双线性变换、最后得到可编程的数字滤波器系数。选择用IIR形式来实现巴特沃斯滤波器，主要是因为在相同的性能要求下，IIR滤波器通常比FIR（有限脉冲响应）滤波器所需的阶数更低，计算效率更高，这对于实时处理至关重要。

2.1 从数字域指标到模拟原型

我们的起点是数字域的指标：采样频率（Fs）、截止频率（Fc）、滤波器阶数（N）。巴特沃斯滤波器的核心在于其模拟原型的传递函数。对于一个N阶低通巴特沃斯滤波器，其模拟传递函数的平方幅度响应是平坦的，其极点均匀分布在S平面左半平面的单位圆上。计算这些极点的角度是第一步关键运算。

这里的一个工程要点是，我们通常先设计一个截止频率为1 rad/s的归一化模拟低通滤波器（称为原型滤波器）。其极点位置仅由阶数N决定，计算公式为：s_k = cos(θ_k) + j * sin(θ_k)，其中θ_k = π * (2k + N - 1) / (2N)，k = 0, 1, ..., N-1。 注意，我们只取实部为负的极点（位于S左半平面），以保证滤波器的稳定性。

2.2 关键转换：双线性变换与预畸变

数字滤波器是在离散时间域工作的，我们需要将设计好的模拟滤波器“数字化”。这里最常用的方法是双线性变换。它通过公式s = (2/T) * (z-1)/(z+1)将S平面映射到Z平面，其中T是采样周期（T = 1/Fs）。

但是，双线性变换有一个致命问题：频率响应会发生扭曲，特别是高频部分。模拟频率Ω和数字频率ω之间的关系是非线性的：Ω = (2/T) * tan(ωT/2)。这意味着，如果我们直接使用数字截止频率ωc（ωc = 2π * Fc / Fs）去设计模拟滤波器，得到的数字滤波器实际截止频率会偏离预期。

解决方案是频率预畸变。我们在设计模拟滤波器时，不使用目标数字截止频率ωc，而是使用一个经过预畸变的模拟截止频率Ωc：Ωc = (2/T) * tan(ωc * T / 2)。这样，在经过双线性变换后，数字滤波器的-3dB点才会精确地落在我们想要的Fc上。这是实现中最容易忽略的一步，也是导致滤波器频率响应不准的首要原因。

2.3 IIR实现结构选型：直接II型

得到数字传递函数H(z)的分子分母多项式系数（b和a）后，我们需要选择一个递归结构来实现它。常见的有直接I型、直接II型（典范型）、级联二阶节（SOS）型。

• 直接I型：结构直观，但延迟单元多，对系数误差较敏感。
• 直接II型（典范型）：使用了最少的延迟单元（等于max(N,M)），是内存效率最高的直接形式实现。我们选择它作为基础实现，因为它在大多数情况下足够稳定且易于编码。
• 级联二阶节型：将高阶滤波器分解为多个二阶节的乘积。这是最稳健的实现方式，能最大限度地减少有限字长效应（如系数量化误差）带来的影响，特别是对于高阶滤波器。它更复杂，但数值稳定性最好。

在本指南中，我们将首先实现直接II型，因为它能最清晰地展示滤波器的差分方程。在高级部分，我们会探讨如何将其转换为更稳定的SOS形式。

3. 核心算法解析与C++类设计

接下来，我们深入到算法和代码层面。我们将设计一个ButterworthFilter类，它能够根据给定的参数（类型、阶数、采样率、截止频率）自动计算系数并完成滤波。

3.1 系数计算：从参数到差分方程

系数计算是整个滤波器的“大脑”。我们将其封装在一个calculateCoefficients函数中。以低通滤波器为例，其步骤如下：

1. 预畸变计算模拟截止频率：

   double T = 1.0 / sampleRate; double wc = 2.0 * M_PI * cutoffFreq; // 数字角频率 double Wc = (2.0 / T) * tan(wc * T / 2.0); // 预畸变后的模拟角频率

2. 计算归一化巴特沃斯极点：

   std::vector<std::complex<double>> poles; for (int k = 0; k < order; ++k) { double angle = M_PI * (2.0 * k + order - 1.0) / (2.0 * order); poles.push_back(std::complex<double>(cos(angle), sin(angle))); } // 注意：此时极点位于S平面单位圆上，实部可能为正或负。

3. 去归一化（频率缩放）： 将归一化极点s_k乘以模拟截止频率Wc，得到实际模拟滤波器的极点：p_k = Wc * s_k。

4. 双线性变换得到数字零极点： 对每一个模拟极点p_k，应用双线性变换公式z_k = (2/T + p_k) / (2/T - p_k)，得到数字极点。对于低通滤波器，我们还会在z = -1处添加N个数字零点（这是双线性变换将模拟无限大频率映射到数字域z=-1的结果）。

5. 由零极点生成传递函数系数： 将所有的数字零点和极点相乘，展开成关于z^{-1}的多项式，即可得到分母系数a（对应极点）和分子系数b（对应零点）。这个过程可以通过多项式乘法（卷积）或使用std::complex进行多项式展开来实现。对于高阶滤波器，使用级联形式是更好的选择，即直接由模拟极点p_k和零点构造二阶节，然后对每个二阶节分别进行双线性变换。

   注意：直接进行多项式展开在高阶时（如N>8）可能因数值精度问题导致系数不准确。工业级实现（如MATLAB的butter）普遍采用零极点配对、排序后转换为二阶节的方法，稳定性远超直接型。

3.2 滤波器类结构与状态管理

我们设计类的公共接口如下：

class ButterworthFilter { public: enum FilterType { LOWPASS, HIGHPASS, BANDPASS, BANDSTOP }; // 初始化滤波器，调用此函数后会计算系数 bool init(FilterType type, int order, double sampleRate, double cutoffFreq1, double cutoffFreq2 = 0.0); // 处理单个数据样本 double process(double input); // 批量处理数据（更高效） void process(const double* input, double* output, size_t numSamples); // 重置滤波器状态（用于新的数据段） void reset(); private: void calculateCoefficients(); // 核心系数计算函数 std::vector<double> b_coeffs; // 分子系数 (b0, b1, ..., bM) std::vector<double> a_coeffs; // 分母系数 (a0, a1, ..., aN)，通常a0=1 std::vector<double> x_history; // 输入延迟线 std::vector<double> y_history; // 输出延迟线 FilterType currentType; int currentOrder; double currentSampleRate; double currentCutoff1, currentCutoff2; };

状态管理（x_history和y_history）：这是IIR滤波器的记忆所在。对于直接II型结构，其差分方程为：y[n] = b0*x[n] + b1*x[n-1] + ... + bM*x[n-M] - a1*y[n-1] - ... - aN*y[n-N]x_history和y_history就是用来存储这些过去的x和y值的循环缓冲区。在process函数中，我们更新这些缓冲区并计算当前输出。

3.3 初始化与参数校验的注意事项

在init函数中，参数校验至关重要：

1. 阶数：必须是正整数。巴特沃斯滤波器阶数越高，过渡带越陡峭，但相位非线性也越严重，计算量越大。
2. 采样率与截止频率：必须满足奈奎斯特采样定理，即截止频率必须小于采样率的一半（Fc < Fs/2）。对于带通/带阻滤波器，需要两个截止频率，且必须满足Fc1 < Fc2。
3. 系数计算有效性：在calculateCoefficients中，如果计算出的极点导致分母多项式根在单位圆外（不稳定），需要返回错误。对于巴特沃斯滤波器，正确的设计流程本身能保证稳定性，但代码中仍应加入检查（如判断任何a系数是否导致极点模长>=1）。

4. 从理论到代码：分步实现指南

现在，让我们将上述设计转化为具体的C++代码。我们将重点关注低通滤波器的实现。

4.1 步骤一：搭建项目框架与基础工具函数

首先，包含必要的头文件，并定义一些数学常量和小工具函数。为了避免依赖大型数学库，我们使用C++标准库中的和。

#include <vector> #include <complex> #include <cmath> #include <algorithm> #include <stdexcept> class ButterworthFilter { public: enum FilterType { LOWPASS, HIGHPASS, BANDPASS, BANDSTOP }; // ... 公共接口声明同上 ... private: // 辅助函数：将复数根向量转换为多项式系数（仅用于低阶或理解原理，生产环境建议用SOS） static std::vector<double> polyFromRoots(const std::vector<std::complex<double>>& roots); // 辅助函数：计算多项式乘法（卷积） static std::vector<double> convolve(const std::vector<double>& a, const std::vector<double>& b); // 其他私有成员... };

4.2 步骤二：实现低通滤波器系数计算

这是最核心的部分。我们采用更稳健的模拟原型→每个二阶节分别进行双线性变换的方法。

void ButterworthFilter::calculateCoefficients() { // 清空旧系数 b_coeffs.clear(); a_coeffs.clear(); if (currentType != LOWPASS) { // 高通、带通、带阻需要不同的处理，此处先实现低通 throw std::runtime_error("Only lowpass is implemented in this example"); } int N = currentOrder; double fs = currentSampleRate; double fc = currentCutoff1; double T = 1.0 / fs; // 1. 预畸变 double wc_digital = 2.0 * M_PI * fc; double wc_analog = (2.0 / T) * tan(wc_digital * T / 2.0); // 2. 计算归一化巴特沃斯极点（仅左半平面） std::vector<std::complex<double>> analog_poles; for (int k = 0; k < N; ++k) { double angle = M_PI * (2.0 * k + N - 1.0) / (2.0 * N); // 巴特沃斯极点位于左半平面单位圆上，实部为负 std::complex<double> pole = std::complex<double>(-sin(angle), cos(angle)); // 去归一化：乘以模拟截止频率 pole *= wc_analog; analog_poles.push_back(pole); } // 3. 将复共轭极点对组合成二阶节，并分别进行双线性变换 // 对于奇数阶，会有一个实极点单独处理 b_coeffs.resize(3 * ((N + 1) / 2)); // 每个二阶节有3个b系数 a_coeffs.resize(3 * ((N + 1) / 2)); // 每个二阶节有3个a系数 int sos_index = 0; for (int i = 0; i < N; i += 2) { std::complex<double> p1 = analog_poles[i]; std::complex<double> p2 = (i + 1 < N) ? analog_poles[i + 1] : std::conj(p1); // 如果是奇数阶最后一个，用其自身共轭（实极点） // 对于实极点，虚部为0，p2 = p1 if (std::abs(p1.imag()) < 1e-10) { p2 = p1; } // 双线性变换应用于模拟传递函数 H(s) = 1 / ((s - p1)(s - p2)) // 经过推导，数字二阶节系数公式如下： double a0 = 1.0; double a1, a2; double b0, b1, b2; if (std::abs(p1.imag()) < 1e-10) { // 处理实极点（一阶节） // 实际上，一阶节可以视为二阶节的特例，其中某些系数为0 // 更简单的做法：将一阶节合并到上一个二阶节，或单独处理。 // 此处为简化，我们假设N为偶数。奇数阶需要额外处理。 // 下面代码按二阶节通用公式编写，对实极点也适用。 } // 通用二阶节双线性变换系数计算 // 参考：数字滤波器设计（Parks & Burrus）或常用公式 // 这里使用一个经过整理的常用公式集： double T2 = T * T; double p1r = p1.real(), p1i = p1.imag(); double p2r = p2.real(), p2i = p2.imag(); // 计算中间变量，用于求解系数 // 注意：此公式假设模拟传递函数为 H(s) = 1 / ( (s-p1)(s-p2) ) // 经过双线性变换 s = (2/T)*(z-1)/(z+1) 后，整理得到数字滤波器系数 // 具体推导过程略，以下是结果系数： double beta1 = - (p1r + p2r); double beta0 = p1r*p2r - p1i*p2i; double A = 4.0 + 2.0*T*beta1 + T2*beta0; double B = -8.0 + 2.0*T2*beta0; double C = 4.0 - 2.0*T*beta1 + T2*beta0; b0 = 1.0 / A; b1 = 2.0 / A; b2 = 1.0 / A; a1 = B / A; a2 = C / A; // 存储系数 [b0, b1, b2, a0, a1, a2]，但通常a0归一化为1 b_coeffs[sos_index * 3 + 0] = b0; b_coeffs[sos_index * 3 + 1] = b1; b_coeffs[sos_index * 3 + 2] = b2; a_coeffs[sos_index * 3 + 0] = 1.0; // a0 a_coeffs[sos_index * 3 + 1] = a1; a_coeffs[sos_index * 3 + 2] = a2; sos_index++; } // 调整大小以匹配实际的二阶节数量 b_coeffs.resize(sos_index * 3); a_coeffs.resize(sos_index * 3); }

重要提示：上面的系数计算公式是经过简化的示例。在实际的工业级实现中，系数计算库（如scipy.signal.butter）使用更为精确和稳定的算法，例如通过将模拟极点/零点直接代入双线性变换公式，并求解线性方程组来得到系数。上述代码旨在展示流程。对于关键应用，建议使用已验证的算法库或直接参考权威教材中的公式。

4.3 步骤三：实现直接II型滤波函数

有了系数，我们就可以实现滤波函数了。我们采用级联二阶节（SOS）的直接II型实现，这是最推荐的方式。

double ButterworthFilter::process(double input) { if (b_coeffs.empty() || a_coeffs.empty()) { throw std::runtime_error("Filter not initialized. Call init() first."); } double output = input; size_t num_sos = b_coeffs.size() / 3; // 假设系数以[b0,b1,b2,a0,a1,a2]...方式存储 // 我们需要为每个二阶节维护独立的状态变量（延迟单元） // 这里简化处理，假设类内部有 `std::vector<std::vector<double>> sos_state;` 存储每个二阶节的w[n-1], w[n-2] // 初始化： if (x_sos_history.size() != num_sos) { x_sos_history.resize(num_sos, std::vector<double>(2, 0.0)); } for (size_t i = 0; i < num_sos; ++i) { double b0 = b_coeffs[i*3]; double b1 = b_coeffs[i*3 + 1]; double b2 = b_coeffs[i*3 + 2]; double a1 = a_coeffs[i*3 + 1]; // a0=1 double a2 = a_coeffs[i*3 + 2]; // 获取当前二阶节的状态 double &w_n1 = sos_state[i][0]; // w[n-1] double &w_n2 = sos_state[i][1]; // w[n-2] // 直接II型（典范型）差分方程： // w[n] = x[n] - a1*w[n-1] - a2*w[n-2] // y[n] = b0*w[n] + b1*w[n-1] + b2*w[n-2] double w_n = output - a1 * w_n1 - a2 * w_n2; double y_n = b0 * w_n + b1 * w_n1 + b2 * w_n2; // 更新状态 w_n2 = w_n1; w_n1 = w_n; output = y_n; // 当前节的输出是下一节的输入 } return output; } void ButterworthFilter::process(const double* input, double* output, size_t numSamples) { for (size_t i = 0; i < numSamples; ++i) { output[i] = process(input[i]); } }

4.4 步骤四：测试与验证

实现完成后，必须进行测试。最直观的方法是生成一个测试信号（如正弦波加白噪声），观察滤波前后的时域波形和频域频谱。

// 简单的测试示例 void testLowpass() { ButterworthFilter filter; int order = 4; double fs = 1000.0; // 采样率 1kHz double fc = 50.0; // 截止频率 50Hz if (!filter.init(ButterworthFilter::LOWPASS, order, fs, fc)) { std::cerr << "Filter initialization failed!" << std::endl; return; } // 生成一个包含10Hz和100Hz成分的信号 size_t numSamples = 1000; std::vector<double> signal(numSamples); std::vector<double> filtered(numSamples); for (size_t i = 0; i < numSamples; ++i) { double t = i / fs; signal[i] = sin(2.0 * M_PI * 10.0 * t) + 0.5 * sin(2.0 * M_PI * 100.0 * t); } // 滤波 filter.process(signal.data(), filtered.data(), numSamples); // 此处可以输出信号和滤波后信号，用绘图工具（如Python的matplotlib）查看 // 预期：100Hz成分被大幅衰减，10Hz成分基本保留。 std::cout << "Original signal sample: " << signal[100] << std::endl; std::cout << "Filtered signal sample: " << filtered[100] << std::endl; }

更专业的验证包括：

1. 频率响应测试：输入单位冲激信号，计算其FFT得到频率响应，绘制幅频和相频特性曲线，检查-3dB点是否在fc处，通带是否平坦。
2. 阶跃响应测试：观察滤波器的瞬态响应和过冲。
3. 稳定性测试：长时间运行，输入随机噪声，观察输出是否发散。

5. 高级话题、优化与避坑指南

一个基础的滤波器实现完成后，要将其用于实际项目，还需要考虑很多工程细节。

5.1 处理高通、带通与带阻滤波器

低通是基础。其他类型的滤波器可以通过频率变换从低通原型得到。

• 高通：在模拟域，将低通传递函数中的s替换为1/s（即倒数频率变换），然后再进行双线性变换。在代码实现上，更常用的方法是：先按低通设计，然后对计算出的数字低通系数进行变换。对于二阶节，高通变换相当于将分子系数[b0, b1, b2]变为[b2, -b1, b0]（同时可能需要对整体增益进行调整）。
• 带通/带阻：涉及从低通到带通的频率变换，更为复杂。通常需要指定中心频率和带宽。一个可靠的方法是先设计一个模拟低通原型，然后通过模拟频率变换得到模拟带通/带阻滤波器，最后再进行双线性变换。强烈建议参考权威的滤波器设计教材（如《Discrete-Time Signal Processing》 by Oppenheim）或直接借鉴成熟开源库（如scipy.signal的iirfilter函数）的变换方法。

5.2 数值稳定性与有限字长效应

这是IIR滤波器在嵌入式或定点DSP中实现时的核心挑战。

• 系数量化误差：双精度系数在存入16位或32位定点数时会被舍入。这可能会轻微改变极点位置，甚至将原本稳定的极点推到单位圆外，导致滤波器振荡发散。
  • 对策：使用级联二阶节（SOS）形式。将高阶滤波器分解为多个二阶节串联，每个二阶节单独量化，能极大降低误差敏感性。在计算系数时，就应直接生成SOS系数。
• 运算溢出：在定点运算中，中间结果可能超出表示范围。
  • 对策：进行定标分析。估算滤波器在最大输入下的输出范围，在运算中适时使用缩放因子。例如，可以将输入信号右移几位（相当于除以2^N）再进行滤波，最后再将输出左移回来。
• 极限环振荡：由于舍入误差，即使输入为零，输出也可能在一个小范围内持续振荡。
  • 对策：使用饱和运算（限制输出范围）和死区处理（当信号绝对值小于某个阈值时强制为零）。

5.3 实时性与性能优化

• 减少计算量：直接II型结构已经是最小延迟单元结构。进一步优化包括：
  • 将系数和状态变量声明为float而非double（如果精度允许），在ARM Cortex-M等平台上能显著提升速度。
  • 使用编译器优化（如-O3），并确保循环内部代码简洁。
  • 对于已知的固定系数滤波器，可以手动展开循环或使用模板元编程在编译期计算部分表达式。
• 批量处理：我们提供的process(const double* in, double* out, size_t n)函数本身就有利于缓存优化，比循环调用单样本process更高效。
• 使用SIMD指令：在x86（SSE/AVX）或ARM（NEON）平台上，可以使用单指令多数据流指令并行处理多个数据或同时计算多个二阶节，能获得数倍的性能提升。但这需要平台相关的内联汇编或 intrinsics 代码。

5.4 常见问题排查速查表

5.5 从浮点到定点：嵌入式实现的额外考量

如果你需要在没有硬件浮点单元（FPU）的MCU上运行，浮点运算会非常缓慢。这时需要定点数实现。

1. 系数定点化：将计算好的浮点系数乘以一个很大的常数（如2^15=32768，即Q15格式），然后四舍五入到最接近的整数。存储为int16_t或int32_t。
2. 运算定点化：所有乘法和加法都使用整数运算。乘法后通常需要右移来保持定点数格式。

   // 示例：Q15格式定点运算 (假设系数和状态都已转为Q15) int32_t w_n = (input << 15) - ((a1_q15 * w_n1) >> 15) - ((a2_q15 * w_n2) >> 15); // 注意移位和精度 int32_t y_n = ((b0_q15 * w_n) >> 15) + ((b1_q15 * w_n1) >> 15) + ((b2_q15 * w_n2) >> 15); // 更新状态，可能需要将w_n饱和处理到16位范围后再存储 w_n2 = w_n1; w_n1 = saturate_q15(w_n); // 自定义饱和函数 output = y_n >> 15; // 将输出转换回原尺度

3. 噪声与精度权衡：定点化会引入量化噪声和舍入误差。需要仔细选择Q格式（如Q15、Q31），在动态范围和精度之间取得平衡。通常，在音频处理中，Q31格式能提供很好的性能。

实现一个工业级的ButterWorth滤波器，从理解原理到写出稳健的C/C++代码，是一个系统工程。它涉及信号处理理论、数值计算和软件工程的交叉。本文提供的指南和代码框架是一个坚实的起点。当你成功地将一个自己实现的滤波器应用于实际项目，并看到嘈杂的信号变得平滑清晰时，那种成就感是调用现成库函数无法比拟的。最重要的是，这个过程赋予了你根据具体需求（计算资源、精度要求、实时性）去调整和优化滤波器的能力，这才是工程师的核心价值所在。如果在实现过程中遇到本文未覆盖的特定问题，深入分析频率响应曲线、检查极点位置、以及用简单的测试信号（如冲激、阶跃）进行验证，永远是调试滤波器最有效的方法。