数字滤波器设计及应用综合实验样本.doc

资料内容仅供您学习参考，如有不当或者侵权，请联系改正或者删除。

数字信号处理实验
实验四、 数字滤波器设计及
应用综合实验

学院: 信息工程学院
班级: 电子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('\n Butterworth Filter Order=%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('digital filter Magnitude Response')
axis([0,1,0,1.1])
subplot(2,1,2);plot(w/pi,db);
title('digital filter Magnitude in DB')
axis([0,1,-40,5]);
结果:
N1 =
5.8858
N =
6
Butterworth Filter Order= 6
OmegaC =
0.7032
b =
0.0000 0.0006 0.0101 0.0161 0.0041 0.0001
a =
1.0000 -3.3635 5.0684 -4.2759 2.1066 -0.5706 0.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('\n Chebyshev Filter Order=%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('digital filter Magnitude Response'); axis([0,1,0,1.1])
subplot(2,1,2);plot(w/pi,db);
title('digital filter Magnitude in DB');
axis([0,1,-40,5]);
Chebyshev Filter Order= 4
b =
0.0018 0.0073 0.0110 0.0073 0.0018
a =
1.0000 -3.0543 3.8290 -2.2925 0.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);
if N>M
b=[b zeros(1,N-M)];
elseif M>N
a=[a zeros(1,M-N)];
else
NM=0;
end
K=floor(N/2); B=zeros(K,3); A=zeros(K,3);
if K*2==N
b=[b 0]; a=[a 0];
end
broots=cplxpair(roots(b));
aroots=cplxpair(roots(a));
for i=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('digital filter Magnitude Response');
axis([0,1,0,1.1])
subplot(2,1,2);plot(w/pi,db);
title('digital filter Magnitude in DB');
axis([0,1,-40,5]);
[b0,B,A]=dir2cas(b,a)
结果:
N1 =
5.3044
N =
6
Butterworth Filter Order= 6
OmegaC =
0.7273
b =
0.0006 0.0035 0.0087 0.0116 0.0087 0.0035 0.0006
a =
1.0000 -3.3143 4.9501 -4.1433 2.0275 -0.5458 0.0628
B =
1.0000 2.0335 1.0338
1.0000 1.9996 1.0000
1.0000 1.9669 0.9673
A =
1.0000 -0.9459 0.2342
1.0000 -1.0541 0.3753
1.0000 -1.3143 0.7149
本实验波形如下:

2. FIR数字滤波器设计
本实验所需M文件如下:
function hd=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.2 0.5 -2 2]); subplot(212)
Fsf=fft(sf,512); AFsf=abs(Fsf);
f=(0:255)*fs/512; plot(f,AFsf(1:256))
xlabel('频率/赫兹');
ylabel('幅度');
title('信号滤波后频域图');
本实验波形如图:

3.利用数字滤波器实现DTMF信号的提取
双音多频( Dual Tone Multi Frequency) 信号是音频电话中的拨号信号, 每一个数字( 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([10 50 0 120])
本实验波形如图:

自己设计实验:
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;
for m=1:6;
d=fix(TN/10^(6-m));
TN=TN-d*10^(6-m);
for p=1:4;
for q=1:4;
if tm(p,q)==abs(d); break,end
end
if tm(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([10 50 0 120])
limit = 80;
for s=5:8;
if val(s) > limit, break, end
end
for r=1:4;
if val(r) > limit, break, end
end
TNr=TNr+tm(r,s-4)*10^(6-m);
end
disp('接收端检测到的号码为: ')
disp(TNr)
实验波形:

三、 实验心得: