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

  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');