数字滤波器设计及应用综合实验.docx

数字滤波器设计及应用综合实验.docx

《数字滤波器设计及应用综合实验.docx》由会员分享，可在线阅读，更多相关《数字滤波器设计及应用综合实验.docx（17页珍藏版）》请在冰豆网上搜索。

数字滤波器设计及应用综合实验

数字信号处理实验

实验四、数字滤波器设计及

应用综合实验

学院：

信息工程学院

班级：

电子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=2000;

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）

实验波形：

三、实验心得：