IIR滤波器设计(调用MATLAB IIR函数来实现)


956人阅读 评论(0) 收藏 举报
转载请注明文章来源 – http://blog.csdn.net/v_hyx ,请勿用于任何商业用途

      对于滤波器设计,以前虽然学过相关的理论(现代数字信号处理和DSP设计),但一直不求甚解,也没用过。趁着最近使用了一下,就来重学一回,温故而知新。
      先来说说IIR滤波器设计,理论与原理参考如下博客,写得简明易懂,不错。
http://blog.csdn.net/thnh169/article/details/9034483  [数字信号处理]IIR滤波器基础
http://blog.csdn.net/thnh169/article/details/9069475  [数字信号处理]IIR滤波器的间接设计(C代码)
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
  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');

效果图








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

参考资料:
1.  http://www.cnblogs.com/sunev/archive/2011/11/23/2260579.html 基于Matlab的FIR滤波器设计与实现
2.  http://blog.csdn.net/thnh169/article/details/9034483  [数字信号处理]IIR滤波器基础
3.  http://blog.csdn.net/thnh169/article/details/9069475  [数字信号处理]IIR滤波器的间接设计(C代码)
4.  http://blog.csdn.net/thnh169/article/details/9076283  [数字信号处理]IIR滤波器的直接设计(C代码)
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值