何为IIR滤波器?无限冲激响应(IIR:Infinite Impulse Response)是一种适用于许多线性时不变系统的属性,这些系统的特征是具有一个冲激响应h(t),该冲激响应h(t)不会在特定点上完全变为零,而是无限期地持续。这与有限冲激响应(FIR:Finite Impulse Response)系统形成对比,在有限冲激响应(FIR)系统中,对于某个有限T,在时间t> T时,冲激响应确实恰好变为零。线性时不变系统的常见示例是大多数电子和数字滤波器。具有此属性的系统称为IIR系统或IIR滤波器。对于什么叫冲激响应,这里就不展开解释了,有兴趣的可以查阅相关书籍。 这是常见的教科书式数学严谨定义,很多人看到这一下就蒙了,能说人话吗? 线性时不变系统理论俗称LTI系统理论,源自应用数学,直接在核磁共振频谱学、地震学、电路、信号处理和控制理论等技术领域运用。它研究的是线性、非时变系统对任意输入信号的响应。虽然这些系统的轨迹通常会随时间变化(例如声学波形)来测量和跟踪,但是应用到图像处理和场论时,LTI系统在空间维度上也有轨迹。因此,这些系统也被称为线性非时变平移,在最一般的范围理论给出此理论。在离散(即采样)系统中对应的术语是线性非时变平移系统。由电阻、电容、电感组成的电路是LTI系统的一个很好的例子。比如一个运放系统在一定频带范围内满足信号的时域叠加,输入一个100Hz和200Hz正弦信号,输出频率是这两种信号的线性叠加。 用数学对LTI系统描述: 线性:输入x1(t),产生响应 y1(t),而输入x2(t),产生相应y2(t) , 那么放缩和加和输入 ax1(t)+bx1(t), 产生放缩、加和的响应ay1(t)+by1(t),其中a和b是标量,对于任意的有: 输入 产生响应为: 时不变性:指如果将系统的输入信号延迟δ秒,那么得到的输出响应也相应延时δ秒。用数学描述,也即如果输入x1(t),产生响应y1(t) ,而输入x1(t+δ) ,产生响应 y1(t+δ)。
这么描述还是不易懂,来个图,有图有真相: 假定一个信号放大电路对100Hz正弦信号放大2倍,则输出为:
上面这么多文字只是为了描述在什么场合可以使用IIR滤波器对信号进行数字滤波。总结而言,就是在线性时不变系统中适用。换言之,在大多数电路系统中我们都可以尝试采用IIR滤波器进行数字滤波。 那么究竟什么是IIR滤波器呢?从数字信号处理的书籍中我们能看到这样的Z变换信号流图: Z的-1次方表示延迟一拍,在数字系统中表示对于输入信号而言,即为上一次采样值,对于输出而言,即为上一次的输出值。 在时域中对于上述流图,用时域描述即为:
上述数字滤波器,如果从编程的角度来看,x(n-1),表示上一次的信号,可能是来自ADC的上次采样,而y(n-1)则为上一次滤波器的输出值,对应就比较好理解x(n-N)就表示前第n次输入样本信号,而y(n-M)则为前第M次滤波器的输出。 说了这么多,只是为了更好的理解概念,只有概念理解正确,才能使用正确。概念理解这对工程师而言,非常之重要。 如何设计呢?MATLAB提供了非常容易使用的FDATool帮助我们设计数字滤波器,真正精彩的地方开始了,让我们拭目以待究竟如何一步一步设计并实施一个IIR滤波器。首先打开MATLAB,在命令行中敲fdatool,然后敲回车 弹出窗体就是fdatool了,如下: A(V)(dB)=20lg(Vo/Vi);电压增益,Vo 为输出电压,Vi为输入电压 A(I)(dB)=20lg(Io/Ii);电流增益,Io 为输出电流,Ii为输入电流 A(p)(dB)=10lg(Po/Pi);功率增益,Po 为输出功率,Pi为输入功率
滤波器类型:这里有Butterworth(巴特沃斯)、Chebyshev Type I,Chebyshev Type II、(切比雪夫)、Elipic 等可选。 巴特沃斯 Butterworth,也被称作最大平坦滤波器。巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦,没有纹波。 切比雪夫 Chebyshev,是在通带或阻带上频率响应幅度等波纹波动的滤波器。切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快,但频率响应的幅频特性不如后者平坦。 椭圆 Elliptic,椭圆滤波器是在通带和阻带等波纹的一种滤波器。 …这里就不一一介绍了,有兴趣可以去查信号处理书籍。
就其特点,这里对其中几种略作介绍: 假设我们需要设计一个IIR滤波器,采样率为32000Hz, 有用信号频率在10000Hz内,设计IIR滤波器对信号进行数字滤波。这里为节省算力,我们指定滤波器的阶数,也即传递函数中N/M中的最大值,一般而言N大于M。 这里指定阶数为8阶,类型指定为巴特沃斯型IIR滤波器,输入阶数8阶,采样率32000Hz,然后点击Design Filter如下图所示: 其相频响应曲线如下: 除此之外,我们还可以将幅频与相频曲线放在一个频率坐标上去看设计结果: 导出滤波器参数,这里我们选择,
然后就得到了一个文件,保存2KHz_LPF.fcf,文件名随你喜欢。 文件内容如下: % Generated by MATLAB(R) 8.4 and the Signal Processing Toolbox 6.22. % Generated on: 27-Mar-2020 21:27:06
% Coefficient Format: Decimal
% Discrete-Time IIR Filter (real) % ------------------------------- % Filter Structure : Direct-Form II, Second-Order Sections % Number of Sections : 4 % Stable : Yes % Linear Phase : No
SOS Matrix: 1 2 1 1 -1.7193929141691948 0.8610574795347461 1 2 1 1 -1.5237898734101736 0.64933827386370635 1 2 1 1 -1.4017399331200424 0.51723237044751591 1 2 1 1 -1.3435020629061745 0.45419615396638446
Scale Values: 0.035416141341387819 0.031387100113383172 0.028873109331868367 0.027673522765052503
至此设计工作就结束了,马上进入滤波器的部署测试阶段。 这里有个概念需要略作解释:什么叫直接II型 SOS 所谓直接II型,SOS(second order section)理解很简单,本质是将IIR Z传递函数分解为上述二阶块的级联形式。
部署测试滤波器到这里,没有经验的朋友可能会说,这么一堆参数我该咋用呢? 需要自己去写前面描述的计算公式吗?当然你也可以这么做,这里就不写了,ARM的CMSIS库已经帮大家设计好了种类繁多的数字信号处理函数实现了,而且经过了测试,这里直接拿来用即可。有兴趣自己写也不难,只要理解Z传递函数概念内涵,非常容易实现。这里我们采用32位浮点实现函数: arm_biquad_cascade_df1_f32。该函数位于: CMSIS\DSP\Source\FilteringFunctions\arm_biquad_cascade_df1_init_f32.c CMSIS\DSP\Source\FilteringFunctions\arm_biquad_cascade_df1_f32.c 我们来看一看这个函数: arm_biquad_cascade_df1_init_f32.c: /* *作用 :初始化滤波器 *S :指向浮点SOS级联结构的实例。 *numStages:滤波器中二阶SOS的数量 *pCoeffs :滤波器参数指针,参数按下列顺序存储 * {b10, b11, b12, a11, a12, b20, b21, b22, a21, a22, ...} *pState :历史状态缓冲区指针 */ void arm_biquad_cascade_df1_init_f32( arm_biquad_casd_df1_inst_f32 * S, uint8_t numStages, const float32_t * pCoeffs, float32_t * pState) { /* Assign filter stages */ S->numStages = numStages;
/* Assign coefficient pointer */ S->pCoeffs = pCoeffs;
/* Clear state buffer and size is always 4 * numStages */ memset(pState, 0, (4U * (uint32_t) numStages) * sizeof(float32_t));
/* Assign state pointer */ S->pState = pState; }
arm_math.h 定义了须用到的结构体,对于本例相关的结构体为arm_biquad_casd_df1_inst_f32 typedef struct { unsigned int numStages; /*2阶节的个数,应为2*numStages. */ float *pState; /*状态系数数组指针,数组长度为4*numStages*/ float *pCoeffs; /*系数数组指针, 数组的长度为5*numStages.*/ } arm_biquad_casd_df1_inst_f32;
滤波器具体滤波函数为: arm_biquad_cascade_df1_f32 /** * *S :指向浮点Biquad级联结构的实例. * *pSrc :指向输入数据块。 * *pDst :指向输出数据块。 * blockSize:每次调用要处理的样本数。 * 返回值 :无. */ void arm_biquad_cascade_df1_f32( const arm_biquad_casd_df1_inst_f32 * S, float * pSrc, float * pDst, unsigned int blockSize) { float *pIn = pSrc; /*源指针 */ float *pOut = pDst; /*目的指针 */ float *pState = S->pState; /*状态指针 */ float *pCoeffs = S->pCoeffs; /*参数指针 */ float acc; /*累加器 */ float b0, b1, b2, a1, a2; /*滤波器参数 */ float Xn1, Xn2, Yn1, Yn2; /*滤波器状态变量*/ float Xn; /*临时输入 */ unsigned int sample, stage = S->numStages; /*循环计数 */
do { /* Reading the coefficients */ b0 = *pCoeffs++; b1 = *pCoeffs++; b2 = *pCoeffs++; a1 = *pCoeffs++; a2 = *pCoeffs++;
Xn1 = pState[0]; Xn2 = pState[1]; Yn1 = pState[2]; Yn2 = pState[3];
sample = blockSize >> 2u;
while(sample > 0u) { /* 读第一个输入 */ Xn = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */ Yn2 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);
/* Store the result in the accumulator in the destination buffer. */ *pOut++ = Yn2;
/* 每次计算输出后,状态都应更新. */ /* 状态应更新为: */ /* Xn2 = Xn1 */ /* Xn1 = Xn */ /* Yn2 = Yn1 */ /* Yn1 = acc */
/* Read the second input */ Xn2 = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */ Yn1 = (b0 * Xn2) + (b1 * Xn) + (b2 * Xn1) + (a1 * Yn2) + (a2 * Yn1);
/* 将结果存储在目标缓冲区的累加器中. */ *pOut++ = Yn1;
/* 每次计算输出后,状态都应更新. */ /* 状态应更新为: */ /* Xn2 = Xn1 */ /* Xn1 = Xn */ /* Yn2 = Yn1 */ /* Yn1 = acc */
/*读第三个输入 */ Xn1 = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */ Yn2 = (b0 * Xn1) + (b1 * Xn2) + (b2 * Xn) + (a1 * Yn1) + (a2 * Yn2);
/* 将结果存储在目标缓冲区的累加器中. */ *pOut++ = Yn2;
/* 每次计算输出后,状态都应更新. */ /* 状态应更新为: */ /* Xn2 = Xn1 */ /* Xn1 = Xn */ /* Yn2 = Yn1 */ /* Yn1 = acc */ /* 读第四个输入 */ Xn = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */ Yn1 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn2) + (a2 * Yn1);
/* 将结果存储在目标缓冲区的累加器中. */ *pOut++ = Yn1;
/* 每次计算输出后,状态都应更新. */ /* 状态应更新为: */ /* Xn2 = Xn1 */ /* Xn1 = Xn */ /* Yn2 = Yn1 */ /* Yn1 = acc */ Xn2 = Xn1; Xn1 = Xn;
/* 递减循环计数器 */ sample--; }
/* 如果blockSize不是4的倍数, *请在此处计算任何剩余的输出样本。 *不使用循环展开. */ sample = blockSize & 0x3u;
while(sample > 0u) { /* 读取输入 */ Xn = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */ acc = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);
/* 将结果存储在目标缓冲区的累加器中. */ *pOut++ = acc;
/* 每次计算输出后,状态都应更新。 */ /* 状态应更新为: */ /* Xn2 = Xn1 */ /* Xn1 = Xn */ /* Yn2 = Yn1 */ /* Yn1 = acc */ Xn2 = Xn1; Xn1 = Xn; Yn2 = Yn1; Yn1 = acc;
/* d递减循环计数器 */ sample--; }
/* 将更新后的状态变量存储回pState数组中 */ *pState++ = Xn1; *pState++ = Xn2; *pState++ = Yn1; *pState++ = Yn2;
/*第一阶段从输入缓冲区到输出缓冲区. */ /*随后的numStages在输出缓冲区中就地发生*/ pIn = pDst;
/* 重置输出指针 */ pOut = pDst;
/* 递减循环计数器 */ stage--;
} while(stage > 0u); }
开始测试: #include <stdio.h>#include <math.h>/* SOS Matrix: 1 2 1 1 -1.7193929141691948 0.8610574795347461 1 2 1 1 -1.5237898734101736 0.64933827386370635 1 2 1 1 -1.4017399331200424 0.51723237044751591 1 2 1 1 -1.3435020629061745 0.45419615396638446
Scale Values: 0.035416141341387819 0.031387100113383172 0.028873109331868367 0.027673522765052503 做如下转换: 1.缩放 [1 2 1] * 0.035416141341387819 [1 2 1] * 0.031387100113383172 [1 2 1] * 0.028873109331868367 [1 2 1] * 0.027673522765052503 得到: [0.035416141341387819 2*0.035416141341387819 0.035416141341387819] [0.031387100113383172 2*0.031387100113383172 0.031387100113383172] [0.028873109331868367 2*0.028873109331868367 0.028873109331868367] [0.027673522765052503 2*0.027673522765052503 0.027673522765052503] 2.舍掉第四列参数 3.将后两列分别乘以-1,即: 0.035416141341387819 2*0.035416141341387819 0.035416141341387819 -1.7193929141691948 0.8610574795347461 0.031387100113383172 2*0.031387100113383172 0.031387100113383172 -1.5237898734101736 0.64933827386370635 0.028873109331868367 2*0.028873109331868367 0.028873109331868367 -1.4017399331200424 0.51723237044751591 0.027673522765052503 2*0.027673522765052503 0.027673522765052503 -1.3435020629061745 0.45419615396638446 这样就得到了滤波器系数组了 */#define IIR_SECTION 4 /*见前面设计输出为4个SOS块*/static float iir_state[4*IIR_SECTION];/*历史状态缓冲区 */const float iir_coeffs[5*IIR_SECTION]={ 0.035416141341387819,2*0.035416141341387819,0.035416141341387819,1.7193929141691948,-0.8610574795347461, 0.031387100113383172,2*0.031387100113383172,0.031387100113383172,1.5237898734101736,-0.64933827386370635, 0.028873109331868367,2*0.028873109331868367,0.028873109331868367,1.4017399331200424,-0.51723237044751591, 0.027673522765052503,2*0.027673522765052503,0.027673522765052503,1.3435020629061745,-0.45419615396638446 }; static arm_biquad_casd_df1_inst_f32 S; /*假定采样512个点*/#define BUF_SIZE 512#define PI 3.1415926#define SAMPLE_RATE 32000 /*32000Hz*/int main() { float raw[BUF_SIZE]; float raw_4k[BUF_SIZE]; float raw_out[BUF_SIZE];
float raw_noise[BUF_SIZE]; float raw_noise_out[BUF_SIZE];
arm_biquad_casd_df1_inst_f32 S; FILE *pFile=fopen("./simulation.csv","wt+"); if(pFile==NULL) { printf("file opened failed"); return -1; }
for(int i=0;i<BUF_SIZE;i++) { /*模拟800Hz正弦幅度171,叠加幅度50随机噪声 */ raw[i] = 0.5*1024.0/3*sin(2*PI*800*i/32000.0f)+rand()%50; raw_4k[i] = 0.5*1024.0/3*sin(2*PI*4000*i/32000.0f); /*模拟800Hz +4000Hz+随机噪声叠加输入 */ raw_noise[i] = raw[i] + raw_4k[i]; } /*初始化滤波器,以及滤波*/ arm_biquad_cascade_df1_init_f32(&S, IIR_SECTION, (float *)&iir_coeffs[0], (float *)&iir_state[0]); arm_biquad_cascade_df1_f32(&S, raw, raw_out, BUF_SIZE);
for(int i=0;i<BUF_SIZE;i++) { fprintf(pFile,"%f,",raw[i]); }
fprintf(pFile,"\n"); for(int i=0;i<BUF_SIZE;i++) { fprintf(pFile,"%f,",raw_4k[i]); } fprintf(pFile,"\n"); for(int i=0;i<BUF_SIZE;i++) { fprintf(pFile,"%f,",raw_out[i]); }
/*初始化滤波器,以及滤波*/ arm_biquad_cascade_df1_init_f32(&S, IIR_SECTION, (float *)&iir_coeffs[0], (float *)&iir_state[0]); arm_biquad_cascade_df1_f32(&S, raw_noise, raw_noise_out, BUF_SIZE);
fprintf(pFile,"\n"); for(int i=0;i<BUF_SIZE;i++) { fprintf(pFile,"%f,",raw_noise[i]); }
fprintf(pFile,"\n"); for(int i=0;i<BUF_SIZE;i++) { fprintf(pFile,"%f,",raw_noise_out[i]); }
fclose(pFile);
return 0; }
利用csv文件,将模拟数据存储,直接用excel打开,将行数据生成曲线图如下:
有兴趣也可以写个界面直接显示,甚至绘制出谱线图,做进一步分析。 第一幅图,为800Hz信号混入随机噪声的波形 第二幅图,为4000Hz信号,对假定系统为无用干扰信号
第三幅图, 为800Hz 混入随机噪声过滤后,已经很好的还原有用信号频率 第四幅图, 为800Hz信号混入随机噪声,同时叠加4000Hz干扰的波形,对系统而言,从时域中,明显可见,有用信号已经完全扭曲 第五幅图,为800Hz信号混入随机噪声,同时叠加4000Hz干扰的输入,经过该低通滤波器后的波形,与第三幅图基本一样,已经非常好的滤除了干扰信号。
总结: IIR滤波器在线性时不变系统中可以很好的解决工程中一般噪声问题 如果需要设计带通、高通滤波器其步骤基本类似,只是滤波器的参数以及SOS块个数可能不一样而已 需要提醒的时,IIR的相频响应不线性,如果系统对相频响应有严格要求,就需要采用其他的数字滤波器拓扑形式了 实际应用中,如果阶数不高时,现在算力强劲的单片机或者DSP以及可以直接使用浮点处理。 如果对处理速度有严格的实时要求,需要在极短时间进行滤波处理,可以考虑降低阶数,或采用定点IIR滤波算法实现。也或者将文中函数进行汇编级优化。
~END~
|