怎样使用matlab去计算出IIR滤波器的系数呢

1.简介
在做音频eq的时候，iir滤波器是必不可少的，所以这节主要学习一下iir滤波器，当然，不会偏于理论，而是着重怎么用。
iir滤波器主要是用来处理信号的频率的，这是第一个最基本的认识。
摆一下基本的公式：
系统函数：

差分方程：

再来说一下该滤波器的结构：
a.直接I型
b.直接II型（典范型）
c.级联型
d.并联型
bcd三种类型都是从直接I型转化而来的，这是第二个最基本的认识。
iir滤波器对于叠加了噪声的信号，也就是信号和噪声的频谱相互混叠是无能为力的，这就是其用在eq上一个基本的原因，这是第三个最基本的认识。
具体再使用一个例子来表述如何使用该公式，由于使用matlab仿真条件下，反馈符号是负号，故需要将
代入一个负号。
假设是二阶滤波器，并且分子，分母如下：


样本点如下：

则计算一下输出值：




使用matlab编写上述程序，如下所示：
clc;
clear all;
format long;
b = ［0.6，0.8，1.0］;
a = ［1，0.5，1.5］;
x = ［1，2.5，1.5，0.5］;
y=filter（b，a，x）
有一点要注意，在做dsp开发的时候，a，b的排序并不一定是按照

，有可能是

，这一点需要查看具体的API即可。
2.关于iir滤波器系数
可以使用matlab来计算出iir滤波器的系数，当然，也有其他的工具可以达成这样的功能，例如：一个直接计算的在线工具
3.二阶滤波器例子
3.1 简单的一段2阶滤波器
再来写一点例子程序，例如c代码：
void test_iir_raw（）
{
int i，j，k;
int len = 4;
float in［4］={0.0，0.1，0.2，0.3};
float geq_para［5］={0.1，0.3，0.5，0，0};
float a［3］，b［3］;
float geq_state［2］={0.0};
float geq_state_in［2］={0.0};
float ym;
float g=1.0;
float tmpdata［4］;
a［0］ = 1;
a［1］ = geq_para［1］;
a［2］ = geq_para［2］;
b［0］ = geq_para［0］;
b［1］ = 0;
b［2］ = -b［0］;
for （k = 0; k 《 len; k++）
{
ym = geq_para［0］ * （in［k］ - geq_state_in［1］） - geq_para［2］ * geq_state［1］ - geq_para［1］ * geq_state［0］;
geq_state［1］ = geq_state［0］;
geq_state［0］ = ym;
geq_state_in［1］ = geq_state_in［0］;
geq_state_in［0］ = in［k］;
tmpdata［k］ += ym * g;
}
for（i = 0 ; i 《 4 ; i ++）{
printf（“tmpdata：%fn”，tmpdata）;
}
}
再例如matlab代码：
clc;
clear all;
format long;
b = ［0.1， 0， -0.1］;
a = ［1， 0.3， 0.5］;
xk=［0.0，0.1，0.2，0.3］;
y=filter（b，a，xk）
可以很简单的分析到：
geq_state可以代表y（n）的临时变量.geq_state_in代表的是x（n）的临时变量。
我再给出结果：
tmpdata:0.000000
tmpdata:0.010000
tmpdata:0.017000
tmpdata:0.009900
最后，在说明一下本例子中的一个临时结果：
geq_state［0］=0.009900，geq_state［1］=0.017000
geq_state_in［0］=0.300000，geq_state_in［1］=0.200000
3.2 简单的二段2阶滤波器
这是一个并联的二阶滤波器，相当于一段音频做两次二阶滤波器效果，对应的是在每一个点上加上上一次做的结果值，我列出c代码：
void test_iir_raw（）
{
int i，j，k;
int len = 4;
float in［4］={0.0，0.1，0.2，0.3};
float geq_para［GEQ_BANDS］［5］={{0.1，0.3，0.5，0，0}，{0.1，0.3，0.5，0，0}};
float a［3］，b［3］;
float geq_state［GEQ_BANDS］［2］={0.0};
float geq_state_in［GEQ_BANDS］［2］={0.0};
float ym;
float g=1.0;
float tmpdata［4］;
for（j = 0 ; j 《 GEQ_BANDS ; j ++）{
for （k = 0; k 《 len; k++）
{
ym = geq_para［j］［0］ * （in［k］ - geq_state_in［j］［1］） - geq_para［j］［2］ * geq_state［j］［1］ - geq_para［j］［1］ * geq_state［j］［0］;
geq_state［j］［1］ = geq_state［j］［0］;
geq_state［j］［0］ = ym;
geq_state_in［j］［1］ = geq_state_in［j］［0］;
geq_state_in［j］［0］ = in［k］;
tmpdata［k］ += ym * g;
}
}
for（i = 0 ; i 《 4 ; i ++）{
printf（“tmpdata：%fn”，tmpdata）;
}
}
对应的matlab结果就是上次给出的值中每一项乘以2即可。至此，这个部分结束。
4.定点
首先给出biquad滤波器的差分函数形式：

注意：
Matlab里的计算就是按照上面的式子计算的，但是STM32F4DSP库里的系数a1，a2是取反的。
再贴出来biquad滤波器的实现代码：
for （k = 0; k 《 len; k++）
{
ym = geq_para［j］［0］ * （in［k］ - geq_state_in［j］［1］） - geq_para［j］［2］ * geq_state［j］［1］ - geq_para［j］［1］ * geq_state［j］［0］;
geq_state［j］［1］ = geq_state［j］［0］;
geq_state［j］［0］ = ym;
geq_state_in［j］［1］ = geq_state_in［j］［0］;
geq_state_in［j］［0］ = in［k］;
tmpdata［k］ += ym * g;
}
为了让定点化来的顺利，可以在第一步的时候只定点化ym，也就是说：
geq_state［j］［1］ = geq_state［j］［0］;
geq_state［j］［0］ = ym;
geq_state_in［j］［1］ = geq_state_in［j］［0］;
geq_state_in［j］［0］ = in［k］;
tmpdata［k］ += ym * g;
这部分不变，而先将：
ym = geq_para［j］［0］ * （in［k］ - geq_state_in［j］［1］） -
geq_para［j］［2］ * geq_state［j］［1］ - geq_para［j］［1］ * geq_state［j］［0］;
各个部分进行定点：
long long int tmp1 = geq_para［j］［0］ * （1 《《 26）;
long long int tmp2 = geq_para［j］［2］ * （1 《《 26）;
long long int tmp3 = geq_para［j］［1］ * （1 《《 26）;
tmp1 = （tmp1 * tmp0） 》》 26;
tmp2 = （tmp2 * tmpp1） 》》 26;
tmp3 = （tmp3 * tmpp2） 》》 26;
ym = （tmp1 - tmp2 - tmp3）*1.0 / （1 《《 26）;
之后再将整体浮点数替换，从而真正实现定点化（64位）。为了方便移植到硬件上，还需要进一步定点化。
可以按照如下方法做进一步处理：
long long int tmp0 = （in［k］）*（1 《《 28） - geq_state_in_int［j］［1］;
long long int tmp1 = geq_para［j］［0］ * （1 《《 26）;
long long int tmp2 = geq_para［j］［2］ * （1 《《 26）;
long long int tmp3 = geq_para［j］［1］ * （1 《《 26）;
long long int tmpp1 = geq_state_int［j］［1］ ;
long long int tmpp2 = geq_state_int［j］［0］ ;
tmp1 = （tmp1 * tmp0） 》》 26;
tmp2 = （tmp2 * tmpp1） 》》 26;
tmp3 = （tmp3 * tmpp2） 》》 26;
ym_int = tmp1 - tmp2 - tmp3;
geq_state_int［j］［1］ = geq_state_int［j］［0］;
geq_state_int［j］［0］ = ym_int;
geq_state_in_int［j］［1］ = geq_state_in_int［j］［0］;
geq_state_in_int［j］［0］ = in［k］ * （1 《《 28）;
tmpdata_int［k］ += ym_int * g;