有限脉冲响应数字滤波器设计实验报告DOC.docx
《有限脉冲响应数字滤波器设计实验报告DOC.docx》由会员分享，可在线阅读，更多相关《有限脉冲响应数字滤波器设计实验报告DOC.docx（28页珍藏版）》请在冰豆网上搜索。
有限脉冲响应数字滤波器设计实验报告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/
升级会员 