IIR滤波器设计（调用MATLAB IIR函数来实现）

ゝ一世哀愁。 2022-08-25 14:53 207阅读 0赞

#     #

2013-08-23 18:20  956人阅读 [评论][Link 1](0) 收藏 [举报][Link 2]

[fft][] [filter][] [MATLAB][] [IIR][]

转载请注明文章来源 – http://blog.csdn.net/v\_hyx ，请勿用于任何商业用途

对于滤波器设计，以前虽然学过相关的理论（现代数字信号处理和DSP设计），但一直不求甚解，也没用过。趁着最近使用了一下，就来重学一回，温故而知新。

先来说说IIR滤波器设计，理论与原理参考如下博客，写得简明易懂，不错。

[http://blog.csdn.net/thnh169/article/details/9034483][http_blog.csdn.net_thnh169_article_details_9034483] \[数字信号处理\]IIR滤波器基础

[http://blog.csdn.net/thnh169/article/details/9069475][http_blog.csdn.net_thnh169_article_details_9069475] \[数字信号处理\]IIR滤波器的间接设计（C代码）

[http://blog.csdn.net/thnh169/article/details/9076283][http_blog.csdn.net_thnh169_article_details_9076283] \[数字信号处理\]IIR滤波器的直接设计(C代码)

一般IIR设计可分三种：间接设计（原型转换设计）、直接设计、直接使用工具软件如MATLAB的IIR函数设计。前2种方法，上面的博客已经写得很清楚，理论比较多，设计还是很复杂的。但在实际工程应用中，多采用MATLAB的IIR函数或者FDATOOL工具进行，非常方便快捷。

OK，来个示例来说明采用MATLAB的IIR函数设计过程，花一会儿的功夫就可以快速入门，so easy！废话不多说，直接上MATLAB IIR.m文件，最后附上效果图。

IIR.m

**\[plain\]** [ view plain][view plain] [copy][view plain] [print][view plain] [?][view plain]

1.  % IIR滤波器设计
2.  % 目的：设计一个采样频率为1000Hz、通带截止频率为50Hz、阻带截止频率为100Hz的低通滤波器，并要求通带最大衰减为1dB，阻带最小衰减为60dB。
3.  
4.  clc;clear;close all;
5.  
6.  % 1. 产生信号（频率为10Hz和100Hz的正弦波叠加）
7.  Fs=1000; % 采样频率1000Hz
8.  t=0:1/Fs:1;
9.  s10=sin(20\*pi\*t); % 产生10Hz正弦波
10. s100=sin(200\*pi\*t); % 产生100Hz正弦波
11. s=s10+s100; % 信号叠加
12. 
13. figure(1); % 画图
14. subplot(2,1,1);plot(s);grid;
15. title('原始信号');
16. 
17. % 2. FFT分析信号频谱
18. len = 512;
19. y=fft(s,len); % 对信号做len点FFT变换
20. f=Fs\*(0:len/2 - 1)/len;
21. 
22. subplot(2,1,2);plot(f,abs(y(1:len/2)));grid;
23. title('原始信号频谱')
24. xlabel('Hz');ylabel('幅值');
25. 
26. % 3. IIR滤波器设计
27. N=0; % 阶数
28. Fp=50; % 通带截止频率50Hz
29. Fc=100; % 阻带截止频率100Hz
30. Rp=1; % 通带波纹最大衰减为1dB
31. Rs=60; % 阻带衰减为60dB
32. 
33. % 3.0 计算最小滤波器阶数
34. na=sqrt(10^(0.1\*Rp)-1);
35. ea=sqrt(10^(0.1\*Rs)-1);
36. N=ceil(log10(ea/na)/log10(Fc/Fp));
37. 
38. % 3.1 巴特沃斯滤波器
39. Wn=Fp\*2/Fs;
40. \[Bb Ba\]=butter(N,Wn,'low'); % 调用MATLAB butter函数快速设计滤波器
41. \[BH,BW\]=freqz(Bb,Ba); % 绘制频率响应曲线
42. Bf=filter(Bb,Ba,s); % 进行低通滤波
43. By=fft(Bf,len); % 对信号f1做len点FFT变换
44. 
45. figure(2); % 画图
46. subplot(2,1,1);plot(t,s10,'blue',t,Bf,'red');grid;
47. legend('10Hz原始信号','巴特沃斯滤波器滤波后');
48. subplot(2,1,2);plot(f,abs(By(1:len/2)));grid;
49. title('巴特沃斯低通滤波后信号频谱');
50. xlabel('Hz');ylabel('幅值');
51. 
52. % 3.2 切比雪夫I型滤波器
53. \[C1b C1a\]=cheby1(N,Rp,Wn,'low'); % 调用MATLAB cheby1函数快速设计低通滤波器
54. \[C1H,C1W\]=freqz(C1b,C1a); % 绘制频率响应曲线
55. C1f=filter(C1b,C1a,s); % 进行低通滤波
56. C1y=fft(C1f,len); % 对滤波后信号做len点FFT变换
57. 
58. figure(3); % 画图
59. subplot(2,1,1);plot(t,s10,'blue',t,C1f,'red');grid;
60. legend('10Hz原始信号','切比雪夫I型滤波器滤波后');
61. subplot(2,1,2);plot(f,abs(C1y(1:len/2)));grid;
62. title('切比雪夫I型滤波后信号频谱');
63. xlabel('Hz');ylabel('幅值');
64. 
65. % 3.3 切比雪夫II型滤波器
66. \[C2b C2a\]=cheby2(N,Rs,Wn,'low'); % 调用MATLAB cheby2函数快速设计低通滤波器
67. \[C2H,C2W\]=freqz(C2b,C2a); % 绘制频率响应曲线
68. C2f=filter(C2b,C2a,s); % 进行低通滤波
69. C2y=fft(C2f,len); % 对滤波后信号做len点FFT变换
70. 
71. figure(4); % 画图
72. subplot(2,1,1);plot(t,s10,'blue',t,C2f,'red');grid;
73. legend('10Hz原始信号','切比雪夫II型滤波器滤波后');
74. subplot(2,1,2);plot(f,abs(C2y(1:len/2)));grid;
75. title('切比雪夫II型滤波后信号频谱');
76. xlabel('Hz');ylabel('幅值');
77. 
78. % 3.4 椭圆滤波器
79. \[Eb Ea\]=ellip(N,Rp,Rs,Wn,'low'); % 调用MATLAB ellip函数快速设计低通滤波器
80. \[EH,EW\]=freqz(Eb,Ea); % 绘制频率响应曲线
81. Ef=filter(Eb,Ea,s); % 进行低通滤波
82. Ey=fft(Ef,len); % 对滤波后信号做len点FFT变换
83. 
84. figure(5); % 画图
85. subplot(2,1,1);plot(t,s10,'blue',t,Ef,'red');grid;
86. legend('10Hz原始信号','椭圆滤波器滤波后');
87. subplot(2,1,2);plot(f,abs(Ey(1:len/2)));grid;
88. title('椭圆滤波后信号频谱');
89. xlabel('Hz');ylabel('幅值');
90. 
91. % 3.5 yulewalk滤波器
92. fyule=\[0 Wn Fc\*2/Fs 1\]; % 在此进行的是简单设计，实际需要多次仿真取最佳值
93. myule=\[1 1 0 0\]; % 在此进行的是简单设计，实际需要多次仿真取最佳值
94. \[Yb Ya\]=yulewalk(N,fyule,myule); % 调用MATLAB yulewalk函数快速设计低通滤波器
95. \[YH,YW\]=freqz(Yb,Ya); % 绘制频率响应曲线
96. Yf=filter(Yb,Ya,s); % 进行低通滤波
97. Yy=fft(Yf,len); % 对滤波后信号做len点FFT变换
98. 
99. figure(6); % 画图
100. subplot(2,1,1);plot(t,s10,'blue',t,Yf,'red');grid;
101. legend('10Hz原始信号','yulewalk滤波器滤波后');
102. subplot(2,1,2);plot(f,abs(Yy(1:len/2)));grid;
103. title('yulewalk滤波后信号频谱');
104. xlabel('Hz');ylabel('幅值');
105. 
106. % 4. 各个滤波器的幅频响应对比分析
107. A1 = BW\*Fs/(2\*pi);
108. A2 = C1W\*Fs/(2\*pi);
109. A3 = C2W\*Fs/(2\*pi);
110. A4 = EW\*Fs/(2\*pi);
111. A5 = YW\*Fs/(2\*pi);
112. 
113. figure(7); % 画图
114. subplot(1,1,1);plot(A1,abs(BH),A2,abs(C1H),A3,abs(C2H),A4,abs(EH),A5,abs(YH));grid;
115. xlabel('频率／Hz');
116. ylabel('频率响应幅度');
117. legend('butter','cheby1','cheby2','ellip','yulewalk');

% IIR滤波器设计
    % 目的：设计一个采样频率为1000Hz、通带截止频率为50Hz、阻带截止频率为100Hz的低通滤波器，并要求通带最大衰减为1dB，阻带最小衰减为60dB。
    
    clc;clear;close all;
    
    % 1. 产生信号（频率为10Hz和100Hz的正弦波叠加）
    Fs=1000; % 采样频率1000Hz
    t=0:1/Fs:1;
    s10=sin(20*pi*t); % 产生10Hz正弦波
    s100=sin(200*pi*t); % 产生100Hz正弦波
    s=s10+s100; % 信号叠加
    
    figure(1); % 画图
    subplot(2,1,1);plot(s);grid;
    title('原始信号');
    
    % 2. FFT分析信号频谱
    len = 512;
    y=fft(s,len);  % 对信号做len点FFT变换
    f=Fs*(0:len/2 - 1)/len;
    
    subplot(2,1,2);plot(f,abs(y(1:len/2)));grid;
    title('原始信号频谱')
    xlabel('Hz');ylabel('幅值');
    
    % 3. IIR滤波器设计
    N=0; % 阶数
    Fp=50; % 通带截止频率50Hz
    Fc=100; % 阻带截止频率100Hz
    Rp=1; % 通带波纹最大衰减为1dB
    Rs=60; % 阻带衰减为60dB
    
    % 3.0 计算最小滤波器阶数
    na=sqrt(10^(0.1*Rp)-1);
    ea=sqrt(10^(0.1*Rs)-1);
    N=ceil(log10(ea/na)/log10(Fc/Fp));
    
    % 3.1 巴特沃斯滤波器
    Wn=Fp*2/Fs;
    [Bb Ba]=butter(N,Wn,'low'); % 调用MATLAB butter函数快速设计滤波器
    [BH,BW]=freqz(Bb,Ba); % 绘制频率响应曲线
    Bf=filter(Bb,Ba,s); % 进行低通滤波
    By=fft(Bf,len);  % 对信号f1做len点FFT变换
    
    figure(2); % 画图
    subplot(2,1,1);plot(t,s10,'blue',t,Bf,'red');grid;
    legend('10Hz原始信号','巴特沃斯滤波器滤波后');
    subplot(2,1,2);plot(f,abs(By(1:len/2)));grid;
    title('巴特沃斯低通滤波后信号频谱');
    xlabel('Hz');ylabel('幅值');
    
    % 3.2 切比雪夫I型滤波器
    [C1b C1a]=cheby1(N,Rp,Wn,'low'); % 调用MATLAB cheby1函数快速设计低通滤波器
    [C1H,C1W]=freqz(C1b,C1a); % 绘制频率响应曲线
    C1f=filter(C1b,C1a,s); % 进行低通滤波
    C1y=fft(C1f,len);  % 对滤波后信号做len点FFT变换
    
    figure(3); % 画图
    subplot(2,1,1);plot(t,s10,'blue',t,C1f,'red');grid;
    legend('10Hz原始信号','切比雪夫I型滤波器滤波后');
    subplot(2,1,2);plot(f,abs(C1y(1:len/2)));grid;
    title('切比雪夫I型滤波后信号频谱');
    xlabel('Hz');ylabel('幅值');
    
    % 3.3 切比雪夫II型滤波器
    [C2b C2a]=cheby2(N,Rs,Wn,'low'); % 调用MATLAB cheby2函数快速设计低通滤波器
    [C2H,C2W]=freqz(C2b,C2a); % 绘制频率响应曲线
    C2f=filter(C2b,C2a,s); % 进行低通滤波
    C2y=fft(C2f,len);  % 对滤波后信号做len点FFT变换
    
    figure(4); % 画图
    subplot(2,1,1);plot(t,s10,'blue',t,C2f,'red');grid;
    legend('10Hz原始信号','切比雪夫II型滤波器滤波后');
    subplot(2,1,2);plot(f,abs(C2y(1:len/2)));grid;
    title('切比雪夫II型滤波后信号频谱');
    xlabel('Hz');ylabel('幅值');
    
    % 3.4 椭圆滤波器
    [Eb Ea]=ellip(N,Rp,Rs,Wn,'low'); % 调用MATLAB ellip函数快速设计低通滤波器
    [EH,EW]=freqz(Eb,Ea); % 绘制频率响应曲线
    Ef=filter(Eb,Ea,s); % 进行低通滤波
    Ey=fft(Ef,len);  % 对滤波后信号做len点FFT变换
    
    figure(5); % 画图
    subplot(2,1,1);plot(t,s10,'blue',t,Ef,'red');grid;
    legend('10Hz原始信号','椭圆滤波器滤波后');
    subplot(2,1,2);plot(f,abs(Ey(1:len/2)));grid;
    title('椭圆滤波后信号频谱');
    xlabel('Hz');ylabel('幅值');
    
    % 3.5 yulewalk滤波器
    fyule=[0 Wn Fc*2/Fs 1]; % 在此进行的是简单设计，实际需要多次仿真取最佳值
    myule=[1 1 0 0]; % 在此进行的是简单设计，实际需要多次仿真取最佳值
    [Yb Ya]=yulewalk(N,fyule,myule); % 调用MATLAB yulewalk函数快速设计低通滤波器
    [YH,YW]=freqz(Yb,Ya); % 绘制频率响应曲线
    Yf=filter(Yb,Ya,s); % 进行低通滤波
    Yy=fft(Yf,len);  % 对滤波后信号做len点FFT变换
    
    figure(6); % 画图
    subplot(2,1,1);plot(t,s10,'blue',t,Yf,'red');grid;
    legend('10Hz原始信号','yulewalk滤波器滤波后');
    subplot(2,1,2);plot(f,abs(Yy(1:len/2)));grid;
    title('yulewalk滤波后信号频谱');
    xlabel('Hz');ylabel('幅值');
    
    % 4. 各个滤波器的幅频响应对比分析
    A1 = BW*Fs/(2*pi);
    A2 = C1W*Fs/(2*pi);
    A3 = C2W*Fs/(2*pi);
    A4 = EW*Fs/(2*pi);
    A5 = YW*Fs/(2*pi);
    
    figure(7); % 画图
    subplot(1,1,1);plot(A1,abs(BH),A2,abs(C1H),A3,abs(C2H),A4,abs(EH),A5,abs(YH));grid;
    xlabel('频率／Hz');
    ylabel('频率响应幅度');
    legend('butter','cheby1','cheby2','ellip','yulewalk');

效果图

![SouthEast][]

![SouthEast 1][]

![SouthEast 2][]

![SouthEast 3][]

![SouthEast 4][]

![SouthEast 5][]

![SouthEast 6][]

转载请注明文章来源 – http://blog.csdn.net/v\_hyx ，请勿用于任何商业用途

参考资料：

1. [http://www.cnblogs.com/sunev/archive/2011/11/23/2260579.html][http_www.cnblogs.com_sunev_archive_2011_11_23_2260579.html] 基于Matlab的FIR滤波器设计与实现

2. [http://blog.csdn.net/thnh169/article/details/9034483][http_blog.csdn.net_thnh169_article_details_9034483] \[数字信号处理\]IIR滤波器基础

3. [http://blog.csdn.net/thnh169/article/details/9069475][http_blog.csdn.net_thnh169_article_details_9069475] \[数字信号处理\]IIR滤波器的间接设计（C代码）

4. [http://blog.csdn.net/thnh169/article/details/9076283][http_blog.csdn.net_thnh169_article_details_9076283] \[数字信号处理\]IIR滤波器的直接设计(C代码)

[Link 1]: http://blog.csdn.net/share_happy_1984/article/details/10236723#comments
[Link 2]: http://blog.csdn.net/share_happy_1984/article/details/10236723#report
[fft]: http://www.csdn.net/tag/fft
[filter]: http://www.csdn.net/tag/filter
[MATLAB]: http://www.csdn.net/tag/MATLAB
[IIR]: http://www.csdn.net/tag/IIR
[http_blog.csdn.net_thnh169_article_details_9034483]: http://blog.csdn.net/thnh169/article/details/9034483
[http_blog.csdn.net_thnh169_article_details_9069475]: http://blog.csdn.net/thnh169/article/details/9069475
[http_blog.csdn.net_thnh169_article_details_9076283]: http://blog.csdn.net/thnh169/article/details/9076283
[view plain]: http://blog.csdn.net/share_happy_1984/article/details/10236723#
[SouthEast]: /images/20220824/e6aecb9fe2524ccd889873f1befda02d.png
[SouthEast 1]: /images/20220824/52eeb30eddc5492f8c9ff2caaf9944cf.png
[SouthEast 2]: /images/20220824/cd1ea6a88291473da777d5a66e9b74d3.png
[SouthEast 3]: /images/20220824/b9e2fd74c65645fe99cd249184603adb.png
[SouthEast 4]: /images/20220824/7742ba37ec894535842d87b7f6cb3b32.png
[SouthEast 5]: /images/20220824/501ae3da979d410889f4667cb66a3dd2.png
[SouthEast 6]: /images/20220824/2da9769f6ea346ec94570096d16cebba.png
[http_www.cnblogs.com_sunev_archive_2011_11_23_2260579.html]: http://www.cnblogs.com/sunev/archive/2011/11/23/2260579.html