DSP实验4巴特沃斯滤波器的设计与实现

实验四巴特沃斯数字滤波器的设计与实现1．数字滤波器的设计参数滤波器的4个重要的通带、阻带参数为：pf：通带截止频率（Hz）sf：阻带起始频率（Hz）pR：通带内波动（dB），即通带内所允许的最大衰减；sR：阻带内最小衰减设采样速率（即奈奎斯特速率）为Nf，将上述参数中的频率参数转化为归一化角频率参数：p：通带截止角频率（rad/s），)2//(Nppff；s：阻带起始角频率（rad/s），)2//(Nssff通过以上参数就可以进行离散滤波器的设计。低通滤波器情况：采样频率为8000Hz，要求通带截止频率为1500Hz，阻带起始频率为2000Hz，通带内波动3dB，阻带内最小衰减为50dB，则p=1500/4000，s=2000/4000，pR=3dB，sR=50dB。高通滤波器情况：采样频率为8000Hz，要求通带截止频率为1500Hz，阻带起始频率为1000Hz，通带内波动3dB，阻带内最小衰减为65dB，则p=1500/4000，s=1000/4000，pR=3dB，sR=65dB。带通滤波器情况：采样频率为8000Hz，要求通带截止频率为[800Hz，1500Hz]，阻带起始频率为[500Hz，1800Hz]，通带内波动3dB，阻带内最小衰减为45dB，则p=[800/4000，1500/4000]，s=[500/4000，1800/4000]，pR=3dB，sR=45dB。带阻滤波器情况：采样频率为8000Hz，要求通带截止频率为[800Hz，1500Hz]，阻带起始频率为[1000Hz，1300Hz]，通带内波动3dB，阻带内最小衰减为55dB，则p=[800/4000，1500/4000]，s=[1000/4000，1300/4000]，pR=3dB，sR=45dB。2．巴特沃斯滤波器设计1）巴特沃斯滤波器阶数的选择：在已知设计参数p，s，pR，sR之后，可利用“buttord”命令可求出所需要的滤波器的阶数和3dB截止频率，其格式为：[n，Wn]=buttord[Wp，Ws，Rp，Rs]，其中Wp，Ws，Rp，Rs分别为通带截止频率、阻带起始频率、通带内波动、阻带内最小衰减。返回值n为滤波器的最低阶数，Wn为3dB截止频率。2）巴特沃斯滤波器系数计算：由巴特沃斯滤波器的阶数n以及3dB截止频率Wn可以计算出对应传递函数H(z)的分子分母系数，MATLAB提供的命令如下：巴特沃斯低通滤波器系数计算：[b，a]=butter(n,Wn)，其中b为H(z)的分子多项式系数，a为H(z)的分母多项式系数巴特沃斯高通滤波器系数计算：[b，a]=butter(n，Wn，’High’)巴特沃斯带通滤波器系数计算：[b，a]=butter(n，[W1，W2])，其中[W1，W2]为截止频率，是2元向量，需要注意的是该函数返回的是2*n阶滤波器系数。巴特沃斯带阻滤波器系数计算：[b，a]=butter(ceil(n/2)，[W1，W2]，’stop’)，其中[W1，W2]为截止频率，是2元向量，需要注意的是该函数返回的也是2*n阶滤波器系数。3．巴特沃斯滤波器设计实例1）采样速率为8000Hz，要求设计一个低通滤波器，pf=2100Hz，sf=2500Hz，pR=3dB，sR=25dB。程序如下：fn=8000;fp=2100;fs=2500;Rp=3;Rs=25;Wp=fp/(fn/2);%计算归一化角频率Ws=fs/(fn/2);[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数[H,F]=freqz(b,a,1000,8000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)subplot(2,1,1)plot(F,20*log10(abs(H)))xlabel('Frequency(Hz)');ylabel('Magnitude(dB)')title('低通滤波器')axis([04000-303]);gridonsubplot(2,1,2)pha=angle(H)*180/pi;plot(F,pha);gridon运行结果如下：2）采样速率为8000Hz，要求设计一个高通滤波器，pf=1000Hz，sf=700Hz，pR=3dB，sR=20dB。程序如下：fn=8000;fp=1000;fs=700;Rp=3;Rs=20;Wp=fp/(fn/2);%计算归一化角频率Ws=fs/(fn/2);[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率[b,a]=butter(n,Wn,’high’);%计算H(z)分子、分母多项式系数[H,F]=freqz(b,a,1000,8000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)subplot(2,1,1)plot(F,20*log10(abs(H)))axis([04000-303])xlabel('Frequency(Hz)')ylabel('Magnitude(dB)')title('高通滤波器')gridonsubplot(2,1,2)pha=angle(H)*180/pi;plot(F,pha)gridon运行结果如下：3）采样速率为10000Hz，要求设计一个带通滤波器，pf=[1000Hz，1500Hz]，sf=[600Hz，1900Hz]，pR=3dB，sR=20dB。程序如下：fn=10000;fp=[1000,1500];fs=[600,1900];Rp=3;Rs=20;Wp=fp/(fn/2);%计算归一化角频率Ws=fs/(fn/2);[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数[H,F]=freqz(b,a,1000,10000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)subplot(2,1,1)plot(F,20*log10(abs(H)))axis([05000-303])xlabel('Frequency(Hz)')ylabel('Magnitude(dB)')title('带通滤波器')gridonsubplot(2,1,2)pha=angle(H)*180/pi;plot(F,pha)gridon运行结果如下：4）采样速率为10000Hz，要求设计一个带通滤波器，pf=[1000Hz，1500Hz]，sf=[1200Hz，1300Hz]，pR=3dB，sR=30dB。程序如下：fn=10000;fp=[1000,1500];fs=[1200,1300];Rp=3;Rs=30;Wp=fp/(fn/2);%计算归一化角频率Ws=fs/(fn/2);[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率[b,a]=butter(n,Wn,'stop');%计算H(z)分子、分母多项式系数[H,F]=freqz(b,a,1000,10000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)subplot(2,1,1)plot(F,20*log10(abs(H)))axis([05000-353])xlabel('Frequency(Hz)')ylabel('Magnitude(dB)')title('带阻滤波器')gridonsubplot(2,1,2)pha=angle(H)*180/pi;plot(F,pha)gridon运行结果如下：5）梳状滤波器的设计使用MATLAB中的函数“iircomb”可以设计出峰值或谷值滤波器H(z)的分子分母多项式系数，其格式为：[num,den]=iircomb(n,bw,ab,’type’)num,den分别为H(z)的分子、分母系数；n为梳状滤波器阶数，在数字归一化频率0～2pi区间梳状滤波器开槽数等于n+1；bw为滤波器开槽的abdB带宽，默认ab=-3dB；type可以是’notch’或’peak’，notch为开槽性梳状滤波器，peak为峰值性梳状滤波器。例：在采样频率为8000Hz的条件下设计一个在500Hz，1000Hz，2000Hz，……，n×500Hz的地方开槽的陷波，陷波带宽（-3dB处）为60Hz。程序如下：Fs=8000;%采样速率Ts=1/Fs;f0=500;%开槽基频率bw=60/(Fs/2);%归一化开槽带宽ab=-3;%开槽带宽位置处的衰减n=Fs/f0;%计算滤波器阶数[num,den]=iircomb(n,bw,ab,'notch');%计算H(z)分子分母多项式系数[H,F]=freqz(num,den,2000,8000);%计算滤波器的幅频响应subplot(2,1,1)plot(F,20*log10(abs(H)))axis([05000-353])xlabel('Frequency(Hz)')ylabel('Magnitude(dB)')title('梳状滤波器')gridonsubplot(2,1,2)pha=angle(H)*180/pi;plot(F,pha)gridon4．滤波器的实现：使用“filter”命令实现在以上设计中求出滤波器H(z)的分子分母系数向量[b,a]之后，可以用“filter”指令实现对输入信号进行滤波。作业：设计滤波器，把输入信号)800(c)400cos()200cos()(tostttx中的三个信号分离出来。附：低通滤波器对输入信号进行分离的程序及结果：Fs=2000;t=(1:100)/Fs;x=cos(2*pi*100*t)+cos(2*pi*200*t)+cos(2*pi*800*t);figure(1)subplot(2,1,1)plot(t,x)axis([00.05-44])X=fft(x);subplot(2,1,2)plot(abs(X))Wp=80/(Fs/2);Ws=120/(Fs/2);Rp=3;Rs=25;[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数[H,F]=freqz(b,a,512,Fs);%计算H(z)的幅频响应,freqz(b,a,计算点数，采样频率)figure(2)subplot(2,1,1)plot(F,20*log10(abs(H)+eps))ylabel('Magnitude(dB)')title('低通滤波器')axis([0,512,-30,3])gridony=filter(b,a,x);subplot(2,1,2)plot(t,y)