IIR滤波器的C实现

本文详细介绍了如何在STM32微控制器上实现并优化IIR滤波器,包括滤波器结构转换、量化处理、结构调整及代码优化等关键步骤,确保滤波器在实际应用中的稳定性和性能。




第一步:
点击菜单中的Edit->Convert Structure 选择Direct Form I ,SOS,(必须是Direct Form I,  II不行)
一般情况下,按照默认设置,fdatool设计都是由二阶部分串联组成的。
这种结构的滤波器稳定性比一个section的要好很多,其他方面的性能也好些。
如果不是的话,点击Convert to second order sections。
这时,滤波器的结构(structure)应该显示为 Direct Form I,second order sections
第二步:
选择quantize filter,精度选择single precision floating point (单精度浮点)
之所以不用定点是因为噪声太大,也不容易稳定。
点击菜单中的Targets -> generate c header ,选择export as:single precision floating point (单精度浮点)
填写变量名称时,把NUM改成IIR_B,DEN改成IIR_A,其他不用动,保存为iir_coefs.h
保存好的文件如下:
//一大堆注释
//然后:
/* General type conversion for MATLAB generated C-code  */
#include "tmwtypes.h"
/* 
* Expected path to tmwtypes.h 
* C:\Program Files\MATLAB\R2010a\extern\include\tmwtypes.h 
*/
/*
* Warning - Filter coefficients were truncated to fit specified data type.  
*   The resulting response may not match generated theoretical response.
*   Use the Filter Design & Analysis Tool to design accurate
*   single-precision filter coefficients.
*/
#define MWSPT_NSEC 9
const int NL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1 };
const real32_T IIR_B[MWSPT_NSEC][3] = {
{
     0.8641357422,              0,              0 
  },
  {
                1,             -2,              1 
  },
  {
     0.9949035645,              0,              0 
  },
  {
                1,   -1.999938965,              1 
  },
  {
     0.9985351563,              0,              0 
  },
  {
                1,    -1.99987793,              1 
  },
{
     0.9996337891,              0,              0 
  },
  {
               1,    -1.99987793,              1 
  },
  {
                1,              0,              0 
  }
};
const int DL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1 };
const real32_T IIR_A[MWSPT_NSEC][3] = {
  {
                1,              0,              0 
  },
  {
                1,   -1.938049316,   0.9401855469 
  },
  {
                1,              0,              0 
  },
  {
                1,   -1.989501953,   0.9900512695 
  },
{
                1,              0,              0 
  },
  {
                1,   -1.996887207,   0.9971923828 
  },
  {
                1,              0,              0 
  },
  {
                1,   -1.999084473,   0.9993286133 
  },
  {
                1,              0,              0 
  }
};
第三步:
打开iir_coefs.h把MWSPT_NSEC替换成IIR_NSEC, 
NL、DL数组删除掉,real32_T改成float ,
其中有一个#include "twmtypes.h",不要它了,删掉
改完的文件如下:
#define IIR_NSEC 9
  //原来叫做MWSPT_NSEC
const float IIR_B[IIR_NSEC][3] = {
  //为什么改为float很明显了吧
  {
    0.8641357422,              0,              0
  },
  {
                1,            -2,              1
  },
  {
    0.9949035645,              0,              0
  },
  {
                1,  -1.999938965,              1
  },
{
    0.9985351563,              0,              0
  },
  {
                1,    -1.99987793,              1
  },
  {
   0.9996337891,              0,              0
  },
  {
                1,    -1.99987793,              1
  },
  {
                1,              0,              0
  }
};
const float IIR_A[IIR_NSEC][3] = {
  {
                1,              0,              0
  },
  {
                1,  -1.938049316,  0.9401855469
  },
  {
                1,              0,              0
  },
  {
                1,  -1.989501953,  0.9900512695
  },
  {
                1,              0,              0
  },
  {
                1,  -1.996887207,  0.9971923828
  },
  {
                1,              0,              0
  },
  {
                1,  -1.999084473,  0.9993286133
  },
  {
                1,              0,              0
  }
};
保存文件,然后使用以下代码进行滤波
这段代码是根据Direct Form I 2阶IIR滤波的差分方程编写的
a0*y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] -a2*y[n-2];
//iir_filter.c
#include "datatype.h"
#include "iir_filter.h"
#include "iir_coefs.h"
static float y[IIR_NSEC][3];
static float x[IIR_NSEC+1][3];
int16 iir_filter(int16 in)
{
      uint16 i;
      x[0][0] = in;
      for(i=0;i<IIR_NSEC;i++)
      {
          y[0] =x[0]*IIR_B[0]+x[1]*IIR_B[1]+x[2]*IIR_B[2]-y[1]*IIR_A[1]-y[2]*IIR_A[2];
          y[0] /= IIR_A[0];
       
          y[2]=y[1];y[1]=y[0];
          x[2]=x[1];x[1]=x[0];
      
          x[i+1][0] = y[0];
      }
      if( x[IIR_NSEC][0]>32767)  x[IIR_NSEC][0]=32767;
      if( x[IIR_NSEC][0]<-32768) x[IIR_NSEC][0]=-32768;
      return  ((int16)x[IIR_NSEC][0]);   
   
}
//复位滤波器
void iir_reset(void)
{
    uint16 i,j;
    for(i=0;i<IIR_NSEC+1;i++)
    {
      for(j=0;j<3;j++)
      {
          x[j]=0;
      }
    }
  
    for(i=0;i<IIR_NSEC;i++)
    {
      for(j=0;j<3;j++)
      {
          y[j]=0;
      }
    }
}
//iir_filter.h
#ifndef _IIR_FILTER_H__
#define _IIR_FILTER_H__
int16 iir_filter(int16 x);
void iir_reset(void);
#endif
使用方法:
首先写好iir_coefs.h,然后调用iir_filter.c对数据流进行滤波
一个伪代码例子:
while(运行中)
{
保存到SD卡(iir_filter(读取ADC采样值()));
}
这个函数比STM32 DSP库中的函数要好很多,DSP库中的2个IIR滤波函数都不能连续处理数据流。
记得在开始滤波之前重置滤波器
iir_reset();

<think>好的,我现在需要帮助用户了解IIR滤波器的C语言实现。首先,我得回忆一下IIR滤波器的基本概念和原理。IIR是无限脉冲响应滤波器,和FIR不同,它使用递归结构,所以具有更高的效率,但可能存在稳定性问题。用户可能已经知道这些基础,但可能需要具体的实现步骤。 接下来,用户需要的是C语言实现方法。我需要考虑如何将差分方程转化为代码。通常,IIR滤波器的差分方程是y[n] = Σa_k*y[n-k] + Σb_m*x[n-m]。这里,系数a和b是关键,需要提前计算或设计好。 然后,我应该解释直接I型或直接II型结构。直接II型更节省内存,因为它减少了延迟单元的数量。可能需要在回答中提到结构的选择,比如使用二阶分段(Biquad)来避免高阶滤波器的不稳定性。 接下来,关于C语言的具体实现,需要定义状态变量数组来保存之前的输入和输出。比如,对于二阶IIR,可能需要两个过去输入和输出的值。然后,编写一个函数,每次处理一个样本,更新状态并计算输出。 还需要注意系数的存储方式,通常使用数组来保存a和b的系数。归一化的问题也要考虑,比如是否将a0归一化为1,这样在计算时可以节省一次乘法。 稳定性问题可能也需要提及,比如在设计滤波器时确保极点位于单位圆内,但在实现时可能还需要量化或缩放系数来防止溢出。 另外,用户可能关心如何测试这个滤波器,比如使用正弦波输入,观察输出是否符合预期,或者用阶跃响应测试。可能还要提到定点数和浮点数的选择,取决于处理器的能力。 最后,示例代码是关键,提供一个简单的二阶IIR滤波器实现,展示如何初始化状态,处理每个样本。同时,提醒用户注意系数设计和结构选择的重要性,以及可能的优化方法,如循环展开或使用移位代替除法。 检查是否有遗漏的部分,比如是否解释了差分方程到代码的转换,是否提到了递归计算和状态保存,以及如何处理实时采样数据。可能需要分步骤解释,先理论后代码,结构清晰。 确保语言通俗易懂,避免过多数学公式,但必要的公式用LaTeX标注。例如,传递函数和差分方程的表示需要正确使用$$格式。同时,注意用户可能对某些术语不熟悉,比如Biquad,需要简要解释。 总结起来,回答的结构应该是:介绍IIR滤波器特点,实现步骤(差分方程、结构选择、状态变量、代码示例),注意事项(稳定性、定点处理),示例代码,以及可能的测试方法。这样用户能逐步理解并实现自己的IIR滤波器。</think>IIR滤波器(无限冲激响应滤波器)的C语言实现可分为以下关键步骤,采用直接Ⅱ型结构(Biquad结构)实现效率最高: $$H(z) = \frac{\sum_{k=0}^{M}b_kz^{-k}}{1 - \sum_{k=1}^{N}a_kz^{-k}}$$ **一、实现步骤** 1. **差分方程转换** 对N阶IIR滤波器,采用二阶分段实现(Biquad): $$y[n] = b_0x[n] + b_1x[n-1] + b_2x[n-2] - a_1y[n-1] - a_2y[n-2]$$ 2. **状态变量定义** ```c typedef struct { float b0, b1, b2; // 分子系数 float a1, a2; // 分母系数(已假设a0=1) float x1, x2; // 输入延迟单元 float y1, y2; // 输出延迟单元 } Biquad; ``` 3. **单次滤波函数** ```c float biquad_process(Biquad *filter, float input) { // 计算当前输出 float output = filter->b0 * input + filter->b1 * filter->x1 + filter->b2 * filter->x2 - filter->a1 * filter->y1 - filter->a2 * filter->y2; // 更新状态变量 filter->x2 = filter->x1; filter->x1 = input; filter->y2 = filter->y1; filter->y1 = output; return output; } ``` **二、完整示例(二阶低通滤波器)** ```c #include <stdio.h> typedef struct { float b0, b1, b2; float a1, a2; float x1, x2; float y1, y2; } Biquad; void biquad_init(Biquad *f, float b0, float b1, float b2, float a1, float a2) { f->b0 = b0; f->b1 = b1; f->b2 = b2; f->a1 = a1; f->a2 = a2; f->x1 = f->x2 = 0; f->y1 = f->y2 = 0; } float biquad_process(Biquad *f, float in) { float out = f->b0*in + f->b1*f->x1 + f->b2*f->x2 - f->a1*f->y1 - f->a2*f->y2; f->x2 = f->x1; f->x1 = in; f->y2 = f->y1; f->y1 = out; return out; } int main() { // 示例:截止频率1kHz,采样率44.1kHz Biquad lp; biquad_init(&lp, 0.02008, 0.04017, 0.02008, // b系数 -1.561, 0.6414 // a系数 ); // 模拟输入信号(示例:直流+噪声) float input[] = {1,1,1,1,1,1,1,1,1,1}; for(int i=0; i<10; i++) { printf("Input: %.2f -> Output: %.4f\n", input[i], biquad_process(&lp, input[i])); } return 0; } ``` **三、关键注意事项** 1. **系数归一化**:设计滤波器时需保证a0=1,系数需预先归一化 2. **稳定性校验**:通过`a1² + a2² < 1`验证极点位置(简化判据) 3. **量化误差**:定点实现时需进行Q格式转换 ```c // 定点数示例(Q15格式) typedef struct { int16_t b0, b1, b2; int16_t a1, a2; int32_t x1, x2; // 扩展精度 int32_t y1, y2; } Biquad_Q15; ``` 4. **级联结构**:高阶滤波器建议采用多个二阶节级联 5. **内存优化**:可通过循环缓冲实现延迟线 **四、性能优化技巧** - 使用ARM CMSIS-DSP库的`arm_biquad_cascade_df1_f32`函数(ARM Cortex-M) - 采用SIMD指令并行计算多个样本 - 对于实时处理,建议采用环形缓冲区 - 使用查表法加速浮点运算(低端MCU) 实际工程中建议先用MATLAB/fdatool生成滤波器系数,再通过自动代码生成工具导出C代码结构。
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值