数字滤波器设计及应用综合实验样本.docx
《数字滤波器设计及应用综合实验样本.docx》由会员分享，可在线阅读，更多相关《数字滤波器设计及应用综合实验样本.docx（16页珍藏版）》请在冰豆网上搜索。
数字滤波器设计及应用综合实验样本
数字信号处理实验
实验四、数字滤波器设计及
应用综合实验
学院:
信息工程学院
班级:
电子101班
姓名:
学号:
一、实验目的
1．熟悉IIR数字滤波器的设计原理及方法。
2．熟悉FIR数字滤波器的设计原理及方法。
3.掌握利用Matlab实现数字滤波器的方法
4.掌握利用数字滤波器进行信号处理的方法。
5.了解基于Simulink的动态仿真实现信号滤波的基本方法。
二、实验内容及要求
实验内容:
综合运用数字滤波器设计的相关知识,根据给定设计方法要求,用脉冲响应不变法和双线性变换法设计IIR数字滤波器;利用窗函数设计法设计FIR数字滤波器。
根据实际信号的频谱特性,分析、确定滤波器设计技术指标,实现对信号的滤波。
1．IIR数字滤波器设计
（1）用脉冲响应不变法设计巴特沃斯数字滤波器。
该实验所需M文件如下:
①、butterworth低通滤波器原型设计函数:
function[b,a]=afd_butt（Wp,Ws,Rp,As）
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（Wp/Ws）））
fprintf（'\nButterworthFilterOrder=%2.0f\n',N）
OmegaC=Wp/（（10^（Rp/10）-1）^（1/（2*N）））
[b,a]=u_buttap（N,OmegaC）
②、非归一化Butterworth模拟低通滤波器设计函数:
function[b,a]=u_buttap（N,Omegac）;
[z,p,k]=buttap（N）;
p=p*Omegac;
k=k*Omegac^N;
B=real（poly（z））;
b=k*B;
a=real（poly（p））;
③、利用脉冲响应不变法从模拟到数字滤波器变换函数:
function[b,a]=imp_invr（c,d,T）
[R,p,k]=residue（c,d）;
p=exp（p*T）;
[b,a]=residuez（R,p,k）;
b=real（b'）;
a=real（a'）;
④、频率响应函数freqz的修正:
function[db,mag,pha,w]=freqz_m（b,a）;
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
501））';
w=（w（1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
本实验程序如下:
wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;
OmegaP=wp/T;OmegaS=ws/T;
[cs,ds]=afd_butt（OmegaP,OmegaS,Rp,As）;
[b,a]=imp_invr（cs,ds,T）
[db,mag,pha,w]=freqz_m（b,a）;
subplot（2,1,1）;plot（w/pi,mag）;
title（'digitalfilterMagnitudeResponse'）
axis（[0,1,0,1.1]）
subplot（2,1,2）;plot（w/pi,db）;
title（'digitalfilterMagnitudeinDB'）
axis（[0,1,-40,5]）;
结果:
N1=
5.8858
N=
6
ButterworthFilterOrder=6
OmegaC=
0.7032
b=
0.00000.00060.01010.01610.00410.0001
a=
1.0000-3.36355.0684-4.27592.1066-0.57060.0661
本实验波形图如下:
（2）用双线性变换法设计切比雪夫数字滤波器。
本实验所需M文件如下:
非归一化切比雪夫I型模拟低通滤波器原型设计:
function[b,a]=u_chb1ap（N,Rp,Omegac）;
[z,p,k]=cheb1ap（N,Rp）;
a=real（poly（p））;
aNn=a（N+1）;
p=p*Omegac;
a=real（poly（p））;
aNu=a（N+1）;
k=k*aNu/aNn;
B=real（poly（z））;
b=k*B;
本实验程序如下:
wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;
OmegaP=（2/T）*tan（wp/2）;
OmegaS=（2/T）*tan（ws/2）;
ep=sqrt（10^（Rp/10）-1）;
Ripple=sqrt（1/（1+ep*ep））;
Attn=1/（10^（As/20））;
A1=1/Attn;a1=sqrt（A1*A1-1）/ep;
a2=OmegaS/OmegaP;
N=ceil（logm（a1+sqrt（a1*a1-1））/logm（a2+sqrt（a2*a2-1）））;
fprintf（'\nChebyshevFilterOrder=%2.0f\n',N）
[cs,ds]=u_chb1ap（N,Rp,OmegaP）;
[b,a]=bilinear（cs,ds,1/T）
[db,mag,pha,w]=freqz_m（b,a）;
subplot（2,1,1）;plot（w/pi,mag）;
title（'digitalfilterMagnitudeResponse'）;axis（[0,1,0,1.1]）
subplot（2,1,2）;plot（w/pi,db）;
title（'digitalfilterMagnitudeinDB'）;
axis（[0,1,-40,5]）;
ChebyshevFilterOrder=4
b=
0.00180.00730.01100.00730.0018
a=
1.0000-3.05433.8290-2.29250.5507
本实验波形如下:
（3）用双线性变换法设计巴特沃斯数字滤波器,并将直接型结构转换成级联型结构。
将直接形式变为级联形式函数:
function[b0,B,A]=dir2cas（b,a）
b0=b
（1）;b=b/b0;a0=a
（1）;a=a/a0;b0=b0/a0;
M=length（b）;N=length（a）;
ifN>M
b=[bzeros（1,N-M）];
elseifM>N
a=[azeros（1,M-N）];
else
NM=0;
end
K=floor（N/2）;B=zeros（K,3）;A=zeros（K,3）;
ifK*2==N
b=[b0];a=[a0];
end
broots=cplxpair（roots（b））;
aroots=cplxpair（roots（a））;
fori=1:
2:
2*K
Brow=broots（i:
1:
i+1,:
）;
Brow=real（poly（Brow））;
B（fix（i+1）/2,:
）=Brow;
Arow=aroots（i:
1:
i+1,:
）;
Arow=real（poly（Arow））;
A（fix（i+1）/2,:
）=Arow;
end
本实验程序如下:
wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;
OmegaP=（2/T）*tan（wp/2）;
OmegaS=（2/T）*tan（ws/2）;
[cs,ds]=afd_butt（OmegaP,OmegaS,Rp,As）;
[b,a]=bilinear（cs,ds,T）
[db,mag,pha,w]=freqz_m（b,a）;
subplot（2,1,1）;plot（w/pi,mag）;
title（'digitalfilterMagnitudeResponse'）;
axis（[0,1,0,1.1]）
subplot（2,1,2）;plot（w/pi,db）;
title（'digitalfilterMagnitudeinDB'）;
axis（[0,1,-40,5]）;
[b0,B,A]=dir2cas（b,a）
结果:
N1=
5.3044
N=
6
ButterworthFilterOrder=6
OmegaC=
0.7273
b=
0.00060.00350.00870.01160.00870.00350.0006
a=
1.0000-3.31434.9501-4.14332.0275-0.54580.0628
B=
1.00002.03351.0338
1.00001.99961.0000
1.00001.96690.9673
A=
1.0000-0.94590.2342
1.0000-1.05410.3753
1.0000-1.31430.7149
本实验波形如下:
2.FIR数字滤波器设计
本实验所需M文件如下:
functionhd=ideal_lp（wc,M）;
alpha=（M-1）/2;
n=[0:
（M-1）];
m=n-alpha+eps;
hd=sin（wc*m）./（pi*m）;
本实验程序如下:
f1=100;f2=200;
fs=;
m=（0.3*f1）/（fs/2）;
M=round（8/m）;
N=M-1;
b=fir1（N,0.5*f2/（fs/2））;
figure
（1）
[h,f]=freqz（b,1,512）;
%[H,W]=freqz（B,A,N）
plot（f*fs/（2*pi）,20*log10（abs（h）））
xlabel（'频率/赫兹'）;
ylabel（'增益/分贝'）;
title（'滤波器的增益响应'）;
figure
（2）
subplot（211）
t=0:
1/fs:
0.5;
s=sin（2*pi*f1*t）+sin（2*pi*f2*t）;
plot（t,s）;
xlabel（'时间/秒'）;
ylabel（'幅度'）;
title（'信号滤波前时域图'）;
subplot（212）
Fs=fft（s,512）;
AFs=abs（Fs）;
f=（0:
255）*fs/512;
plot（f,AFs（1:
256））;
xlabel（'频率/赫兹'）;ylabel（'幅度'）;
title（'信号滤波前频域图'）;
figure（3）
sf=filter（b,1,s）;subplot（211）
plot（t,sf）xlabel（'时间/秒'）;
ylabel（'幅度'）;
title（'信号滤波后时域图'）;
axis（[0.20.5-22]）;subplot（212）
Fsf=fft（sf,512）;AFsf=abs（Fsf）;
f=（0:
255）*fs/512;plot（f,AFsf（1:
256））
xlabel（'频率/赫兹'）;
ylabel（'幅度'）;
title（'信号滤波后频域图'）;
本实验波形如图:
3.利用数字滤波器实现DTMF信号的提取
双音多频（DualToneMultiFrequency）信号是音频电话中的拨号信号,每一个数字（0-9）由两个不同频率单音组成（每个单音用正弦信号表示）,所用频率分为高频带和低频带两组,每个数字由高、低频带中各一个频率组成,例如数字9使用852Hz和1477Hz两个频率。
数字与符号对应频率关系见表3-1所示。
表3-1双频拨号的频率分配
1209Hz
1336Hz
1477Hz
1633Hz
697Hz
1
2
3
A
770Hz
4
5
6
B
852Hz
7
8
9
C
941Hz
*
0
#
D
实验内容:
电话中DTMF信号的产生于检测方法:
在电话中,数字0~9的中每一个都用两个不同的单音频传输,所用的8个频率分成高频带和低频带两组,低频带有四个频率:
679Hz,770Hz,852Hz和941Hz;高频带也有四个频率:
1209Hz,1336Hz,1477Hz和1633Hz.。
每一个数字均由高、低频带中各一个频率构成,例如1用697Hz和1209Hz两个频率,信号用
表示,其中
。
这样8个频率形成16种不同的双频信号。
具体号码以及符号对应的频率如表3-2所示。
表中最后一列在电话中暂时未用。
列
行
1209Hz
1336Hz
1477Hz
633Hz
697Hz
1
2
3
A
770Hz
4
5
6
B
852Hz
7
8
9
C
942Hz
*
0
#
D
表3-2双频拨号的频率分配
DTMF信号在电话中有两种作用,一个是用拨号信号去控制交换机接通被叫的用户电话机,另一个作用是控制电话机的各种动作,如播放留言、语音信箱等。
自己选择一个数字（0-9）,经过数字方法产生该数字的双频信号;设采样频率为8000Hz。
选择数字1,产生双频信号:
N=205;K=[18,20,22,24,31,34,38,42];
f1=[697,770,852,941];
f2=[1209,1336,1477,1633];
n=0:
1023;
x=sin（2*pi*n*f1
（1）/8000）+sin（2*pi*n*f2
（1）/8000）;
X=goertzel（x（1:
N）,K+1）;
val=abs（X）;
subplot
stem（K,val,'.'）;grid;xlabel（'k'）;ylabel（'|X（k）|'）;
axis（[10500120]）
本实验波形如图:
自己设计实验:
DTMF双频拨号信号的生成和检测程序
tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68];
N=205;K=[18,20,22,24,31,34,38,42];
f1=[697,770,852,941];
f2=[1209,1336,1477,1633];
TN=input（'键入6位电话号码='）;
TNr=0;
form=1:
6;
d=fix（TN/10^（6-m））;
TN=TN-d*10^（6-m）;
forp=1:
4;
forq=1:
4;
iftm（p,q）==abs（d）;break,end
end
iftm（p,q）==abs（d）;break,end
end
n=0:
1023;
x=sin（2*pi*n*f1（p）/8000）+sin（2*pi*n*f2（q）/8000）;
sound（x,8000）;
pause（0.1）
%接收检测端的程序
X=goertzel（x（1:
N）,K+1）;
val=abs（X）;
subplot（3,2,m）;
stem（K,val,'.'）;grid;xlabel（'k'）;ylabel（'|X（k）|'）
axis（[10500120]）
limit=80;
fors=5:
8;
ifval（s）>limit,break,end
end
forr=1:
4;
ifval（r）>limit,break,end
end
TNr=TNr+tm（r,s-4）*10^（6-m）;
end
disp（'接收端检测到的号码为:
'）
disp（TNr）
实验波形:
三、实验心得:
升级会员 