★ FEATURED ARTICLE 简介这份资源是一套面向 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::vectordouble频域统一用std::vectorstd::complexdouble。有人会问能不能用 float 省内存——可以但只在滤波器阶数很低、采样率也不高的时候安全低频 IIR 用 float 很容易出问题后面第 5 章会细说。采样率不塞进数据容器而是作为Signal对象的成员单独存放。理由是同一个缓冲区可能被不同采样率复用比如上采样前的原始数据把 fs 和数据绑死反而别扭。索引约定统一成k 0是直流k N/2是奈奎斯特频率第 k 个 bin 对应k * fs / N赫兹。这条一旦定了就不能改否则单边谱折叠的时候必出错。常见接口大致是这样一张表方法输入输出关键约定timeDomain()—const std::vectordouble只读长度 帧长spectrum(win)窗类型vectorcomplexdouble双边谱未归一化amplitudeSpectrum(win)窗类型vectordouble单边谱已按相干增益标定powerDb(win, ref)窗类型、参考值vectordouble20*log10(A/ref)firFilter(h, mode)系数、模式vectordouble模式含Full/Same/ValidsosFilter(sos, mode)二阶节数组vectordouble级联双二阶抗溢出envelope()—vectordouble希尔伯特包络长度不变resample(newFs)目标采样率vectordouble内部先低通再插值提示接口里凡是返回新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::vectordouble x(n); for (std::size_t i 0; i n; i) { const double t static_castdouble(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。窗主瓣宽度bin最高旁瓣幅度误差典型场合矩形2-13 dB最大整周期采样、瞬态分析汉宁4-31 dB中通用谱分析默认选它海明4-43 dB中语音、需要压近旁瓣布莱克曼6-58 dB较大强弱信号差 60 dB 以上平顶约 5-90 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::vectorstd::complexdouble 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_castdouble(len) * (inverse ? 1.0 : -1.0); const std::complexdouble wlen(std::cos(ang), std::sin(ang)); for (std::size_t i 0; i n; i len) { std::complexdouble 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_castdouble(n); } } // 生成对称窗n 为窗长 std::vectordouble makeWindow(int type, std::size_t n) { std::vectordouble w(n, 1.0); if (type 0) return w; // 矩形 const double M static_castdouble(n - 1); for (std::size_t i 0; i n; i) { const double x 2.0 * kPi * static_castdouble(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::vectordouble amplitudeSpectrum(const std::vectordouble x, double fs, int winType) { std::size_t n 1; while (n x.size()) n 1; // 补零到 2 的幂 std::vectorstd::complexdouble buf(n, {0.0, 0.0}); // 去均值避免直流泄漏污染低频段 double mean 0.0; for (double v : x) mean v; mean / static_castdouble(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::complexdouble((x[i] - mean) * w[i], 0.0); cg w[i]; } fft(buf, false); std::vectordouble 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; } // 转 dBref 是参考幅值填 1.0 得到 dBFS 之外的绝对 dB std::vectordouble toDb(const std::vectordouble amp, double ref 1.0) { std::vectordouble 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[k1]) / (A[k-1] - 2*A[k] A[k1]) f_true (k delta) * binWidth A_true A[k] - 0.25 * (A[k-1] - A[k1]) * delta这个公式对汉宁窗近似得很好在信噪比尚可的场景下频率估计精度能到 0.1 个 bin 以内。我在做转速跟踪时基本都靠它比直接取最大 bin 靠谱得多。4. 滤波器设计FIR 窗函数法、双二阶 IIR 与零相位滤波4.1 FIR 还是 IIR先看相位要求再看算力这个问题没有标准答案取决于你更在乎什么。FIR 的优点是严格线性相位群延迟是常数(N-1)/2个采样点所有频率成分延迟一样波形形状不会被拉歪。缺点是同样的过渡带指标需要高得多的阶数通常几十到几百阶乘加运算量大。IIR 的优点是阶数低效率高四阶巴特沃斯就能做出很陡的过渡带适合高采样率、低截止频率这种归一化频率很小的场景。缺点是非线性相位而且极点靠近单位圆时对系数量化极其敏感。我的判断顺序是先把信号波形本身当结果比如要看冲击波形的形状、要做峰值检测就选 FIR只关心幅值谱或者有效值这类统计量就选 IIR。如果既想要 IIR 的低阶数又想要线性相位那就用零相位滤波代价是延迟翻倍、不能实时流式处理。还有一个常被忽略的点采样率高、截止频率低的时候比如 48 kHz 采样滤 50 HzFIR 需要的阶数会大到不现实——过渡带宽度除以采样率乘以阶数粗略估一下就知道要上千阶。这种情况老老实实用 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和阻带衰减AdBnumTaps ≈ (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::vectordouble designLowpassFIR(int numTaps, double fcNorm, int winType 2) { assert(numTaps 2 (numTaps % 2) 1); assert(fcNorm 0.0 fcNorm 0.5); std::vectordouble h(numTaps); const int M numTaps - 1; const auto w makeWindow(winType, static_caststd::size_t(numTaps)); for (int n 0; n numTaps; n) { const double m static_castdouble(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的长式子。高阶多项式的根对系数扰动极其敏感系数量化误差会让极点跑出单位圆输出直接发散。正确做法是分解成一堆二阶节SOSsecond-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.0011 - 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::vectordouble filtfilt(const std::vectordouble h, const std::vectordouble x, std::size_t pad 0) { const std::size_t p pad ? pad : h.size(); // 默认延拓长度 抽头数 const std::size_t n x.size(); // 1) 奇对称延拓抑制边界瞬态 std::vectordouble 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::vectordouble s) { std::vectordouble 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::vectordouble(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|/CGk 0和k N/2用|X|/CGCG按实际窗长求和。写完立刻用整周期正弦验证一遍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-4w0极小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::vectordouble envelope(const std::vectordouble x) { std::size_t n 1; while (n x.size()) n 1; // 补零到 2 的幂 std::vectorstd::complexdouble 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::vectordouble 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(...)) / fsdenom那个判断不能省。互相关峰如果很平坦信号带宽窄、信噪比低抛物线拟合会给出荒谬的偏移量宁可返回 0 也别硬插。我一般再加一道门限偏移量绝对值超过 0.5 就直接本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
