零依赖C++信号处理类：FFT单边谱标定、滤波与包络实战

简介：这份资源是一套面向 C++ 开发者与系统编程学习者的信号处理类实现代码，聚焦进程如何响应操作系统发送的各类信号这一实际问题，适合已掌握 C++ 基础语法、希望深入理解进程生命周期管理与异常恢复的中级读者参考。压缩包内共 2 个文件，分别为 1 个头文件和 1 个源文件，整体约 3KB，个头虽小但结构完整：头文件承担类声明与接口定义，源文件给出信号注册、注销与状态查询等成员函数的具体实现。代码围绕常见信号展开，涉及中断、终止、挂断、段错误等典型场景，并展示信号处理器安装、撤销与活跃状态判断的写法，同时触及 sigaction 相较 signal 在标志位控制与系统调用重启方面的差异，以及多线程环境下信号归属与线程安全的考量。已有 373 人学习，可将其作为编写健壮程序时的信号处理参考模板与调试对照素材。

1. 一帧 4096 点的振动数据，为什么值得单独养一个 C++ 信号处理类

产线上的一台旋转设备，加速度传感器接采集卡，每 10 ms 往上层丢一帧 4096 点的 int16。上位机得在半帧时间内算出 1 倍频、2 倍频幅值和总有效值，超阈值立刻报警。这种活儿真正的难点从来不是公式本身——FFT、FIR、包络都是几十年前的老东西——而是你不可能为了几个函数去背一整套第三方库的编译链，更不敢在实时线程里让某个库偷偷 new 一块内存出来。

这份 C++ 信号处理类走的是「一个头文件、一组自由函数、零外部依赖」的路子：把时域到频域的变换、常用窗函数、FIR 与双二阶 IIR 滤波、卷积与相关、重采样、希尔伯特包络，收进一个可读性还过得去的类里。适合做嵌入式上位机、测试测量软件、振动与声学分析工具的人，也适合刚学完数字信号处理、想把课本公式变成能跑起来的 C++ 的人。下文按「结构怎么搭 → 谱怎么算准 → 滤波器怎么设 → 哪儿容易翻车 → 还能怎么用」的顺序拆。

2. 资源的骨架：类结构、数据约定与 20 行接入代码

2.1 为什么不直接拉一个大库进来

先说要绕开的问题。做信号处理，市面上的选择其实不少：通用 FFT 库性能很强，一些数值计算库也自带滤波和重采样，还有专门的音频处理框架。但它们共同的代价是引入一个不小的依赖面——要么编译链复杂，要么许可证要单独评估，要么在嵌入式交叉编译时折腾半天。

我一般按三个维度做取舍。第一看帧长：帧长超过 64K、并且是离线批处理，那确实该用成熟 FFT 库，它的混合基和 SIMD 优化不是手写能追上的；帧长在 256 到 8192 之间、实时性要求高，自己写基 2 加一点循环展开就够用，还省掉一次跨库调用。第二看内存行为：实时线程里最怕的是不可预测的分配，自己掌控的类可以在构造时把缓冲区一次性 reserve 好。第三看部署：单头文件往工程里一丢就完事，这种便利在客户现场是刚需。

所以这份类的定位很明确——不做通用数值计算，只做「传感器数据处理这条流水线上最常伸的那几只手」。它不管线程调度、不管文件 IO、不管绘图，边界清楚，读起来才不累。

2.2 数据约定：三件事必须在写第一行代码前定死

新手最容易在接口上反复返工，本质是三个约定没定死：用什么精度、采样率放哪、索引怎么对应频率。

精度上，时域统一用std::vector<double>，频域统一用std::vector<std::complex<double>>。有人会问能不能用 float 省内存——可以，但只在滤波器阶数很低、采样率也不高的时候安全，低频 IIR 用 float 很容易出问题，后面第 5 章会细说。

采样率不塞进数据容器，而是作为Signal对象的成员单独存放。理由是同一个缓冲区可能被不同采样率复用（比如上采样前的原始数据），把 fs 和数据绑死反而别扭。

索引约定统一成：k = 0是直流，k = N/2是奈奎斯特频率，第 k 个 bin 对应k * fs / N赫兹。这条一旦定了就不能改，否则单边谱折叠的时候必出错。

常见接口大致是这样一张表：

提示：接口里凡是返回新vector的函数，都在内部做了一次分配。实时路径上建议改成传入输出缓冲区的重载版本，把分配挪到初始化阶段。

2.3 接入工程：CMake 与最小可运行例子

类做成头文件库最省事，不需要单独编译，交叉编译时也不会有链接问题。CMake 里这样挂：

cmake_minimum_required(VERSION 3.16) project(dsp_demo LANGUAGES CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) # 头文件库：源码全在 include 目录下，不产生 .a/.so add_library(dsp INTERFACE) target_include_directories(dsp INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/include) # MSVC 下换成 /W4 target_compile_options(dsp INTERFACE -Wall -Wextra) add_executable(demo main.cpp) target_link_libraries(demo PRIVATE dsp)

INTERFACE库的意思是「只传递头文件路径和编译选项，自己不产生二进制产物」，这是纯头文件库的标准接法。如果你的工程还要被别人当子模块引用，把add_subdirectory那层的路径写相对，别写死绝对路径。

跑通它的最小例子只有二十来行，构造一个已知成分的信号，看看谱上是不是对得上：

#include "dsp/signal.hpp" #include <cmath> #include <cstdio> #include <vector> int main() { constexpr double kPi = 3.14159265358979323846; const double fs = 48000.0; const std::size_t n = 4800; // 100 ms 数据 // 1 kHz 幅值 1.0 的正弦 + 3 kHz 幅值 0.3 的正弦 std::vector<double> x(n); for (std::size_t i = 0; i < n; ++i) { const double t = static_cast<double>(i) / fs; x[i] = 1.0 * std::sin(2.0 * kPi * 1000.0 * t) + 0.3 * std::sin(2.0 * kPi * 3000.0 * t); } dsp::Signal sig(x, fs); const auto amp = sig.amplitudeSpectrum(dsp::Window::Hann); // 打印两个峰值 bin 的频率和幅值 std::printf("bin width = %.3f Hz, bins = %zu\n", sig.binWidth(), amp.size()); for (std::size_t k = 1; k < amp.size(); ++k) { if (amp[k] > 0.05) { std::printf("f = %8.2f Hz, A = %.4f\n", k * sig.binWidth(), amp[k]); } } return 0; }

逻辑说明：Signal构造时只保存数据的引用或拷贝，不做任何变换；amplitudeSpectrum内部依次完成「去均值 → 乘窗 → FFT → 取模 → 单边折叠 → 按窗的相干增益归一化」六步。binWidth()就是fs / N，所有频率换算都基于它，改帧长不用改调用方代码。

参数说明：dsp::Window::Hann可以换成Rect、Hamming、Blackman、FlatTop；帧长 n 如果不是 2 的幂，amplitudeSpectrum内部会补零到下一个 2 的幂，但binWidth会按补零后的长度算——这一点后面第 3.4 节要专门讲清楚。

跑通这一步，环境就算通了。接下来才是真正决定结果对不对的部分。

3. 频谱分析落地：FFT 单边谱标定与窗函数怎么选

3.1 为什么你的峰值总是差一半

这是最高频的翻车点。有人拿一段幅值 1.0 的正弦做 FFT，谱上峰值读到 0.5，第一反应是"FFT 算错了"，其实十有八九是归一化没做对。

先把推导摆清楚。对长度 N 的序列x[n]，DFT 定义为X[k] = Σ x[n]·e^(-j2πkn/N)。对一个幅值 A、频率落在 bin k0 上的实正弦，理想情况下|X[k0]| = A·N/2（另一半能量在共轭的N-k0上）。所以双边谱里单看一边是 A/2，这就是"差一半"的来源。

单边谱的做法是把负频率那部分能量折回正频率，于是0 < k < N/2的幅度为2|X[k]|/N，而k = 0（直流）和k = N/2（奈奎斯特）本身没有共轭配对，不能再乘 2。

加了窗之后还要除以窗的相干增益。定义窗序列w[n]的相干增益为CG = Σw[n]，那么标定后的单边幅度谱是：

A[k] = |X[k]| / CG , k = 0 或 k = N/2 A[k] = 2·|X[k]| / CG , 0 < k < N/2

汉宁窗的CG ≈ N/2，代进去A[k] = 4|X[k]|/N，跟很多教科书的系数对上了。矩形窗CG = N，退化回2|X[k]|/N。

注意：如果你的用途是噪声功率谱密度而不是正弦幅值，归一化因子要换成等效噪声带宽ENBW = N·Σw²/(Σw)²，两者差一个窗形状因子。混用会让噪声底偏差 1.5 dB 左右，做趋势对比时看不出来，做绝对量级判定时会出事。

3.2 窗函数怎么选：主瓣宽度与旁瓣泄漏的取舍

不加窗（矩形窗）等价于把信号乘了一个矩形，频域上就是跟 sinc 做卷积。后果是强信号的能量会"漏"到邻近的弱信号 bin 上，如果两个频率相隔不到几个 bin，弱信号直接被淹没。

窗函数的本质是拿主瓣变宽换旁瓣变低。选型就看两件事：你要分辨的两个频率隔多远，以及你最弱的信号比最强的小多少 dB。

平顶窗要单独提一句：它的主瓣很宽，频率分辨能力差，但幅度平坦度极好。如果你的核心诉求是"这个频率上的幅值到底是多少"，而不是"这附近有几个频率成分"，平顶窗是正确选择，代价是相邻成分会糊成一片。

我一般的默认配置是：通用巡检用汉宁，做幅值精确测量切平顶，做整周期同步采样（比如 fs/N 正好是工频整数倍）才敢用矩形。

3.3 从时域数组到 dB 谱的完整代码

下面这段是可以直接编译验证的最小实现，包含基 2 FFT、窗生成和单边谱标定，不依赖任何外部库。

#include <complex> #include <vector> #include <cmath> #include <algorithm> #include <cassert> constexpr double kPi = 3.14159265358979323846; // 原地基 2 FFT，长度必须是 2 的幂 void fft(std::vector<std::complex<double>>& a, bool inverse = false) { const std::size_t n = a.size(); assert(n >= 2 && (n & (n - 1)) == 0); // 1) 位反转置换 for (std::size_t i = 1, j = 0; i < n; ++i) { std::size_t bit = n >> 1; for (; j & bit; bit >>= 1) j ^= bit; j ^= bit; if (i < j) std::swap(a[i], a[j]); } // 2) 蝶形运算，逐级合并 for (std::size_t len = 2; len <= n; len <<= 1) { const double ang = 2.0 * kPi / static_cast<double>(len) * (inverse ? 1.0 : -1.0); const std::complex<double> wlen(std::cos(ang), std::sin(ang)); for (std::size_t i = 0; i < n; i += len) { std::complex<double> w(1.0, 0.0); for (std::size_t k = 0; k < len / 2; ++k) { const auto u = a[i + k]; const auto v = a[i + k + len / 2] * w; a[i + k] = u + v; a[i + k + len / 2] = u - v; w *= wlen; } } } // 3) 逆变换最后除以 n if (inverse) { for (auto& v : a) v /= static_cast<double>(n); } } // 生成对称窗，n 为窗长 std::vector<double> makeWindow(int type, std::size_t n) { std::vector<double> w(n, 1.0); if (type == 0) return w; // 矩形 const double M = static_cast<double>(n - 1); for (std::size_t i = 0; i < n; ++i) { const double x = 2.0 * kPi * static_cast<double>(i) / M; switch (type) { case 1: w[i] = 0.5 - 0.5 * std::cos(x); break; // 汉宁 case 2: w[i] = 0.54 - 0.46 * std::cos(x); break; // 海明 case 3: w[i] = 0.42 - 0.5 * std::cos(x) + 0.08 * std::cos(2 * x); break; // 布莱克曼 default: break; } } return w; } // 单边幅度谱，已按窗的相干增益标定 std::vector<double> amplitudeSpectrum(const std::vector<double>& x, double fs, int winType) { std::size_t n = 1; while (n < x.size()) n <<= 1; // 补零到 2 的幂 std::vector<std::complex<double>> buf(n, {0.0, 0.0}); // 去均值，避免直流泄漏污染低频段 double mean = 0.0; for (double v : x) mean += v; mean /= static_cast<double>(x.size()); const auto w = makeWindow(winType, x.size()); double cg = 0.0; // 相干增益 for (std::size_t i = 0; i < x.size(); ++i) { buf[i] = std::complex<double>((x[i] - mean) * w[i], 0.0); cg += w[i]; } fft(buf, false); std::vector<double> amp(n / 2 + 1); amp[0] = std::abs(buf[0]) / cg; // 直流不乘 2 for (std::size_t k = 1; k < n / 2; ++k) amp[k] = 2.0 * std::abs(buf[k]) / cg; // 正频率乘 2 amp[n / 2] = std::abs(buf[n / 2]) / cg; // 奈奎斯特不乘 2 return amp; } // 转 dB，ref 是参考幅值，填 1.0 得到 dBFS 之外的绝对 dB std::vector<double> toDb(const std::vector<double>& amp, double ref = 1.0) { std::vector<double> db(amp.size()); for (std::size_t k = 0; k < amp.size(); ++k) db[k] = 20.0 * std::log10(std::max(amp[k], 1e-12) / ref); return db; }

逻辑说明：fft用标准的迭代 Cooley-Tukey，先做位反转置换，再逐级蝶形合并，全程原地操作、无额外分配。amplitudeSpectrum先把帧长补到 2 的幂，逐点减去均值后再乘窗——去均值这一步很多人省掉，结果是强直流泄漏把 0 到 3 个 bin 全糊住，低频分析直接失效。

参数说明：cg是窗的相干增益，必须用实际参与运算的x.size()个窗点求和，不能拿补零后的 n 去算，否则幅度会整体偏小。ref一般填传感器的灵敏度标定值（比如加速度计 100 mV/g 就填 0.1），这样 dB 值直接对应物理量。

想验证标定对不对，就构造一个幅值 1.0 的整周期正弦：fs = 4800、N = 4800、频率 1000 Hz，此时 1000 Hz 正好落在第 1000 个 bin 上，主瓣能量集中，汉宁窗下读出来的amp[1000]应该在 1.00 附近，误差在 0.5% 以内。如果读到 0.5 就是忘了乘 2，读到 0.25 就是既忘了乘 2 又用了1/N而非1/CG。

3.4 补零、分辨率与栅栏效应

补零能提高谱线的密度，但不能提高真实频率分辨率。真实分辨率只由观测时长决定，等于1/T = fs/N_original。补零只是在原有信息上做插值，让峰值位置看起来更细，谱线之间的间隔变成fs/N_padded。

这个区别在工程上很重要。比如你用 4800 点、48 kHz 采样，观测时长 100 ms，真实分辨率 10 Hz。补零到 65536 点之后binWidth变成 0.73 Hz，看起来精细多了，但两个相隔 5 Hz 的成分依然分不开——它们的主瓣本来就重叠在一起。

栅栏效应说的是另一件事：当信号频率不落在 bin 中心时，峰值会被"削"掉一部分。汉宁窗最坏情况下幅度误差约 1.4 dB，矩形窗能到 3.9 dB。要消除它，要么用平顶窗，要么在峰值附近做三点抛物线插值反推真实频率和幅值：

delta = 0.5 * (A[k-1] - A[k+1]) / (A[k-1] - 2*A[k] + A[k+1]) f_true = (k + delta) * binWidth A_true = A[k] - 0.25 * (A[k-1] - A[k+1]) * delta

这个公式对汉宁窗近似得很好，在信噪比尚可的场景下频率估计精度能到 0.1 个 bin 以内。我在做转速跟踪时基本都靠它，比直接取最大 bin 靠谱得多。

4. 滤波器设计：FIR 窗函数法、双二阶 IIR 与零相位滤波

4.1 FIR 还是 IIR：先看相位要求，再看算力

这个问题没有标准答案，取决于你更在乎什么。

FIR 的优点是严格线性相位，群延迟是常数(N-1)/2个采样点，所有频率成分延迟一样，波形形状不会被拉歪。缺点是同样的过渡带指标需要高得多的阶数，通常几十到几百阶，乘加运算量大。

IIR 的优点是阶数低效率高，四阶巴特沃斯就能做出很陡的过渡带，适合高采样率、低截止频率这种"归一化频率很小"的场景。缺点是非线性相位，而且极点靠近单位圆时对系数量化极其敏感。

我的判断顺序是：先把信号波形本身当结果（比如要看冲击波形的形状、要做峰值检测）就选 FIR；只关心幅值谱或者有效值这类统计量，就选 IIR。如果既想要 IIR 的低阶数又想要线性相位，那就用零相位滤波，代价是延迟翻倍、不能实时流式处理。

还有一个常被忽略的点：采样率高、截止频率低的时候（比如 48 kHz 采样滤 50 Hz），FIR 需要的阶数会大到不现实——过渡带宽度除以采样率乘以阶数，粗略估一下就知道要上千阶。这种情况老老实实用 IIR 级联。

4.2 窗函数法设计 FIR：阶数怎么估

窗函数法的思路很直白：理想低通的冲激响应是 sinc 函数，无限长且不可实现，那就截断它。截断等价于乘矩形窗，会在频域产生吉布斯振铃，所以再把矩形窗换成旁瓣更低的窗，用过渡带变宽换取阻带衰减。

截止频率为fc、采样率为fs的理想低通，其冲激响应为：

h_ideal[n] = 2·fcNorm·sinc(2·fcNorm·(n - M/2)), fcNorm = fc/fs

其中sinc(x) = sin(πx)/(πx)，M = numTaps - 1。

阶数估算有个经验公式，给定过渡带宽度Δf（归一化到 fs）和阻带衰减A（dB）：

numTaps ≈ (A - 7.95) / (2.285 · 2π · Δf) // 凯泽窗公式

海明窗固定约 53 dB 阻带衰减，过渡带宽度约3.3/numTaps（归一化），所以反过来numTaps ≈ 3.3/Δf_norm。比如 48 kHz 采样，过渡带要求 2 kHz，归一化是 0.0417，那numTaps ≈ 79，向上取到 81 保持奇数（奇数阶是 I 型线性相位，适合做低通和高通）。

// 窗函数法设计线性相位低通 FIR // numTaps: 抽头数，建议奇数；fcNorm: 归一化截止频率，取值 (0, 0.5) std::vector<double> designLowpassFIR(int numTaps, double fcNorm, int winType = 2) { assert(numTaps > 2 && (numTaps % 2) == 1); assert(fcNorm > 0.0 && fcNorm < 0.5); std::vector<double> h(numTaps); const int M = numTaps - 1; const auto w = makeWindow(winType, static_cast<std::size_t>(numTaps)); for (int n = 0; n < numTaps; ++n) { const double m = static_cast<double>(n) - M / 2.0; double sinc; if (std::abs(m) < 1e-12) sinc = 2.0 * fcNorm; // 极限值 else sinc = std::sin(2.0 * kPi * fcNorm * m) / (kPi * m); h[n] = sinc * w[n]; } // 归一化直流增益为 1，保证通带内不放大也不衰减 double sum = 0.0; for (double v : h) sum += v; for (double& v : h) v /= sum; return h; }

逻辑说明：先算出无限长理想低通的采样值，再逐点乘窗完成截断；m = 0处 sinc 是 0/0 型，直接取极限值2·fcNorm，这一行用1e-12做阈值判断，不要写成m == 0.0，浮点比较靠不住。

参数说明：fcNorm是fc/fs，不是fc/(fs/2)，很多人在这里差一倍，滤出来的截止频率直接对不上。归一化直流增益那一步不能省，否则由于窗的相干增益损失，通带会整体掉几个 dB。想设计高通、带通，用谱反转法：全通减去低通得到高通，两个低通相减得到带通，改动量很小。

4.3 双二阶 IIR 的级联与系数稳定性

高阶 IIR 千万不要直接用高阶差分方程的实现形式，也就是那个y[n] = Σb·x - Σa·y的长式子。高阶多项式的根对系数扰动极其敏感，系数量化误差会让极点跑出单位圆，输出直接发散。

正确做法是分解成一堆二阶节（SOS，second-order sections）级联，每节独立算，中间结果用 double 暂存。二阶节的系数用 RBJ 那套公式算最省事，低通的形式是：

w0 = 2π·f0/fs alpha = sin(w0)/(2Q) b0 = (1 - cos(w0))/2 b1 = 1 - cos(w0) b2 = (1 - cos(w0))/2 a0 = 1 + alpha a1 = -2·cos(w0) a2 = 1 - alpha

全部系数除以a0归一化后使用。

struct Biquad { double b0 = 1.0, b1 = 0.0, b2 = 0.0; double a1 = 0.0, a2 = 0.0; double z1 = 0.0, z2 = 0.0; // 直接 II 型转置结构的两个状态 // 单样本处理：每个样本 5 次乘法，适合实时流式 double process(double x) { const double y = b0 * x + z1; z1 = b1 * x - a1 * y + z2; z2 = b2 * x - a2 * y; return y; } void reset() { z1 = z2 = 0.0; } }; // RBJ 低通系数，f0 为截止频率，Q 为品质因数 Biquad makeLowpass(double f0, double fs, double Q = 0.7071067811865476) { const double w0 = 2.0 * kPi * f0 / fs; const double cs = std::cos(w0); const double alpha = std::sin(w0) / (2.0 * Q); const double a0 = 1.0 + alpha; Biquad bq; bq.b0 = (1.0 - cs) / 2.0 / a0; bq.b1 = (1.0 - cs) / a0; bq.b2 = (1.0 - cs) / 2.0 / a0; bq.a1 = -2.0 * cs / a0; bq.a2 = (1.0 - alpha) / a0; return bq; }

逻辑说明：这里用的是直接 II 型转置结构，状态变量z1、z2只有两个，比直接 I 型少一半存储，数值特性也更好。process每样本只做 5 次乘加，在 Cortex-M4 这类带 FPU 的 MCU 上跑 48 kHz 双声道毫无压力。

参数说明：Q = 0.7071对应巴特沃斯响应，通带最平坦；要做切比雪夫就把 Q 调高，代价是通带出现纹波。四阶滤波器就是两个二阶节级联，两级 Q 要按巴特沃斯极点分布配：Q1 = 0.5412、Q2 = 1.3066，不要两节都填 0.7071，那样响应就不对了。

还有一个坑：f0/fs很小的时候（比如 50/48000 ≈ 0.001），1 - cos(w0)会出现严重的有效位丢失，w0小到一定程度后cos(w0)几乎等于 1，减完只剩几位有效数字。这种情况用1 - cos(w0) = 2·sin²(w0/2)替代，误差能降几个数量级。这个替换我在所有用到 RBJ 公式的地方都默认开着。

4.4 零相位滤波 filtfilt：代价与延拓方式

FIR 有(N-1)/2个采样点的群延迟，IIR 的相位失真更麻烦。如果你只是想要"滤完之后波形还在原来的位置上"，那就滤两遍：正向一遍，反转再滤一遍。两次相位失真符号相反、正好抵消，总幅频响应变成单次滤波的平方。

代价有两个。第一，幅频响应被平方了，原本 -3 dB 的截止点变成 -6 dB，要补偿的话把截止频率往上调一点。第二，滤波后的有效数据要丢掉前后各一段（FIR 是numTaps-1个点），否则边界的瞬态会污染结果，不能用于严格的实时流式处理。

真正的难点在边界处理。简单反转再滤，起点的初始状态是零，会在开头产生一段明显的瞬态。稳妥的做法是做端点延拓——常用奇对称延拓（2*x[0] - x[k]反着拼在前面），或者直接复制端点值。延拓长度取滤波器阶数的 3 倍左右比较保险。

// 零相位滤波：正反各滤一遍，前后各丢 pad 个点 std::vector<double> filtfilt(const std::vector<double>& h, const std::vector<double>& x, std::size_t pad = 0) { const std::size_t p = pad ? pad : h.size(); // 默认延拓长度 = 抽头数 const std::size_t n = x.size(); // 1) 奇对称延拓，抑制边界瞬态 std::vector<double> ext; ext.reserve(n + 2 * p); for (std::size_t i = p; i > 0; --i) ext.push_back(2.0 * x[0] - x[std::min(i, n - 1)]); ext.insert(ext.end(), x.begin(), x.end()); for (std::size_t i = 1; i <= p; ++i) ext.push_back(2.0 * x[n - 1] - x[n - 1 - std::min(i, n - 1)]); // 2) 正向卷积（'same' 模式，长度与输入一致） auto convSame = [&](const std::vector<double>& s) { std::vector<double> y(s.size(), 0.0); for (std::size_t i = 0; i < s.size(); ++i) { double acc = 0.0; for (std::size_t k = 0; k < h.size(); ++k) { if (i >= k) acc += h[k] * s[i - k]; } y[i] = acc; } return y; }; auto y1 = convSame(ext); std::reverse(y1.begin(), y1.end()); // 反转 auto y2 = convSame(y1); std::reverse(y2.begin(), y2.end()); // 再反转回来 // 3) 裁掉延拓区 return std::vector<double>(y2.begin() + p, y2.begin() + p + n); }

逻辑说明：延拓让滤波器在进入有效数据区之前就已经"热"起来，状态接近稳态；正反两次卷积把相位抵消；最后裁掉两端的延拓区，剩下来的就是零相位的有效结果。

参数说明：pad默认取抽头数，实际用的时候可以加大到3 * h.size()更保险。卷积那层是 O(N·M) 的朴素实现，长序列要换成 FFT 卷积，把复杂度降到 O(N log N)，切换点大概在numTaps > 64的时候——这个阈值跟 CPU 缓存大小有关，实测一下再定。

5. 避坑与排查：五条现场血泪经验

5.1 现象：谱上峰值只有理论值的一半

现象很典型：给一个已知幅值 1.0 V 的校准正弦，谱上读出来 0.5，换了台机器还是 0.5，排除了数据源的问题。

原因只有一个方向——归一化。要么单边谱折叠时忘了给正频率乘 2，要么用1/N直接归一化而没考虑窗的相干增益CG，要么两者都犯了。还有一种隐蔽情况：CG用了补零后的长度去算，而窗其实只作用在原始长度上，导致幅度按N_orig/N_padded的比例整体缩小。

解决就是老老实实按 3.1 节的公式写：0 < k < N/2用2|X|/CG，k = 0和k = N/2用|X|/CG，CG按实际窗长求和。写完立刻用整周期正弦验证一遍，fs = 4800、N = 4800、f = 1000这组参数最省事，一个 bin 都不偏。

5.2 现象：滤波后开头两百个点全是振铃，报警逻辑误触发

现象是离线回放时前段数据一塌糊涂，后面就正常了；实时跑的时候表现为设备刚启动那几秒频繁误报。

原因有两层。一是 IIR 从零状态起步，第一个样本进来时内部状态跟稳态差得远，冲激响应还没衰减完，输出就会有一段指数型瞬态。二是因果滤波本身有群延迟，输出相对于输入整体后移了(N-1)/2个点，如果你拿滤波后的数据和原始时间戳对齐，位置就错开了。

解决分场景：离线分析直接用filtfilt，顺带把群延迟也消掉；实时流式处理就在启动阶段丢掉前5~10倍时间常数的输出，或者用稳态初值预置状态（IIR 可以用lfilter_zi那套思路，把初始状态直接设成对应第一个输入值的稳态解）。FIR 的话，丢掉前numTaps-1个点是最简单粗暴也最可靠的办法。

5.3 现象：帧长从 4096 改成 5000，程序直接崩或者结果全是 NaN

这个坑我在两个不同项目上踩过。一次是数组越界读到野内存，一次是蝶形运算的循环条件写错导致除零。

原因是基 2 FFT 对长度有硬性要求，必须是 2 的幂。5000 不是，如果代码里只写了assert而 Release 构建把断言关掉了，那就变成静默的越界写。另一类实现用递归分治，长度不对时会算出非整数索引，同样完蛋。

解决分两步。第一步是把断言改成运行时检查，长度不合法就返回错误码或者抛异常，不依赖assert——发布版本里断言是失效的。第二步是自动补零到下一个 2 的幂，同时把原始长度和补零后长度都暴露出来，让调用方知道binWidth是按哪个算的。如果确实需要任意长度的高精度变换，就上 Bluestein 算法，代价是实现复杂度和内存占用都上一个台阶。

5.4 现象：低频段滤波后输出发散，或者出现一块平的"地板"

现象是高采样率、低截止频率的场景，比如 48 kHz 采样滤 20 Hz，跑几秒之后输出开始缓慢漂移，幅度越来越大，最后变成 inf 或者 NaN。

根因是系数的有效位丢失加上极点贴近单位圆。f0/fs = 20/48000 ≈ 4.2e-4，w0极小，cos(w0)和1的差只有 8.6e-8 量级，用 float 算下来几乎全是舍入噪声带走。系数的微小误差会让极点跑到单位圆外，系统就不稳定了。

解决办法按优先级排：第一，全程用 double，不要用 float，这个改动能解决八成问题；第二，用 SOS 级联而不是高阶直连，二阶节的极点分布对量化误差容忍度高得多；第三，1 - cos(w0)换成2·sin²(w0/2)计算；第四，如果必须在 MCU 上用单精度，那就把采样率先做一次抽取降下来，别在 48 kHz 上直接滤 20 Hz。

5.5 现象：多通道数据拼在一起，时延估计总差一个采样点

做两路传感器的时延差估计，理论上应该差 3.7 个采样点，实测每次都在 3.0 或者 4.0 附近跳。换了几次互相关算法都一样。

原因往往不在互相关本身，而在前面的预处理。两路如果用了不同阶数的滤波器，群延迟就不一样；如果一路做了重采样另一路没做，相位差更是乱的；还有就是重采样之前没做抗混叠低通，混进来的高频成分会在互相关峰附近制造伪峰。

解决办法是统一预处理链：两路走完全相同的滤波器（同一个h数组，同一份 SOS 系数），需要重采样就都重采样到同一个目标采样率，重采样前先低通到新奈奎斯特频率的 0.8 倍以下。互相关峰值取整数点之后再做抛物线插值，能把精度从 1 个采样点提到 0.1 个采样点左右，这是第 6.2 节要讲的内容。

6. 进阶用法：包络解调、时延估计与一套可复现的自检流程

6.1 希尔伯特包络：把调制信息从载波里剥出来

旋转设备的早期故障特征往往不是某个频率幅值升高，而是高频载波被低频冲击调制，形成边频带。这种情况下直接看谱很难判断，得先做包络解调，把调制信号还原出来再对它做一次 FFT，看包络谱里的特征频率。

实现方式就是希尔伯特变换：对信号做 FFT，把负频率全部置零、正频率乘 2（直流和奈奎斯特保持不变），再逆变换回时域，取模就是包络。

// 希尔伯特包络，输入输出等长 std::vector<double> envelope(const std::vector<double>& x) { std::size_t n = 1; while (n < x.size()) n <<= 1; // 补零到 2 的幂 std::vector<std::complex<double>> z(n, {0.0, 0.0}); for (std::size_t i = 0; i < x.size(); ++i) z[i] = x[i]; fft(z, false); // 正频率乘 2，负频率置零，DC 与 Nyquist 不动 for (std::size_t k = 1; k < n / 2; ++k) z[k] *= 2.0; for (std::size_t k = n / 2 + 1; k < n; ++k) z[k] = 0.0; fft(z, true); std::vector<double> env(x.size()); for (std::size_t i = 0; i < x.size(); ++i) env[i] = std::abs(z[i]); return env; }

参数上唯一要留意的是n/2那两行边界：正频率范围是[1, n/2-1]，k = n/2是奈奎斯特，既不乘 2 也不置零。写成k <= n/2会让包络整体偏大一倍。另外包络的频率分辨率依然是fs/N_original，补零不会让它变好。

6.2 互相关时延估计与抛物线插值

两路信号的时延用互相关找峰最直接。整数采样点的峰值精度只有 ±0.5 个采样点，做抛物线插值能到 ±0.05 个采样点，代价只是三行代码：

// 在整点峰 k 附近做三点抛物线插值，返回亚采样偏移 double parabolicOffset(double ym, double y0, double yp) { const double denom = ym - 2.0 * y0 + yp; if (std::abs(denom) < 1e-12) return 0.0; // 峰太钝，放弃插值 return 0.5 * (ym - yp) / denom; } // 实际时延 = (k + parabolicOffset(...)) / fs

denom那个判断不能省。互相关峰如果很平坦（信号带宽窄、信噪比低），抛物线拟合会给出荒谬的偏移量，宁可返回 0 也别硬插。我一般再加一道门限：偏移量绝对值超过 0.5 就直接

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