有限脉冲响应数字滤波器设计实验报告DOC.docx
《有限脉冲响应数字滤波器设计实验报告DOC.docx》由会员分享,可在线阅读,更多相关《有限脉冲响应数字滤波器设计实验报告DOC.docx(28页珍藏版)》请在冰豆网上搜索。
![有限脉冲响应数字滤波器设计实验报告DOC.docx](https://file1.bdocx.com/fileroot1/2023-2/27/95819485-e1f4-4f6d-8bb9-83c7d0556c95/95819485-e1f4-4f6d-8bb9-83c7d0556c951.gif)
有限脉冲响应数字滤波器设计实验报告DOC
成绩:
《数字信号处理》
作业与上机实验
(第二章)
班级:
学号:
姓名:
任课老师:
完成时间:
信息与通信工程学院
2014—2015学年第1学期
第7章有限脉冲响应数字滤波器设计
1、教材p238:
19.设信号x(t)=s(t)+v(t),其中v(t)是干扰,s(t)与v(t)的频谱不混叠,其幅度谱如题19图所示。
要求设计数字滤波器,将干扰滤除,指标是允许|s(f)|在0≤f≤15kHz频率范围中幅度失真为±2%(δ1=0.02);f>20kHz,衰减大于40dB(δ2=0.01);希望分别设计性价比最高的FIR和IIR两种滤波器进行滤除干扰。
请选择合适的滤波器类型和设计方法进行设计,最后比较两种滤波器的幅频特性、相频特性和阶数。
题19图
(1)matlab代码:
%基于双线性变换法直接设计IIR数字滤波器
Fs=80000;
fp=15000;fs=20000;rs=40;
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
Rp=-20*log10(1-0.02);As=40;
[N1,wp1]=ellipord(wp/pi,ws/pi,Rp,As);
[B,A]=ellip(N1,Rp,As,wp1);
[Hk,wk1]=freqz(B,A,1000);
mag=abs(Hk);pah=angle(Hk);
%窗函数法设计FIR数字滤波器
Bt=ws-wp;
alph=0.5842*(rs-21)^0.4+0.07886*(rs-21);
N=ceil((rs-8)/2.285/Bt);
wc=(wp+ws)/2/pi;
hn=fir1(N,wc,kaiser(N+1,alph));
M=1024;
Hk=fft(hn,M);
k=0:
M/2-1;
wk=(2*pi/M)*k;
%画出各种比较结果图
figure
(2);
plot(wk/pi,20*log10(abs(Hk(k+1))),':
','linewidth',2.5);
holdon
plot(wk1/pi,20*log10(mag),'linewidth',2);
holdoff
legend('FIR滤波器','IIR滤波器');
axis([0,1,-80,5]);xlabel('w/\pi');ylabel('幅度/dB');
title('损耗函数');
figure(3)
plot(wk/pi,angle(Hk(k+1))/pi,':
','linewidth',2.5);
holdon
plot(wk1/pi,pah/pi,'linewidth',2);
holdoff
legend('FIR滤波器','IIR滤波器');
xlabel('w/\pi');ylabel('相位/\pi');
title('相频特性曲线');
(2)两种数字滤波器的损耗函数和相频特性的比较分别如图1、2所示:
图1损耗函数比较图图2相频特性比较图
(3)IIR数字滤波器阶数:
N=5
FIR数字滤波器阶数:
N=36
(4)运行结果分析:
由图2及阶数可见,IIR阶数低得多,但相位特性存在非线性失真,FIR具有线性相位特性。
20.调用MATLAB工具箱函数fir1设计线性相位低通FIR滤波器,要求希望逼近的理想低通滤波器通带截止频率ωc=π/4rad,滤波器长度N=21。
分别选用矩形窗、Hanning窗、Hamming窗和Blackman窗进行设计,绘制用每种窗函数设计的单位脉冲响应h(n)及其损耗函数曲线,并进行比较,观察各种窗函数的设计性能。
(1)matlab代码:
wc=pi/4;N=21;
hn_boxcar=fir1(N-1,wc/pi,boxcar(N));
hn_hanning=fir1(N-1,wc/pi,hanning(N));
hn_hamming=fir1(N-1,wc/pi,hamming(N));
hn_blackman=fir1(N-1,wc/pi,blackman(N));
n=0:
N-1;
plot(n,hn_boxcar);
holdon
plot(n,hn_hanning,':
','linewidth',2);
plot(n,hn_hamming,'+','linewidth',2);
plot(n,hn_blackman,'o');
holdoff
xlabel('n');ylabel('h(n)');
legend('矩形窗','汉宁窗','哈明窗','布莱克曼窗');
title('单位冲激响应');
M=1024;Hk=fft(hn_boxcar,M);k=0:
M/2-1;
wk=(2*pi/M)*k;
figure();
plot(wk/pi,20*log10(abs(Hk(k+1))),'linewidth',2);
Hk=fft(hn_hanning,M)
holdon
plot(wk/pi,20*log10(abs(Hk(k+1))),':
','linewidth',3);
Hk=fft(hn_hamming,M)
plot(wk/pi,20*log10(abs(Hk(k+1))),'o');
Hk=fft(hn_blackman,M)
plot(wk/pi,20*log10(abs(Hk(k+1))),'*');
holdoff
legend('矩形窗','汉宁窗','哈明窗','布莱克曼窗');
axis([0,1,-80,5]);xlabel('w/\pi');ylabel('幅度/dB');
title('损耗函数');
(2)四种窗函数设计的单位脉冲响应的比较如图3所示:
图3单位脉冲响应比较图
(3)四种窗函数设计的损耗函数的比较如图4所示:
图4损耗函数比较图
(4)运行结果分析:
由图4可见,当滤波器长度N不变时,矩形窗设计的滤波器的过渡带最窄,阻带最小衰减最小;布莱克曼窗设计的滤波器的过渡带最宽,同时阻带最小衰减最大。
21.将要求改成设计线性相位高通FIR滤波器,重作题20。
(1)matlab代码:
wc=pi/4;N=21;
hn_boxcar=fir1(N-1,wc/pi,'high',boxcar(N));
hn_hanning=fir1(N-1,wc/pi,'high',hanning(N));
hn_hamming=fir1(N-1,wc/pi,'high',hamming(N));
hn_blackman=fir1(N-1,wc/pi,'high',blackman(N));
n=0:
N-1;
plot(n,hn_boxcar);
holdon
plot(n,hn_hanning,':
','linewidth',2);
plot(n,hn_hamming,'+','linewidth',2);
plot(n,hn_blackman,'o');
holdoff
xlabel('n');ylabel('h(n)');
legend('矩形窗','汉宁窗','哈明窗','布莱克曼窗');
title('单位冲激响应');
M=1024;Hk=fft(hn_boxcar,M);k=0:
M/2-1;
wk=(2*pi/M)*k;
figure();
plot(wk/pi,20*log10(abs(Hk(k+1))),'linewidth',2);
Hk=fft(hn_hanning,M)
holdon
plot(wk/pi,20*log10(abs(Hk(k+1))),':
','linewidth',3);
Hk=fft(hn_hamming,M)
plot(wk/pi,20*log10(abs(Hk(k+1))),'o');
Hk=fft(hn_blackman,M)
plot(wk/pi,20*log10(abs(Hk(k+1))),'*');
holdoff
legend('矩形窗','汉宁窗','哈明窗','布莱克曼窗');
axis([0,1,-80,5]);xlabel('w/\pi');ylabel('幅度/dB');
title('损耗函数');
(2)四种窗函数设计的单位脉冲响应的比较如图5所示:
图5单位脉冲响应比较图
(3)四种窗函数设计的损耗函数的比较如图6所示:
图6损耗函数比较图
(5)运行结果分析:
由图6可见,当滤波器长度N不变时,矩形窗设计的滤波器的过渡带最窄,阻带最小衰减最小;布莱克曼窗设计的滤波器的过渡带最宽,同时阻带最小衰减最大。
25.调用MATLAB工具箱函数fir1设计线性相位高通FIR滤波器。
要求通带截止频率为0.6πrad,阻带截止频率为0.45π,通带最大衰减为0.2dB,阻带最小衰减为45dB。
显示所设计的单位脉冲响应h(n)的数据,并画出损耗函数曲线。
(1)matlab代码:
wp=0.6*pi;ws=0.45*pi;
Bt=wp-ws;
N0=ceil(6.6*pi/Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
hn=fir1(N-1,wc,'high',hamming(N));
M=1024;
Hk=fft(hn,M);
n=0:
N-1;
stem(n,hn);
xlabel('n');ylabel('h(n)');
title('单位冲激响应');
k=0:
M/2-1;
wk=(2*pi/M)*k;
figure
(2);
plot(wk/pi,20*log10(abs(Hk(k+1))));
axis([0,1,-80,5]);xlabel('w/\pi');ylabel('幅度/dB');
title('损耗函数');gridon
(2)高通FIR滤波器的单位脉冲响应、损耗函数如图7、8所示:
图7单位脉冲响应图8损耗函数
26.调用MATLAB工具箱函数fir1设计线性相位带通FIR滤波器。
要求通带截止频率为0.55πrad和0.7πrad,阻带截止频率为0.45πrad和0.8πrad,通带最大衰减为0.15dB,阻带最小衰减为40dB。
显示所设计的单位脉冲响应h(n)的数据,并画出损耗函数曲线。
(1)matlab代码:
wp1=0.55*pi;wp2=0.7*pi;ws1=0.45*pi;ws2=0.8*pi;
Bt=wp2-wp1;
N=ceil(6.2*pi/Bt);
wc=[(wp1+ws1)/2/pi,(ws2+wp2)/2/pi];
hn=fir1(N-1,wc,hanning(N));
M=1024;
Hk=fft(hn,M);
n=0:
N-1;
stem(n,hn);
xlabel('n');ylabel('h(n)');
title('单位冲激响应');
k=0:
M/2-1;
wk=(2*pi/M)*k;
figure
(2);
plot(wk/pi,20*log10(abs(Hk(k+1))));
axis([0,1,-80,5]);xlabel('w/\pi');ylabel('幅度/dB');
title('损耗函数');gridon
(2)带通FIR滤波器的单位脉冲响应、损耗函数如图9、10所示:
图9单位脉冲响应图10损耗函数
2、某信号
为:
,其中
设计最低阶FIR数字滤波器,按下图所示对
进行数字滤波处理,实现:
1)将
频率分量以高于50dB的衰减抑制,同时以低于2dB的衰减通过
和
频率分量;
一、基于窗函数法设计FIR数字滤波器:
(1)matlab代码:
Fs=3800;
fp=130;
fs=600;
rs=50;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Bt=ws-wp;
alph=0.5842*(rs-21)^0.4+0.07886*(rs-21);
N=ceil((rs-8)/2.285/Bt);
wc=(wp+ws)/2/pi;
hn=fir1(N,wc,kaiser(N+1,alph));
M=1024;
Hk=fft(hn,M);
k=0:
M/2-1;
wk=(2*pi/M)*k;
figure
(2);
plot(wk/pi,20*log10(abs(Hk(k+1))));
axis([0,1,-80,5]);
xlabel('w/\pi');
ylabel('幅度/dB');
title('损耗函数');gridon
figure(3)
plot(wk/pi,angle(Hk(k+1))/pi);gridon
xlabel('w/\pi');
ylabel('相位/\pi');
title('相频特性曲线');
(2)数字滤波器的损耗函数和相频特性分别如图11、12所示:
图11损耗函数曲线图12相频特性曲线
二、按直接型网络结构编程编写滤波程序:
(1)matlab代码:
N=500;n=0:
N-1;
f=2800;T=1/f;t=n*T;
x=0.5*cos(2*pi*100*t)+0.7*cos(2*pi*130*t+0.1*pi)+0.4*cos(2*pi*600*t);
m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;m7=0;m8=0;m9=0;m10=0;m11=0;m12=0;m13=0;m14=0;
m15=0;m16=0;m17=0;m18=0;m19=0;m20=0;m21=0;m22=0;m23=0;m24=0;
form=1:
length(x)
y(m)=0.0012*x(m)+m1*0.0011-m2*0.0014-m3*0.0072-m4*0.0147-m5*0.0193-...
m6*0.0145+m7*0.0055+m8*0.0423+m9*0.0910+m10*0.1410+m11*0.1786+m12*0.1926+...
m13*0.1786+m14*0.1410+m15*0.0910+m16*0.0423+m17*0.0055-m18*0.0145-...
m19*0.0193-m20*0.0147-m21*0.0072-m22*0.0014+m23*0.0011+m24*0.0012;
m24=m23;m23=m22;m22=m21;m21=m20;m20=m19;m19=m18;m18=m17;m17=m16;m16=m15;
m15=m14;m14=m13;m13=m12;m12=m11;m11=m10;m10=m9;m9=m8;m8=m7;m7=m6;m6=m5;
m5=m4;m4=m3;m3=m2;m2=m1;m1=x(m);
end
plot(n,x);
title('信号x(n)');ylabel('幅值');xlabel('n');
S=fft(x,N);fs=n/(N*T);
figure
(2)
plot(fs,abs(S));
axis([0,1500,0,180]);
title('原信号幅度频谱(采样点数为500)');
xlabel('频率/Hz');ylabel('幅值');
figure(3)
plot(n,y);title('信号y(n)');
ylabel('幅值');xlabel('n');
S=fft(y,N);fs=n/(N*T);
figure(4)
plot(fs,abs(S));
axis([0,1500,0,160]);
title('幅度频谱');xlabel('频率/Hz');ylabel('幅值');
(2)
原信号及其幅度频谱分别如图13、14所示:
图13信号x(n)波形图14幅度频谱
(3)滤波后信号y(n)及其幅度频谱分别如图15、16所示:
图15信号y(n)波形图16幅度频谱
2)将
和
频率分量以高于50dB的衰减抑制,同时以低于2dB的衰减通过
频率分量;
一、基于频率采样法设计FIR数字滤波器:
(1)matlab代码:
T=0.48;
Fs=3800;
fp=600;fs=100;
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
datB=wp-ws;wc=wp;
m=1;
N=ceil((m+1)*2*pi/datB+1);
N=N+mod(N+1,2);
Np=fix(wc/(2*pi/N));Ns=N-2*Np-1;
Ak=[zeros(1,Np+1),ones(1,Ns),zeros(1,Np)];
Ak(Np+2)=T;Ak(N-Np)=T;
thetak=-pi*(N-1)*(0:
N-1)/N;
Hk=Ak.*exp(1j*thetak);
hn=real(ifft(Hk));
M=1024;
Hk=fft(hn,M);
k=0:
M/2-1;
wk=(2*pi/M)*k;
figure
(2);
plot(wk/pi,20*log10(abs(Hk(k+1))));
axis([0,1,-80,5]);xlabel('w/\pi');ylabel('幅度/dB');
title('损耗函数');gridon
figure(3)
plot(wk/pi,angle(Hk(k+1))/pi);gridon
xlabel('w/\pi');ylabel('相位/\pi');
title('相频特性曲线');
(2)数字滤波器的损耗函数和相频特性分别如图17、18所示:
图17损耗函数曲线图18相频特性曲线
二、按直接型网络结构编程编写滤波程序:
(1)matlab代码:
N=500;n=0:
N-1;
f=2800;T=1/f;t=n*T;
x=0.5*cos(2*pi*100*t)+0.7*cos(2*pi*130*t+0.1*pi)+0.4*cos(2*pi*600*t);
m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;m7=0;m8=0;m9=0;m10=0;m11=0;m12=0;m13=0;m14=0;
m15=0;m16=0;
form=1:
length(x)
y(m)=-0.0009*x(m)-m1*0.0169-m2*0.0128+m3*0.0282+m4*0.0627+m5*0.0198-...
m6*0.1198-m7*0.2827+m8*0.6447-m9*0.2827-m10*0.1198+m11*0.0198+m12*0.0627+...
m13*0.0282-m14*0.0128-m15*0.0169-m16*0.0009;
m16=m15;m15=m14;m14=m13;m13=m12;m12=m11;m11=m10;m10=m9;m9=m8;m8=m7;
m7=m6;m6=m5;m5=m4;m4=m3;m3=m2;m2=m1;m1=x(m);
end
plot(n,x);
title('信号x(n)');ylabel('幅值');xlabel('n');
S=fft(x,N);fs=n/(N*T);
figure
(2)
plot(fs,abs(S));
axis([0,1500,0,180]);
title('原信号幅度频谱(采样点数为500)');
xlabel('频率/Hz');ylabel('幅值');
figure(3)
plot(n,y);
title('信号y(n)');
ylabel('幅值');xlabel('n');
S=fft(y,N);fs=n/(N*T);
figure(4)
plot(fs,abs(S));
axis([0,1500,0,160]);
title('幅度频谱');
xlabel('频率/Hz');ylabel('幅值');
(4)
原信号及其幅度频谱分别如图19、20所示:
图19信号x(n)波形图20幅度频谱
(5)
滤波后信号y(n)及其幅度频谱分别如图21、22所示:
图21信号y(n)波形图22幅度频谱
要求:
按数字滤波器直接型结构图编写滤波程序,求得
;1)中的FIR滤波器采用窗函数法设计;2)中的FIR滤波器采用频率采样法设计。
画出所设计的滤波器频率特性图、信号时域图;给出滤波器设计的MATLAB代码与滤波器实现的代码;选择合适的信号采样周期T。
3)与第6章作业2的IIR滤波方法进行比较研究。
一、低通滤波器部分:
(1)matlab代码:
Fs=3800;
fp=130;fs=600;rs=50;
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
Bt=ws-wp;
alph=0.5842*(rs-21)^0.4+0.07886*(rs-21);
N=ceil((rs-8)/2.285/Bt);
wc=(wp+ws)/2/pi;
hn=fir1(N,wc,kaiser(N+1,alph));
M=1024;
Hk=fft(hn,M);
k=0:
M/2-1;
wk=(2*pi/M)*k;
wp2=2*fp/Fs;ws2=2*fs/Fs;
Rp=2;As=50;
[N1,wp1]=ellipord(wp2,ws2,Rp,As);
[B,A]=ellip(N1,Rp,As,wp1)
[Hk1,wk1]=freqz(B,A);
mag=abs(Hk1);pah=angle(Hk1);
plot(wk1/pi,20*log10(mag));gridon
holdon
plot(wk/pi,20*log10(abs(Hk(k+1))),':
','linewidth',3);gridon
holdoff
xlabel('w/\pi');ylabel('幅度/dB')
title('损耗函数曲线');
legend('IIR','FIR');
figure
(2)
plot(wk1/pi,pah/pi);gridon
holdon
plot(wk/pi,angle(Hk(k+1))/pi,':
','linewidth',3);gridon
holdoff
xlabel('w/\pi');ylabel('相位/\pi');
title('相频特性曲线');
legend('IIR','FIR');
(2)两种滤波器的损耗函数、相频特性的比较图见图23、24:
图23损耗函数比较图图24相频特性比较图
(3)IIR滤波器的阶数:
N1=3
FIR滤波器的阶数:
N=17
二、高通滤波器部分:
(1)matlab代码:
T=0.48;
Fs=3800;
fp=600;fs=100;
wp=2*pi*fp/Fs;ws=2*pi*fs/