实验原理IIR数字滤波器一、脉冲响应不变法变换原理将模拟滤波器的s平面变换成数字滤波器的平面,而将模拟滤波器映射成数字滤波器。
MATLAB 信号处理工具箱中提供了IIR数字滤波器设计的函数,常用的函数:IIR滤波器的阶数选择Buttord----巴特沃斯滤波器阶数选择cheb1ord-----切比雪夫I型滤波器阶数选择cheb2ord ----切比雪夫Ⅱ型滤波器阶数选择IIR滤波器的设计Buttrer-----巴特沃斯滤波器设计cheby1-----切比雪夫I型滤波器设计cheby2 ----切比雪夫Ⅱ型滤波器设计maxflat----通用的巴特沃斯低通滤波器设计二、巴特沃斯滤波器设计巴特沃斯滤波器是通带,阻带都单调衰减的滤波器。
(1)调用buttord函数确定巴特沃斯滤波器的阶数,格式[N,Wc]=buttord(Wp,Ws,Ap,As)(2)调用butter函数设计巴特沃斯滤波器,格式[b,a]=butter(N,Wc,options)利用以上两个函数可以设计出模拟滤波器,格式为[N,Wc]=buttord(Wp,Ws,Ap,As,’s’)[b,a]=butter(N,Wc,options.’s’)三、切比雪夫I型滤波器的设计切比雪夫I型滤波器为通带波纹控制器;在通带呈现波纹特性,阻带单调衰减。
[N,Wc]= cheb1ord (Wp,Ws,Ap,As)[b,a]= cheby1 (N,Ap,Wc,options)四、切比雪夫Ⅱ型滤波器的设计切比雪夫Ⅱ型滤波器为阻带波纹控制器;在阻带呈现波纹特性,通带单调衰减。
[N,Wc]= cheb2ord (Wp,Ws,Ap,As)[b,a]= cheby2 (N,As,Wc,options)已知模拟滤波器。
可以利用脉冲响应不变法转换函数impinvar将其变为数字滤波器,调用格式为[bz,az]=impinvar(b,a,Fs).五、双线性变换原理采用非线性频率压缩方法,克服了脉冲响应不变法产生频率响应的混叠失真,使是s平面与z平面建立了一一对应的单值关系,消除多值变换性,频谱混叠现象。
已知模拟滤波器,可以利用双线性变换函数bilinear将其变换为数字滤波器,调用格式为[bz,az]=bilinear(b,a,Fs).FIR窗函数法设计FIR数字滤波器FIR滤波器设计需使频率响应H(ejw)逼近所要求的理想频率响应,窗函数法设计FIR数字滤波器是在时域中进行的,用窗函数截取无限长的hd(n),这样H(ejw)逼近与理想值。
一旦选取了窗函数,其指标就是给定的,所以由窗函数设计FIR滤波器就是有阻带衰减指标确定用什么窗,由过渡带宽估计窗函数的长达N。
常用的有:hd=boxcar(N) ht=triang(N) hd=hanning(N)hd=hamming(N) hd=blackman(N) hd=kaiser(N,β)MATLAB中提供的fir可以用来设计FIR滤波器,调用格式为h=fir1(M,Wc,’ftype’,window)实验内容一1、要求通带截止频率fp=3KHz,通带最大衰减ap=1dB,阻带截止fs=4.5kHz,阻带最小衰减as=15dB,采样频率fc=30kHz,用脉冲响应不变法设计一个切比雪夫低通滤波器,并图示滤波器的振幅特性,检验wp,ws对应的衰减。
(1)>> wp=6*pi*10^3;ws=9*pi*10^3;ap=1;as=15;>> Fs=30*10^3;>> wp1=wp/Fs;ws1=ws/Fs; %参数设置>> [N,WC]=cheb1ord(wp,ws,ap,as,'s');>> [b,a]=cheby1(N,ap,WC,'s'); %采用切比雪夫I型滤波器设计>> [bz,az]=impinvar(b,a,Fs); %采用脉冲响应不变法>> w0=[wp1,ws1];>> Hx=freqz(bz,az,w0);>> [H,W]=freqz(bz,az);>> dbHx=-20*log10(abs(Hx)/max(abs(H))); %显示幅频特性>> plot(W,abs(H));>> xlabel('相对频率');ylabel('幅频');>> grid>>得到的结果如下>> bzbz =-0.0000 0.0054 0.0181 0.0040 0>> azaz =1.0000 -3.0591 3.8323 -2.2919 0.5495>> dbHxdbHx =1.0005 21.5790结果显示ap 略大于1,as 符合要求,得到的图形文件如下所示0.511.522.533.500.20.40.60.811.21.4相对频率幅频(2)>> wp=6*pi*10^3;ws=9*pi*10^3;ap=1;as=15; >> Fs=30*10^3;>> wp1=wp/Fs;ws1=ws/Fs; %参数设置 >> [N,WC]=cheb2ord(wp,ws,ap,as,'s');>> [b,a]=cheby2(N,as,WC,'s'); %采用切比雪夫Ⅱ型滤波器设计 >> [bz,az]=impinvar(b,a,Fs); %采用脉冲响应不变法 >> w0=[wp1,ws1];>> Hx=freqz(bz,az,w0); >> [H,W]=freqz(bz,az);>> dbHx=-20*log10(abs(Hx)/max(abs(H))); %显示幅频特性 >> plot(W,abs(H));>> xlabel('相对频率');ylabel('幅频'); >> grid >> bz得到的结果及图形文件如下 bz =-0.3682 0.9359 -0.8472 0.4014 0.0000>> azaz =1.0000 -1.9388 1.7388 -0.6968 0.1261>> dbHxdbHx =0.0449 8.796700.511.522.533.50.40.50.60.70.80.911.11.21.3相对频率幅频显然,这种设计不符合要求.2、 用双线性变换法设计一个切比雪夫I 型数字高通滤波器,技术指标为:采样频率fc=2kHz,通带截止频率fp=700Hz,通带最大衰减ap<=1dB,阻带边缘频率fs=500Hz,阻带最小衰减as>=32Db. 程序如下>> wp=2000*pi;ws=500*pi;ap=1;as=32; >> wp=1400*pi;ws=1000*pi;ap=1;as=32; >> Fs=2000;>> wp1=wp/Fs;ws1=ws/Fs;>> omp1=2*Fs*tan(wp1/2);omps=2*Fs*tan(ws1/2);>> [N,WC]=cheb1ord(omp1,omps,ap,as,'s'); %采用切比雪夫Ⅰ型滤波器设计 >> [b,a]=cheby1(N,ap,WC,'high','s');>> [bz,az]=bilinear(b,a,Fs); %采用双线性变换法 Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 8.974782e-024. > In bilinear at 8900.511.522.533.5相对频率幅频Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 8.974782e-024. > In bilinear at 90 >> w0=[wp1,ws1]; >> Hx=freqz(bz,az,w0); >> [H,W]=freqz(bz,az);>> dbHx=-20*log10(abs(Hx)/max(abs(H))); >> plot(W,abs(H));>> xlabel('相对频率');ylabel('幅频'); >> grid 得到的结果及图形文件如下>> bzbz =0.0084 -0.0335 0.0502 -0.0335 0.0084 >>>> azaz =1.00002.3741 2.7057 1.5917 0.4103>> dbHxdbHx =1.0000 33.1098 可见,ap,as 符合设计要求二1、窗函数法设计低通滤波器,w=0.4,(1)N=26,分别利用矩形窗、汉宁窗和blackman 窗设计该滤波器,且滤波器具有线性相位。
绘出脉冲响应h(n)及滤波器的频率响应。
>> N=26; >> wc=0.4; >> nn=[0:25];>> h1=fir1(25,wc,boxcar(N)); %利用矩形窗 >> [H,W]=freqz(h1,1); >> subplot(311),plot(nn,h1) >> title('矩形窗频率响应') >> xlabel('nn'),ylabel('h1')510152025矩形窗频率响应nnh 10.050.10.150.20.250.30.350.40.450.5矩形窗幅频特性响应wa b s矩形窗相频特性响应wa n g l e0510152025blackman 窗频率响应nnh 300.050.10.150.20.250.30.350.40.450.5blackman 窗幅频特性响应wa b sblackman 窗相频特性响应wa n g l e510152025汉宁窗频率响应nnh 20.050.10.150.20.250.30.350.40.450.5汉宁窗幅频特性响应wa b s0.050.10.150.20.250.30.350.40.450.5汉宁窗相频特性响应wa n g l e>> subplot(312),plot(W/2/pi,abs(H)) >> title('矩形窗幅频特性响应') >> xlabel('w'),ylabel('abs')>> subplot(313),plot(W/2/pi,angle(H)*180/pi) >> title('矩形窗相频特性响应') >>xlabel('w'),ylabel('angle')>> N=26; %利用汉宁窗 >> wc=0.4; >> nn=[0:25];>> h2=fir1(25,wc,hanning(N)); >> [H,W]=freqz(h2,1); >> subplot(311),plot(nn,h2) >> title('汉宁窗频率响应') >> xlabel('nn'),ylabel('h2')>> subplot(312),plot(W/2/pi,abs(H)) >> title('汉宁窗幅频特性响应') >> xlabel('w'),ylabel('abs')>> subplot(313),plot(W/2/pi,angle(H)*180/pi) >> title('汉宁窗相频特性响应') >> xlabel('w'),ylabel('angle')>> N=26; %利用blackman 窗>> wc=0.4;>> nn=[0:25];>> h3=fir1(25,wc,blackman(N)); >> [H,W]=freqz(h3,1);>> subplot(311),plot(nn,h3)>> title('blackman 窗频率响应')>> xlabel('nn'),ylabel('h3') >> subplot(312),plot(W/2/pi,abs(H))>> title('blackman 窗幅频特性响应')>> xlabel('w'),ylabel('abs') >> subplot(313),plot(W/2/pi,angle(H)*180/pi)(2)增加N 至N=66,观察过渡带和最大尖峰值的变化。