对于fir滤波器，已经在前面的文章中记录了（https://blog.csdn.net/suiji2442/article/details/112394026POWER-Z仿制DIY&关于MATLAB中滤波器设计工具的使用心得记录），其设计和实现都非常简单。如果在嵌入式系统中可以满足且有必要实时iir运算，那么使用iir滤波器相对fir滤波器可以在使用更小的阶数的情况下实现更好的效果。实验证明，可能20阶的iir效果堪比500阶左右的fir滤波器效果。

首先放出iir的matlab仿真代码：

%本程序为直接2型iir滤波器的实现
%当阶数较高或fc与fs相差多个数量级时，建议使用基本二阶结的级联形式，避免计算出来的b和a系数非常小或者非常大，因运算中有限字长导致运算出错

NN=8000;
NNN=20+1;%在这里写iir滤波器的阶数+1


t=0:1/NN:1;%采样1秒的信号
%y= 10+10*cos(2*pi*5*t-pi*30/180)+5*cos(2*pi*50*t-pi*30/180);
%y= 10+10*cos(2*pi*50*t-pi*30/180)+10*cos(2*pi*200*t-pi*30/180);
%y= 10+10*cos(2*pi*5*t-pi*30/180)+10*cos(2*pi*50*t-pi*30/180)+2*cos(2*pi*200*t-pi*30/180);
y= 10+2*cos(2*pi*5*t-pi*30/180)+2*cos(2*pi*30*t-pi*30/180)+10*cos(2*pi*100*t-pi*30/180)+100*cos(2*pi*1000*t-pi*30/180);
N=length(t); %样点个数
figure(1);
plot(t,y);
fs=NN;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率
Y=fft(y(1:N))/N*2;%真实的幅值
figure(2)
plot(f(1:floor(N/2)),abs(Y(1:floor(N/2))));


[b,a]=sos2tf(SOS,G);

pre_x = zeros(1,NNN);% 前几个时刻的输入值
pre_y = zeros(1,NNN);% 前几个时刻的输出值
res_iir = zeros(1,N);

for i=1:N
    pre_x(1) = y(i);
    % 差分方程
    % res_iir(i) = b(1)*pre_x(1)+b(2)*pre_x(2)+b(3)*pre_x(3)+b(4)*pre_x(4)-a(2)pre_y(2)-a(3)pre_y(3)-a(4)pre_y(4);
    bx = 0;
    ay = 0;
    for k=1:NNN
        bx = bx + b(k) * pre_x(k);
    end
    for k=2:NNN
        ay = ay + a(k) * pre_y(k);  
    end  
    res_iir(i) = bx - ay;
          
    pre_y(1) = res_iir(i);
    
    for j = NNN : -1 : 2  %更新历史数据
        pre_x(j) = pre_x(j-1);
        pre_y(j) = pre_y(j-1);
    end
end


figure(3);
plot(t,res_iir);
fs=NN;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率
Y2=fft(res_iir(1:N))/N*2;%真实的幅值
figure(4)
plot(f(1:floor(N/2)),abs(Y2(1:floor(N/2))));

上面代码中，需要使用fdatool生成对应阶数的iir系数，并将系数导出到matlab的工作空间，即导出SOS矩阵和G数组。上述代码是20阶IIR滤波器的实现代码。

上面的代码是直接结构的IIR实现，首先放出一个出问题的情况：
其噪声信号为：y= 10+2cos(2pi5t-pi30/180)+2cos(2pi30t-pi30/180)+10cos(2pi100t-pi30/180)+100cos(2pi1000t-pi30/180);
使用fdatool设计20阶iir，采样频率8000Hz，低通截止频率10Hz（与8000Hz差别很大）。
代码执行效果如下：

图1 实验1的原始含噪声信号的时域及频域波形

图2 实验1的20阶iir滤波器滤波之后信号的时域及频域波形

可以从上图看到，滤波结果明显出错。
那么问题在哪呢？问题就是数字滤波器的运算中有限字长效应。可以看到，将fdatool生成的10个基本二阶结的参数转换为直接形式，得到下面的b和a系数



可以看到，其中b系数的值非常非常小，基本都是10的-几十次方，而a系数的值相对正常。那么，会导致什么问题呢？即便是使用双精度浮点进行运算，也会在运算的过程中出现运算字长不足从而使得运算出错的现象，或者iir滤波器的极点出现在单位圆上，出现极限环振荡现象，或者iir滤波器的极点直接出了单位圆，出现不稳定现象，上面的现象便是iir滤波器不稳定，不能正常工作。

那么有什么解决办法呢？办法有三种：
一、降低iir滤波器的阶数
因为在级联的过程中，每一个基本二阶结是级联相乘的，所以当iir的阶数较高时计算出来的直接结构的b、a参数很病态（非常小或者非常大），某些情况下即使使用双精度浮点计算，也是远远不够的，会出现稳定性问题。
但是当阶数较少的时候就没问题了，例如，还是上面的输入信号，将iir的阶数改为2阶iir，即一个基本二阶结，得到如下的效果（注意将matlab代码中的NNN改成2+1）：

图3 实验2的原始含噪声信号的时域及频域波形

图4 实验2的2阶iir滤波器滤波之后信号的时域及频域波形

可以看到，上面的iir输出波形，已经明显将高频的30Hz、100Hz、1000Hz的高频噪声滤除，只剩下直流分量10，和5Hz的交流低频分量，滤波效果还可以接受，试想，这只是一个2阶的iir呀，效果真好了，使用得当足以秒杀fir滤波器的效果。

但是上面的2阶iir效果还不是很好，那么再试一下4阶的iir吧！见下图：

图5 实验3的原始含噪声信号的时域及频域波形

图6 实验3的4阶iir滤波器滤波之后信号的时域及频域波形

看到没，iir滤波器的效果就是这么棒，仅用了一个4阶的iir，上面已经完全把高频分量滤除了，5Hz的原始信号是多么光滑。更加印证了iir滤波器的厉害之处。
但是，小心，万事都有其两面性。不要沉醉在iir滤波器使用较小的阶数便可以获得较好的结果的优势上，想想，我们为什么要对2阶和4阶iir进行实验？
因为高阶的iir直接实现形式不稳定！！！
此时再看，同样的滤波器参数，使用4阶的时候，计算出来的b，a参数：


可以看到，其实b参数已经特别小了，达到了10的-10次方数量级，这样的精度在计算过程中其实也比较难了，也可能会出现稳定性问题，但是相对前面20阶的iir，已经好很多了，因此上面的滤波结果也可以正常运行。

可以看到，降低iir滤波器的阶数，可以帮助稳定iir滤波器。
那么还有没有其它办法呢？
有！

二、使用直接形式的iir滤波器时，fc与fs不要相差多个数量级

上面的实验中，fs是8000Hz，而低通滤波器的截止频率仅为10Hz，相差800倍，太大了，因此使用直接形式的iir计算出来的b和a系数可能会很病态。
下面进行实验：
fs为8000Hz，而iir低通滤波器的截止频率为1000Hz，iir的阶数为20阶，使用直接形式进行计算，看能不能稳定。
此时输入信号改为：y= 10+2cos(2pi10t-pi30/180)+2cos(2pi1500t-pi30/180)+10cos(2pi2000t-pi30/180)+100cos(2pi3000t-pi30/180);

实验结果如下：

图7 实验4的原始含噪声信号的时域及频域波形

图8 实验4的20阶iir滤波器滤波之后信号的时域及频域波形

如上图，可以看到，上面的滤波效果也非常好，将1000Hz以上的高频分量已经滤除。而且即便使用20阶的iir的直接形式，其输出也能基本稳定，原因就是不让fs和fc差别特别大即可。
此时，计算得到的直接形式的iir滤波器的b，a参数如下：


可以看到，只要让fs和fc别差别很大，那么计算得到的b，a参数就不会很病态，此时iir滤波器就能正常运行，当然，前提是使用双精度浮点进行运算是这样的，如果使用单精度浮点进行运算，也不一定能稳定。。。。。。这里顺带再夸一句fir滤波器！！！yyds

上面，也是可以稳定的一种方法，但是，第一中方法没法做高阶的iir，第二种方法又不得不改变iir的截止频率fc和fs之间的关系。。。。显然，在某些情况下，这两种方法都不能被同时接受，那么，还有没有更好的办法呢？

劳动人民（此处指代研究人员=打工人）的智慧是无穷的！！！

三、使用级联结构实现高阶的iir滤波器
一般情况下，使用的都是多个基本二阶结形式进行级联从而实现高阶的iir滤波器，这样的话，便可以让每个基本二阶结的系数不那么病态了，可以避免因为计算字长的限制而导致的iir滤波器不稳定情况。
下面便是基本二阶结结构实现的高阶iir形式：


下面是基本二阶结级联形式的代码（需要SOS和G矩阵，与上面相同）：

%本程序为多个基本二阶结形式的iir滤波器的实现
%当阶数较高或fc与fs相差多个数量级时，建议使用基本二阶结的级联形式，避免计算出来的b和a系数非常小或者非常大，因运算中有限字长导致运算出错

NN=8000;
NNN=20+1;%在这里写iir滤波器的阶数+1
Sections = floor(NNN/2);%得到iir滤波器的级联形式基本二阶结的级数

t=0:1/NN:1;%采样1秒的信号
%y= 10+10*cos(2*pi*5*t-pi*30/180)+5*cos(2*pi*50*t-pi*30/180);
%y= 10+10*cos(2*pi*50*t-pi*30/180)+10*cos(2*pi*200*t-pi*30/180);
%y= 10+10*cos(2*pi*5*t-pi*30/180)+10*cos(2*pi*50*t-pi*30/180)+2*cos(2*pi*200*t-pi*30/180);
y= 1+0.5*cos(2*pi*5*t-pi*30/180)+2*cos(2*pi*15*t-pi*30/180)+1*cos(2*pi*50*t-pi*30/180)+5*cos(2*pi*2000*t-pi*30/180);
N=length(t); %样点个数
figure(1);
plot(t,y);
fs=NN;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率
Y=fft(y(1:N))/N*2;%真实的幅值
figure(2)
plot(f(1:floor(N/2)),abs(Y(1:floor(N/2))));


%[b,a]=sos2tf(SOS,G);%不需要再计算直接形式的b和a的系数了


%首先需要先定义存放数据的空间
data = zeros(Sections,6);%有3个维度。第一个维度存放基本二阶结的信息，第二个维度是当前和历史x,y信息
%data(k,1)表示第k级的x
%data(k,2)表示第k级的x（n-1）
%data(k,3)表示第k级的x（n-2）
%data(k,4)表示第k级的y
%data(k,5)表示第k级的y（n-1）
%data(k,6)表示第k级的y（n-2）

res_iir = zeros(1,N);%存放滤波的结果

for i=1:N
    data(1,1) = y(i);%将当前点的信息放进来
    for k=1:Sections
        data(k,4) =  G(k) * (SOS(k,1)*data(k,1)+SOS(k,2)*data(k,2)+SOS(k,3)*data(k,3)) -SOS(k,5)*data(k,5)-SOS(k,6)*data(k,6);
        
        data(k,6) = data(k,5);
        data(k,5) = data(k,4);
        data(k,3) = data(k,2);
        data(k,2) = data(k,1);

        if(k~=Sections)%如果不是计算到了最后一级，就将当前的输出送给下一级输入
            data(k+1,1) = data(k,4);
        else
             res_iir(i) = data(k,4);%如果是计算到了最后一级，就将当前得到结果输出
        end
    end
end



figure(3);
plot(t,res_iir);
fs=NN;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率
Y2=fft(res_iir(1:N))/N*2;%真实的幅值
figure(4)
plot(f(1:floor(N/2)),abs(Y2(1:floor(N/2))));

设计一个20阶的iir滤波器，同样，采样频率为8000Hz，低通截止频率为10Hz。信号为：y= 1+0.5cos(2pi5t-pi30/180)+2cos(2pi15t-pi30/180)+1cos(2pi50t-pi30/180)+5cos(2pi2000t-pi30/180);效果如下：

图9 实验5的原始含噪声信号的时域及频域波形

图10 实验5的20阶iir级联滤波器滤波之后信号的时域及频域波形

可以见到上图，使用级联的模式的情况下，便可以使用高阶的iir滤波器对采样的数据进行滤波，即使fs和fc相差非常大也不会带来稳定性问题。

另外，再进行一个小实验，试一下20阶和50阶的iir带阻滤波器（陷波器）效果怎么样，其中信号中含有幅度非常巨大的49、50、51Hz干扰，输入信号如下：y= 1+2cos(2pi5t-pi30/180)+100cos(2pi49t-pi30/180)+100cos(2pi50t)+100cos(2pi51t-pi30/180)+0.5cos(2pi70t-pi30/180)+5cos(2pi2000t-pi*30/180);

实验中，因为iir阶数变大，所以将数据时间长度调大到了5秒的数据长度。
陷波器阻带为45Hz~55Hz。采样频率8000Hz。

下面是实验结果：

图11 实验6的原始含噪声信号的时域及频域波形

图12 实验6的20阶iir级联带阻滤波器滤波之后信号的时域及频域波形

可以看到仅使用20阶iir滤波效果就非常好了。

下面开始50阶iir滤波器的验证：

图13 实验7的原始含噪声信号的时域及频域波形

图14 实验7的50阶iir级联带阻滤波器滤波之后信号的时域及频域波形

可以看到上面20阶和50阶iir滤波器的对比，发现50阶iir滤波器需要更长的时间才能稳定，即其反馈环路更长，然而他们的实际滤波效果并差不多，因此，给我们的启示是：在实际的系统中，应当合理选择iir滤波器的阶数，并不是越高越好。

下面展示一下两个滤波器（20阶和50阶）的滤波后稳定波形的细节放大图：

图15 20阶iir级联带阻滤波器滤波之后信号的时域波形细节图

图16 50阶iir级联带阻滤波器滤波之后信号的时域波形细节图

对比图15和图16，发现无论是20阶还是50阶iir滤波器稳定之后的滤波效果都非常好，输出的波形里面只包含了幅度为1的直流分量、幅度为2的5Hz低频分量和幅度为0.5的70Hz高频分量。两个图的波形基本一致。然而50阶的iir滤波器需要更长时间稳定。

至此，无论是fir滤波器还是iir滤波器，我们均从其具体实现方法方面进行了一些实验，发现、解决了一些问题，对于里面的代码，可以非常方便地修改为C语言代码应用于嵌入式开发中。

这两篇小文档作为自己进行这些实验的一个小结和心得吧，可以留给以后的自己去看。

但是，这两种滤波器仅当噪声的频率与实际信号频率之间不重叠的时候有用。对于一些随机时间信号噪声，这种经典滤波器是无能为力的，然后就需要一些现代统计最优滤波器去对这些噪声进行滤波，例如维纳滤波器、卡尔曼滤波器、自适应滤波等。这些统计最优滤波器的相关实验等待后面更新吧~
最后，再次感谢本学期“信号处理与数据分析”课程的邱老师，这门课程真的很有用，特别是应用于实践之后。