ARTICLE DETAIL

资讯详情

深耕网站建设、视觉设计与SEO优化的一线实战洞察。

嵌入式开发常用滤波算法C语言实现:从原理到代码优化

嵌入式开发常用滤波算法C语言实现:从原理到代码优化 1. 项目缘起为什么嵌入式开发绕不开滤波算法干了这么多年嵌入式从8位单片机玩到ARM Cortex-M/A系列再到各种RTOS和Linux我越来越觉得算法是嵌入式的灵魂而滤波器绝对是灵魂里最基础、最常用、也最考验功力的那一部分。你可能觉得不就是几个公式、几行代码吗但真到了项目里一个滤波器的选型、参数设计和代码实现直接关系到整个系统的稳定性、响应速度和精度。比如你用ADC采集一个传感器的温度值那读数是不是总在跳你的电机控制PWM输出是不是偶尔会有毛刺导致抖动这些问题的背后往往都需要一个合适的滤波器来“熨平”数据。网上关于滤波器的理论文章、MATLAB仿真一大堆但一落到嵌入式C语言实现上很多新手就懵了。是直接用浮点数算还是定点数用数组缓存历史数据会不会太占内存中断服务函数里调用滤波函数计算时间会不会超这些问题教科书不会细讲但却是我们一线开发者每天都要面对的。所以今天我就结合自己踩过的坑和项目经验把嵌入式里最常用的几种滤波器的C语言实现从原理到代码再到优化技巧给你掰开揉碎了讲清楚。无论你是做智能小车、物联网终端还是工业控制这篇文章里的“干货”你都能直接拿去用。2. 滤波器的“第一性原理”从需求到选型在动手写代码之前我们必须搞清楚一个核心问题我到底需要什么样的滤波器这不是拍脑袋决定的而是由你的信号特性和系统需求决定的。2.1 信号与噪声你需要滤掉什么嵌入式系统处理的信号大体可以分为两类缓变信号和快变信号。缓变信号比如温度、压力、电池电压快变信号比如音频、振动、电机转速。噪声也分很多种工频干扰50/60Hz、高频开关噪声、随机白噪声还有因为传感器本身或电路引起的脉冲干扰偶尔出现的尖峰。选型的第一步就是分析你的信号和噪声在频率上的关系。如果噪声的频率明显高于有用信号比如高频开关噪声叠加在直流电压上那么一个低通滤波器就是你的首选。如果噪声是某个特定的频率比如工频干扰那么陷波滤波器可能更合适。如果你需要的是一个平滑的、反应“趋势”的值而不是每一个跳变的瞬时值那么各种平滑滤波器也叫软件滤波器就派上用场了。2.2 经典滤波器家族IIR与FIR的抉择理论上的滤波器分为无限脉冲响应和有限脉冲响应两大类。对于嵌入式开发者你可以这样简单理解IIR滤波器它的输出不仅与输入有关还与过去的输出有关有反馈。优点是效率高用较少的阶数就能实现很陡的衰减特性比如巴特沃斯、切比雪夫滤波器。缺点是相位响应是非线性的可能会引起信号失真且稳定性需要仔细设计。FIR滤波器它的输出只与有限个过去的输入有关无反馈。优点是绝对稳定且可以设计成具有线性相位保证信号波形不失真。缺点是要达到好的滤波效果通常需要更高的阶数计算量更大。在资源紧张的MCU上对于实时性要求高、且对相位要求不严苛的场景比如单纯的去噪IIR往往是更实际的选择。而对于音频处理、通信解调等对波形保真度要求高的场景即使计算量大也常选择FIR。2.3 嵌入式实现的特殊约束理论很美但嵌入式很“骨感”。我们必须考虑三大约束计算资源MCU的主频、是否有硬件乘法器、DSP指令集、FPU浮点单元。这决定了你能跑多复杂的算法。内存资源RAM大小。滤波器需要存储历史数据状态变量阶数越高内存占用越大。实时性采样频率是固定的你必须保证在两次采样间隔内完成本次滤波计算否则数据就会堆积系统崩溃。基于这些约束我们在C语言实现时会做出很多妥协和优化比如大量使用定点数运算、查找表、牺牲精度换速度等。下面我们就进入实战环节看看这些滤波器在C语言里究竟长什么样。3. 轻量级平滑滤波器新手必备的“三板斧”这些滤波器简单、计算量小非常适合新手入门和在资源极其有限的8位/16位MCU上使用用于处理缓变信号。3.1 限幅滤波法对付“野点”的硬汉这可能是最简单的滤波了。原理就是根据经验设定一个最大允许偏差值。如果本次采样值与前一次有效值的差值超过了这个偏差就认为本次是干扰脉冲舍弃它用上次值代替。#define A 10 // 最大允许偏差值 int value; // 上次有效值 int new_value; // 本次采样值 int LimiterFilter(int new_sample) { static int last_valid_value 0; int ret; if (last_valid_value 0) { // 第一次调用直接使用 last_valid_value new_sample; return new_sample; } if ((new_sample - last_valid_value A) || (last_valid_value - new_sample A)) { ret last_valid_value; // 超出限幅取上次值 } else { ret new_sample; // 正常值更新 last_valid_value new_sample; } return ret; }注意这个方法的“A”值很关键。设小了会把正常的信号变化也滤掉设大了又滤不掉干扰。它只能滤除偶然出现的、幅度大的脉冲干扰对于周期性或小幅噪声无能为力。3.2 中位值滤波法排序找“中坚”连续采样N次N取奇数把这N个值从小到大排序取中间的那个值作为本次有效值。这种方法对脉冲干扰和偶然的抖动有很好的抑制效果。#define N 11 // 采样次数建议为奇数 int MedianFilter(int new_sample) { static int data_buf[N]; static int data_index 0; int temp_buf[N]; int i, j, temp; // 1. 循环存储新数据 data_buf[data_index] new_sample; data_index (data_index 1) % N; // 2. 复制数据到临时数组进行排序避免破坏原数组 for (i 0; i N; i) { temp_buf[i] data_buf[i]; } // 3. 使用冒泡排序简单N小时可用 for (i 0; i N - 1; i) { for (j 0; j N - 1 - i; j) { if (temp_buf[j] temp_buf[j 1]) { temp temp_buf[j]; temp_buf[j] temp_buf[j 1]; temp_buf[j 1] temp; } } } // 4. 返回中值 return temp_buf[N / 2]; }实操心得N的选取很重要。N越大滤波效果越好但滞后也越严重且排序耗时剧增冒泡排序时间复杂度O(N²)。在实时性要求高的中断里如果N较大建议使用更高效的排序算法如希尔排序或者将排序过程放在主循环中异步进行。我曾在一个电机电流采样的项目里用中值滤波滤除PWM开关引起的尖峰效果立竿见影但N选了5因为排序5个数很快。3.3 递推平均滤波法滑动平均滤波平滑的代价连续取N个采样值作为一个队列每次采样到一个新数据就放入队尾并扔掉队首的一个老数据然后计算当前队列里N个数据的算术平均值。它对于周期性干扰有良好的抑制作用平滑度高但灵敏度低对偶然出现的脉冲干扰抑制作用差。#define N 12 // 队列长度 int MovingAverageFilter(int new_sample) { static int sum 0; static int data_buf[N]; static int data_index 0; // 减去即将被移出的旧数据 sum - data_buf[data_index]; // 存入新数据 data_buf[data_index] new_sample; // 加上新数据 sum new_sample; // 更新索引 data_index (data_index 1) % N; // 返回平均值 return sum / N; }技巧与坑点这是滑动平均的高效实现。传统做法是每次求和N个数计算量是O(N)。而这个实现维护了一个总和sum每次更新只需做一次减法和一次加法计算量是O(1)非常适合在中断中调用。但要注意初始化和sum的溢出问题。sum和data_buf最好用int32_t即使采样值是int16_t因为累加N次可能溢出。队列长度N决定了平滑度和滞后性需要权衡。4. 一阶低通数字滤波器从模拟到数字的桥梁这是嵌入式中最最常用的IIR滤波器没有之一。它计算简单效果直观广泛应用于信号去噪和软件抗混叠。4.1 公式推导与物理意义一阶低通滤波器的传递函数来自模拟RC电路H(s) 1 / (τs 1)其中τRC是时间常数。通过后向差分法离散化我们可以得到其数字形式的差分方程Y(n) α * X(n) (1 - α) * Y(n-1)其中X(n)是本次采样值输入。Y(n)是本次滤波输出值。Y(n-1)是上次滤波输出值。α是滤波系数α ΔT / (τ ΔT)ΔT是采样周期。这个公式的物理意义非常直观本次的输出是本次输入和上次输出的加权平均。α越大本次输入的权重越大滤波器响应越快但平滑效果差α越小上次输出的权重越大滤波效果越平滑但滞后越严重。4.2 定点数实现MCU上的速度与精度平衡在MCU上做浮点乘法非常慢尤其是没有FPU的情况下。我们必须使用定点数运算。思路是把小数α放大为整数运算完后再缩小。假设我们使用Q格式定点数。例如使用Q15格式16位有符号整数小数点在第15位之后那么1.0就用327670x7FFF来表示。#define SAMPLE_TIME_MS 10 // 采样周期10ms #define TIME_CONSTANT_MS 100 // 时间常数100ms #define ALPHA_Q15 (int16_t)((float)(SAMPLE_TIME_MS) / (SAMPLE_TIME_MS TIME_CONSTANT_MS) * 32768) // 计算α的Q15值 int16_t FirstOrderLowPassFilter_Q15(int16_t new_sample) { static int32_t y_prev_q15 0; // 上次输出值用32位存储中间结果防溢出 int32_t y_current_q15; // Y(n) α * X(n) (1-α) * Y(n-1) // 注意 (1-α) 在Q15下等于 32768 - ALPHA_Q15 y_current_q15 ((int32_t)ALPHA_Q15 * new_sample (32768 - ALPHA_Q15) * y_prev_q15) 15; y_prev_q15 y_current_q15; return (int16_t)(y_current_q15); // 转回16位输出 }关键解析ALPHA_Q15是提前计算好的常数避免了运行时做浮点除法。乘法ALPHA_Q15 * new_sample的结果是32位的Q30格式。(32768 - ALPHA_Q15)计算的是(1-α)的Q15表示。注意这里是32768因为1.0在Q15下是32768不对这里有个细节在Q15中能表示的最大数略小于10.999969...即32767。但(1-α)的运算我们通常在整数域做(115) - ALPHA_Q15即32768 - ALPHA_Q15得到的是Q15格式的(1-α)。两个Q15数相乘结果是Q30右移15位变回Q15。使用int32_t存储中间结果至关重要防止乘法累加溢出。这个函数只有两次乘法、一次加法和一次移位在大多数MCU上都能在几个微秒内完成。4.3 参数整定α不是随便设的α的选取直接决定滤波器的截止频率fc。它们的关系是α ≈ 2πfcΔT当fc远小于采样频率时。例如采样率100HzΔT0.01s想要截止频率10Hz那么α ≈ 2*3.14*10*0.01 0.628。这个值很大滤波效果很弱。通常我们会让截止频率远低于采样频率。比如采样率100Hz想滤掉10Hz以上的噪声可以设fc5Hz则α≈0.314。更实用的方法是试凑法。在真实系统或仿真中给一个阶跃信号或叠加了噪声的信号观察不同α下输出的平滑度和响应速度选择一个平衡点。我的经验是对于慢变信号温度α可以小到0.01~0.05对于稍快的信号速度α可以在0.1~0.3之间。5. 二阶IIR滤波器更陡的衰减与实现陷阱当一阶滤波器的性能不能满足要求时比如阻带衰减不够我们就需要二阶甚至更高阶的IIR滤波器。最经典的就是二阶巴特沃斯低通滤波器。5.1 直接I型与直接II型结构一个二阶IIR滤波器的差分方程是y[n] b0*x[n] b1*x[n-1] b2*x[n-2] - a1*y[n-1] - a2*y[n-2]这里有5个系数b0, b1, b2, a1, a2和4个状态变量x[n-1], x[n-2], y[n-1], y[n-2]。实现结构有两种直接I型严格按照差分方程实现需要4个状态变量。结构直观但数值精度可能稍差。直接II型典范型通过转换只需要2个状态变量是更常用、更节省内存的形式。下面我们以实现一个二阶巴特沃斯低通滤波器为例系数可以通过MATLAB、Pythonscipy.signal或在线工具计算得到。5.2 C语言实现直接II型假设我们已经设计好了一个截止频率为fc归一化频率fc 截止频率/采样频率的二阶巴特沃斯低通滤波器并得到了其系数通常是浮点数。我们需要将其转换为定点数。// 假设采样率Fs1000Hz, 截止频率Fc50Hz // 使用工具计算得到的浮点系数 // b00.2066, b10.4132, b20.2066 // a1-0.3695, a20.1958 // 注意差分方程形式为 y[n] b0*x[n] b1*x[n-1] b2*x[n-2] - a1*y[n-1] - a2*y[n-2] // 很多工具输出的是 a0, a1, a2且方程为 a0*y[n] b0*x[n]...我们的a1,a2是这里的负号。 #define Q 14 // 使用Q14定点格式 #define B0_Q14 (int16_t)(0.2066 * (1Q)) // 3382 #define B1_Q14 (int16_t)(0.4132 * (1Q)) // 6767 #define B2_Q14 (int16_t)(0.2066 * (1Q)) // 3382 #define A1_Q14 (int16_t)(-0.3695 * (1Q)) // -6053 #define A2_Q14 (int16_t)(0.1958 * (1Q)) // 3208 int16_t SecondOrderButterworthLPF_DirectII(int16_t x_new) { // 直接II型典范型状态变量 static int32_t w1 0, w2 0; // 注意用32位存储防溢出 int32_t w0; // 中间变量 int32_t y_out_q; // 输出Q14格式 // 典范型计算公式 // w0 x[n] - a1*w1 - a2*w2 // y[n] b0*w0 b1*w1 b2*w2 // 然后更新状态 w2 w1; w1 w0; w0 (int32_t)x_new - ((A1_Q14 * w1) Q) - ((A2_Q14 * w2) Q); y_out_q ((B0_Q14 * w0) Q) ((B1_Q14 * w1) Q) ((B2_Q14 * w2) Q); // 更新状态 w2 w1; w1 w0; return (int16_t)(y_out_q); }严重警告与调试技巧系数转换的符号这是最大的坑不同的工具、不同的传递函数形式如tf或zpk生成的系数其对应的差分方程符号可能不同。务必根据你所用工具的文档明确其输出的系数对应哪个方程。最稳妥的方法是用工具生成系数后在MATLAB或Python里用浮点仿真正确性再移植到C。定点化精度与溢出系数定点化时* (1Q)会引入误差。Q值越大精度越高但乘法越容易溢出。这里状态变量w1, w2和中间结果w0、y_out_q必须使用int32_t。你需要根据输入信号的范围和系数大小估算中间结果的范围确保不会溢出。有时需要牺牲一些Q值降低精度来保证不溢出。稳定性IIR滤波器可能因为系数定点化误差或数值舍入而变得不稳定输出发散。如果发现输出越来越大直至饱和很可能是这个问题。解决方法包括使用更高精度的定点数如Q值更大的int32_t系数、采用更稳定的结构如级联的二阶节、或者在算法中加入饱和处理。6. FIR滤波器实现线性相位的代价与优化当你的应用对相位有严格要求时比如通信中的匹配滤波、图像处理FIR滤波器是更好的选择。一个N阶FIR滤波器就是一个N-1阶的滑动平均的加权版本。6.1 基本原理与直接型实现FIR滤波器的差分方程很简单y[n] b0*x[n] b1*x[n-1] ... b_{N-1}*x[n-(N-1)]它没有反馈所以绝对稳定。实现起来就是一个乘累加操作。#define FIR_ORDER 32 // 滤波器阶数 static int16_t fir_coeffs[FIR_ORDER] { /* ... 你的系数 ... */ }; // Q15格式 static int16_t fir_buffer[FIR_ORDER] {0}; // 循环缓冲区 static int buffer_index 0; int16_t FIR_Filter_Direct(int16_t new_sample) { int32_t acc 0; // 累加器32位 int i, j; // 1. 新数据存入缓冲区 fir_buffer[buffer_index] new_sample; // 2. 乘累加计算 j buffer_index; for (i 0; i FIR_ORDER; i) { acc (int32_t)fir_coeffs[i] * fir_buffer[j]; j--; if (j 0) { j FIR_ORDER - 1; // 循环缓冲区回绕 } } // 3. 更新缓冲区索引 buffer_index; if (buffer_index FIR_ORDER) { buffer_index 0; } // 4. 输出结果 (Q15 * Q15 Q30, 右移15位变回Q15) return (int16_t)(acc 15); }这个实现使用了循环缓冲区避免了每次计算后都要移动大量数据是标准做法。但计算量是O(N)N很大时比如128阶在低速MCU上可能无法实时完成。6.2 针对ARM Cortex-M的优化利用SIMD和DSP库对于ARM Cortex-M4/M7/M33等带有DSP指令集的芯片我们可以利用其单指令多数据和乘累加指令来极大加速FIR运算。CMSIS-DSP库提供了高度优化的函数。#include arm_math.h #define FIR_ORDER 32 static float32_t fir_coeffs[FIR_ORDER] { /* ... 浮点系数 ... */ }; static float32_t fir_state[FIR_ORDER BLOCK_SIZE - 1]; // 状态缓冲区 arm_fir_instance_f32 S; // FIR滤波器实例结构体 void FIR_Filter_Init(void) { arm_fir_init_f32(S, FIR_ORDER, fir_coeffs, fir_state, BLOCK_SIZE); } void FIR_Filter_Process(float32_t *pSrc, float32_t *pDst, uint32_t blockSize) { arm_fir_f32(S, pSrc, pDst, blockSize); }使用CMSIS-DSP库你只需要关心系数设计和初始化。arm_fir_f32函数内部会使用SIMD指令进行并行乘累加速度比纯C实现快一个数量级。如果你的芯片有FPU直接使用浮点系数和运算会更方便。6.3 FIR系数设计窗函数法如何得到fir_coeffs这些系数最常用的方法是窗函数法。你首先确定想要的频率响应如低通、高通、带通然后计算其理想脉冲响应的无限长序列最后用一个有限长的窗函数去截断它。常用的窗有矩形窗、汉宁窗、汉明窗、布莱克曼窗等它们是在主瓣宽度过渡带和旁瓣衰减阻带抑制之间的权衡。你可以用MATLAB的fir1函数或者Python SciPy的scipy.signal.firwin函数来轻松设计。例如设计一个32阶、截止频率0.2倍奈奎斯特频率的低通FIR滤波器汉明窗import scipy.signal as signal numtaps 32 cutoff 0.2 coeffs signal.firwin(numtaps, cutoff, windowhamming)然后将得到的浮点系数数组coeffs定点化后放入你的C代码中。7. 进阶话题滤波器在真实项目中的调优与测试理论实现只是第一步让滤波器在真实系统中稳定、高效地工作才是真正的挑战。7.1 多速率采样与抗混叠如果你的信号频率成分很宽直接以系统最高频率采样然后进行数字滤波计算负担会很大。这时可以采用多速率信号处理先用一个简单的模拟滤波器或数字滤波器进行抗混叠防止高频噪声混叠到低频然后进行降采样再在较低的采样率下进行更复杂的滤波处理。这能大幅降低对处理器性能的要求。在软件上降采样不是简单丢弃数据通常需要配合一个抗混叠的低通滤波器。7.2 动态调整滤波参数自适应滤波在一些场景下信号的特性或噪声的特性是变化的。例如一个运动传感器在静止时需要很强的平滑滤波而在快速运动时需要快速响应。这时就需要动态调整滤波系数。比如可以根据信号的方差或差分值来动态调整一阶低通滤波器的α值。这属于自适应滤波的简单应用实现起来逻辑不复杂但效果显著。7.3 测试与验证时域与频域如何验证你的滤波器工作正常时域测试输入一个阶跃信号比如从0突变到1000观察输出信号的上升时间和超调量。输入一个正弦波观察输出幅度的衰减和相位的延迟。这些可以直观反映滤波器的动态特性。频域测试更专业如果你有条件可以输入一个扫频信号测量系统在不同频率下的增益和相位绘制伯德图。或者采集一段带噪声的真实信号在PC上用MATLAB或Python进行FFT分析对比滤波前后频谱的变化这是最有力的证据。在资源受限的嵌入式设备上我们通常做时域测试。一个实用的技巧是通过串口或DAC将滤波前后的波形实时输出用示波器或上位机软件观察一目了然。7.4 资源与性能的永恒权衡最后嵌入式滤波器的实现永远是在效果、速度、内存、精度之间做权衡。效果 vs 速度/内存高阶滤波器、FIR滤波器效果更好但计算量大、内存占用多。精度 vs 速度浮点数精度高但速度慢定点数速度快但需要精心设计Q值防止溢出和精度损失。通用 vs 专用通用的滤波器函数灵活但可能有冗余计算针对特定传感器写的专用滤波代码往往极其精简高效。我的习惯是在项目初期先用PC工具Python/MATLAB设计并仿真滤波器验证算法有效性。然后在嵌入式端先用浮点数实现原型确保逻辑正确。最后根据MCU的资源和实时性要求进行定点化优化、内存优化甚至汇编级优化。记住没有最好的滤波器只有最适合你当前项目约束的滤波器。
返回列表