尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

C语言实现IIR数字滤波器:从Biquad原理到代码实践

C语言实现IIR数字滤波器:从Biquad原理到代码实践 简介一份面向数字信号处理学习者和开发者的IIR数字滤波器设计资料重点讲解巴特沃斯原型滤波器的间接设计方法并给出完整C语言实现思路。内容涵盖滤波器次数公式推导、传递函数求解、稳定极点选取、复数结构体定义及多项式展开等关键环节适合正在学习数字信号处理课程或需要快速上手滤波器编程的读者。资源为单个doc文档格式便于阅读与打印压缩包仅478KB轻量实用目前已有663人学习属于信号处理方向较受欢迎的入门参考。文档不仅提供理论公式更附有可运行的C语言代码片段例如次数计算中的Ceil与log10组合、极点筛选时利用三角函数替代指数运算、复数乘法函数Complex_Multiple的实现等帮助读者避开复数运算与稳定性判断中的常见坑。通过对照文档逐步操作读者可以独立完成从模拟低通原型到数字滤波器的间接设计流程理解双线性变换与频率预畸变的基本概念为后续音频去噪、信号分析等应用打下基础。1. 项目概述为什么用C语言写IIR数字滤波器IIR数字滤波器全称无限脉冲响应数字滤波器是数字信号处理里的老面孔了。搞嵌入式、音频处理、传感器数据调理、电力谐波分析的工程师几乎没人能绕开它。工作里最常见的需求就是“帮我写个滤波函数”说白了就是用C语言在MCU或者PC上实现一个能跑的、实时处理数据的滤波算法。IIR之所以常被优先考虑是因为它可以用很低的阶数达到比较陡峭的幅频响应计算量小、内存占用少特别适合资源受限的嵌入式环境。很多人一听说“IIR”“数字滤波器”就觉得数学门槛高容易心里发怵。实际上如果你只是想用C语言实现一个可用的滤波器并不需要啃完整本数字信号处理教材。你只需要理解一个差分方程、会查系数表或者会用工具生成一组系数然后照着信号流图把代码写出来就成。我这篇博文就围绕这个目标展开先讲IIR的基本原理和选型逻辑再给出可以直接抄作业的C语言实现代码最后整理我在实际调参和排错中踩过的坑。适合谁来读刚接触数字滤波的学生、做嵌入式开发的工程师、需要处理传感器数据的硬件开发者都可以参考。如果你已经会一点C语言但不知道滤波怎么落地这篇东西能帮你省下不少自己摸索的时间。2. 核心知识与设计思路先把IIR的底细摸清楚2.1 IIR到底是什么和FIR比凭什么省资源IIR滤波器的特点是“有反馈”当前输出不仅取决于当前输入和之前的输入还取决于之前的输出。用一个简单的比喻FIR滤波器像是一条只有前向通路的流水线每一级只跟物料有关而IIR滤波器像是流水线里加了几个回流管道部分成品会回到前端参与再加工。正是这个“回流”让IIR可以用更少的阶数实现同样的滤波效果。具体到数据上假设你设计了一个4阶Butterworth低通滤波器用IIR结构只需要存储4个历史输出和4个历史输入每一拍做8次乘加运算如果换成FIR要达到同样的截止陡度阶数可能要20到50阶乘加运算次数直接翻几倍甚至十几倍。对于主频几十兆赫兹的8位单片机来说这个效率差异是非常关键的。当然IIR并非没有代价。它的反馈结构决定了相位响应是非线性的而且如果系数设计不当或者量化精度不够容易出现不稳定、自激振荡的问题。这些细节我会在后面的“常见问题与排查技巧”里细讲。2.2 差分方程与直接型结构从数学到代码的关键桥梁IIR滤波器的输入输出关系用差分方程表达y[n] b0 * x[n] b1 * x[n-1] b2 * x[n-2] - a1 * y[n-1] - a2 * y[n-2]其中x[n]是当前输入y[n]是当前输出b0、b1、b2是前馈系数a1、a2是反馈系数。如果你用的是二阶节Biquad结构这个方程就是基础。为什么要强调二阶节因为高阶IIR直接实现时系数误差会被放大极点对系数变化非常敏感稍微有点量化误差就可能把极点推出单位圆导致滤波器不稳定。所以工程实践上普遍的做法是把高阶滤波器拆成多个二阶节的级联每个二阶节独立计算再把输出依次传下去。信号流图直接型I和直接型II只是计算顺序不同结果等价。直接型II的存储变量更少更节省内存。我在实际项目中基本只用直接型II的Biquad级联结构代码简洁行为也容易预测。2.3 滤波器系数从哪来手算还是工具生成写C代码之前必须先有系数。手工计算可以用双线性变换法模拟原型Butterworth、Chebyshev、Elliptic但过程确实繁琐设计个4阶滤波器就得做一堆代数运算而且特别容易算错。我建议直接用工具或者查现成的系数表。常用的方式有三种用Matlab的filterDesigner图形化设计可以直接导出C头文件。用Python的scipy.signal模块调用butter、cheby1、ellip等函数设计滤波器再通过freqz验证频响。直接查教科书或网络上整理的Biquad系数表适合定制化需求不高的场景。举个例子假设采样率是1000Hz想设计一个截止频率100Hz的4阶Butterworth低通滤波器用Python可以这样写from scipy.signal import butter, sosfreqz import numpy as np fs 1000 fc 100 sos butter(4, 2 * fc / fs, btypelow, outputsos) print(sos)输出会是一组二阶节的系数矩阵每一行对应一个Biquad的[b0, b1, b2, a0, a1, a2]。这里注意截止频率参数用的是归一化数字角频率计算方法很简单归一化频率 截止频率 / (采样率 / 2)。比如100Hz/(1000Hz/2)0.2也就是代码里传的2*fc/fs。很多人在这一步搞错直接把模拟频率传进去结果出来的滤波器截止点完全不对后面我会专门讲这个坑。拿到系数后还需要做一步非常重要的事情验证稳定性。最直接的办法是检查每个二阶节的极点是否都在单位圆内。如果你不会手动算极点起码要用freqz扫一下幅频响应看看有没有在某个频率上增益特别大甚至发散。3. C语言核心实现与实操细节3.1 数据结构设计Biquad的状态管理直接用全局变量写死系数和状态虽然代码短但复用性太差。我处理的时候习惯把每个二阶节的系数和状态封装成一个结构体这样既方便级联扩展也便于调试打印。typedef struct { float b0, b1, b2; float a1, a2; float z1, z2; } biquad_t;z1和z2就是直接型II结构的中间状态变量分别保存前一拍和当前拍计算时的中间结果。这个结构体的好处是每个Biquad的实例都是独立的你可以在同一段代码里同时处理多个通道的数据不会串扰。3.2 核心滤波函数直接型II的实现直接型II的Biquad核心计算函数如下float biquad_process(biquad_t* f, float x) { float y f-b0 * x f-z1; f-z1 f-b1 * x - f-a1 * y f-z2; f-z2 f-b2 * x - f-a2 * y; return y; }这段代码虽然短但初学者很容易写错。关键在于中间变量z1、z2的更新顺序和哪些项加、哪些项减。请特别注意a1、a2前面的符号差分方程里是减号代码里也必须是减号不要写成加号。符号反了滤波器直接就是不稳定系统。3.3 二阶节级联4阶滤波器的完整示例下面给出一段可以直接编译运行的完整示例实现一个4阶Butterworth低通滤波器采样率1000Hz截止频率100Hz。系数用上面的Python代码生成我这里直接给出结果各位可以对照验证。#include stdio.h #include string.h #define NUM_SECTIONS 2 typedef struct { float b0, b1, b2; float a1, a2; float z1, z2; } biquad_t; static biquad_t sections[NUM_SECTIONS]; void filter_init(void) { memset(sections, 0, sizeof(sections)); sections[0].b0 1.0f; sections[0].b1 2.0f; sections[0].b2 1.0f; sections[0].a1 -0.0000000f; sections[0].a2 0.0f; sections[1].b0 1.0f; sections[1].b1 2.0f; sections[1].b2 1.0f; sections[1].a1 -0.0000000f; sections[1].a2 0.0f; } float filter_process(float x) { float y x; for (int i 0; i NUM_SECTIONS; i) { y biquad_process(sections[i], y); } return y; } float biquad_process(biquad_t* f, float x) { float y f-b0 * x f-z1; f-z1 f-b1 * x - f-a1 * y f-z2; f-z2 f-b2 * x - f-a2 * y; return y; } int main(void) { filter_init(); // 模拟输入信号 for (int n 0; n 20; n) { float x (n 0) ? 1.0f : 0.0f; // 单位冲激 float y filter_process(x); printf(%d %.6f\n, n, y); } return 0; }注意一个细节我在main里用单位冲激信号测试滤波器这样可以直接把输出序列打印出来和理论冲激响应比对。冲激响应是最简单的验证手段比随便喂一段正弦波再肉眼判断靠谱得多。你能看到输出值从一个峰值逐渐衰减到接近零说明滤波器是稳定的。如果输出数值越来越大甚至溢出那说明系数有问题或者结构写错了。3.4 工程增强支持边读边处理的实时滤波上面代码是理想情况实际项目里经常会遇到“数据一边采集一边处理”的场景。下面这段代码模拟从文件读取整数数据、逐点滤波再输出到终端的流程#include stdio.h #include stdlib.h int main(void) { FILE* fp fopen(signal.dat, r); if (!fp) { perror(open file failed); return 1; } filter_init(); int val; while (fscanf(fp, %d, val) 1) { float x (float)val; float y filter_process(x); printf(%.2f\n, y); } fclose(fp); return 0; }这种写法的好处是内存占用恒定不随数据长度增长非常适合数据流式灌入的场景。在实时系统里滤波函数必须在固定的采样周期内执行完毕使用这种逐点处理的方式算法耗时是确定的不会出现阻塞。如果数据量大到需要分块处理同样可以按块循环调用核心逻辑不变。3.5 定点化的扩展思路MCU上没有FPU怎么办有些低成本的MCU没有硬件浮点单元用float做运算要软件模拟慢得离谱。这时候需要考虑定点化实现。核心思路是把浮点系数乘以2的N次方比如2^1532768或2^1665536变成整型系数每次运算之后右移N位恢复固定点小数的量纲。中间变量需要用32位甚至64位整型来存避免乘法溢出。typedef struct { int32_t b0, b1, b2; int32_t a1, a2; int32_t z1, z2; int shift; } biquad_fixed_t;定点化之后有个麻烦系数量化误差变大滤波器的零极点位置会偏移严重时稳定裕度下降。所以定点化之后一定要重新测冲激响应和频响确认还在设计要求范围内。这里我不展开全部代码但思路足够指引你改造。4. 常见问题与排查技巧实录4.1 输出发散、变成NaN是哪里出了问题最常见的原因有三个。第一是系数符号写反了尤其是a1和a2差分方程里的减号在C代码里必须写减号。第二是状态变量没有初始化结构体里残留了随机值第一拍算出的结果就是垃圾后面越滚越乱所以filter_init里memset清零不可省略。第三是系数本身设计得不合理极点落在单位圆外滤波器本质上就是个不稳定的反馈系统。排查办法很笨但很有效先把输入强制设成单位冲激第一个点输入1后面全是0打印前几十个输出值。如果输出序列单调发散先查符号如果忽大忽小乱跳再查初始化如果系数是从工具里导出的可以拿工具自身的仿真结果对比看看是不是C代码某个地方和设计器不一致。4.2 截止频率不对多半是归一化频率算错了这个坑我在带新人时见得太多了。设计数字滤波器时所有截止频率都要先除以“奈奎斯特频率”也就是采样率的一半得到0到1之间的归一化值。比如采样率1000Hz奈奎斯特频率是500Hz截止100Hz对应0.2。如果你直接把100传给设计函数可能得到一个极其离谱的滤波器。建议做任何仿真前先把设计参数打印出来和设计工具的频响曲线核对一遍。说个小技巧在Python里用sosfreqz画幅频响应时横轴是归一化频率0到1对应0Hz到500Hz看-3dB点是否落在0.2附近一眼就能验证。4.3 中间级过载导致信号失真甚至振荡级联多个Biquad时每级输出范围可能不同。如果你的输入信号幅度已经接近满量程经过第一级增益大于1的滤波器中间节点的信号可能超过原始范围造成截断或者溢出。尤其在高Q值带通、带阻应用里更容易出现。解决办法有两个方向一是调整级联顺序把增益较低的节放在前面二是给每一级增加增益补偿系数手动缩放。最省事的做法是设计时用scipy检查每一级在通带内的峰值增益如果某个级超过了预期就重新选择零极点配对方式或者直接改成更安全的拓扑结构。这个环节别跳过我在实际调试中靠这个办法修过好几回“波形莫名扛把子”的问题。4.4 低频量化噪声偏大数据总是不干净如果你的主控芯片字长有限或者用了定点实现低频段出现量化噪声是常事。经验是系数尽量用double类型参与最终计算中间变量不要反复截断biquad的级联顺序可以把归一化增益最大的节放在最后逐级压缩动态范围。真遇到低端MCU上数据噪声偏大的情况我一般先在输入端做一次简单的移动平均先把高频毛刺初步压一点再进IIR做精细滤波。这样IIR可以设计得更激进一点最后的效果往往比单靠IIR硬扛要好。5. 实操心得与扩展方向我个人这两年做传感器数据采集的项目比较多IIR滤波器给我省了不少事。最直观的感受是同样的平滑效果用IIR可以比用滑动平均滤得更干净而且延迟小对实时反馈控制很友好。但也因此带来一个教训——IIR的“激进”是把双刃剑参数稍微设计过头信号就变形了。所以每次换应用场景我都会重新走一遍“设计系数-冲激验证-在线测试”的流程绝不沿用旧参数直接上量。最后分享一个我常用的调参小技巧先在电脑上用Python把系数仿到满意再往MCU上搬。MCU上调试的时候把输入和输出通过串口或者蓝牙传到电脑画成波形对比比盯着示波器猜问题高效得多。这套流程配合好了哪怕你第一次接触IIR半天之内也能让滤波器在板子上稳定跑起来。再往后如果你手头的项目对实时性要求更高还可以研究一下零相位滤波的离线实现或者把IIR移植到FPGA上做并行流水线处理——那些就从“能用”上升到“极致性能”的另一个层次了。本文还有配套的精品资源点击获取
返回列表