使用matlab设计iir滤波器并自行编写代码实现iir滤波器(可对应于C语言应用在嵌入式系统中) |
您所在的位置:网站首页 › iir滤波器群延迟 › 使用matlab设计iir滤波器并自行编写代码实现iir滤波器(可对应于C语言应用在嵌入式系统中) |
对于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滤波器滤波之后信号的时域及频域波形 那么有什么解决办法呢?办法有三种: 一、降低iir滤波器的阶数 因为在级联的过程中,每一个基本二阶结是级联相乘的,所以当iir的阶数较高时计算出来的直接结构的b、a参数很病态(非常小或者非常大),某些情况下即使使用双精度浮点计算,也是远远不够的,会出现稳定性问题。 但是当阶数较少的时候就没问题了,例如,还是上面的输入信号,将iir的阶数改为2阶iir,即一个基本二阶结,得到如下的效果(注意将matlab代码中的NNN改成2+1): 图3 实验2的原始含噪声信号的时域及频域波形 但是上面的2阶iir效果还不是很好,那么再试一下4阶的iir吧!见下图: 图5 实验3的原始含噪声信号的时域及频域波形 图6 实验3的4阶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的原始含噪声信号的时域及频域波形 上面,也是可以稳定的一种方法,但是,第一中方法没法做高阶的iir,第二种方法又不得不改变iir的截止频率fc和fs之间的关系。。。。显然,在某些情况下,这两种方法都不能被同时接受,那么,还有没有更好的办法呢? 劳动人民(此处指代研究人员=打工人)的智慧是无穷的!!! 三、使用级联结构实现高阶的iir滤波器 一般情况下,使用的都是多个基本二阶结形式进行级联从而实现高阶的iir滤波器,这样的话,便可以让每个基本二阶结的系数不那么病态了,可以避免因为计算字长的限制而导致的iir滤波器不稳定情况。 下面便是基本二阶结结构实现的高阶iir形式: 设计一个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的原始含噪声信号的时域及频域波形 另外,再进行一个小实验,试一下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的原始含噪声信号的时域及频域波形 下面开始50阶iir滤波器的验证: 图13 实验7的原始含噪声信号的时域及频域波形 图14 实验7的50阶iir级联带阻滤波器滤波之后信号的时域及频域波形 可以看到上面20阶和50阶iir滤波器的对比,发现50阶iir滤波器需要更长的时间才能稳定,即其反馈环路更长,然而他们的实际滤波效果并差不多,因此,给我们的启示是:在实际的系统中,应当合理选择iir滤波器的阶数,并不是越高越好。 下面展示一下两个滤波器(20阶和50阶)的滤波后稳定波形的细节放大图: 图15 20阶iir级联带阻滤波器滤波之后信号的时域波形细节图 至此,无论是fir滤波器还是iir滤波器,我们均从其具体实现方法方面进行了一些实验,发现、解决了一些问题,对于里面的代码,可以非常方便地修改为C语言代码应用于嵌入式开发中。 这两篇小文档作为自己进行这些实验的一个小结和心得吧,可以留给以后的自己去看。 但是,这两种滤波器仅当噪声的频率与实际信号频率之间不重叠的时候有用。对于一些随机时间信号噪声,这种经典滤波器是无能为力的,然后就需要一些现代统计最优滤波器去对这些噪声进行滤波,例如维纳滤波器、卡尔曼滤波器、自适应滤波等。这些统计最优滤波器的相关实验等待后面更新吧~ 最后,再次感谢本学期“信号处理与数据分析”课程的邱老师,这门课程真的很有用,特别是应用于实践之后。 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |