实验1 离散系统的时域分析.docx

实验1 离散系统的时域分析.docx

《实验1 离散系统的时域分析.docx》由会员分享，可在线阅读，更多相关《实验1 离散系统的时域分析.docx（26页珍藏版）》请在冰豆网上搜索。

实验1离散系统的时域分析

实验1离散系统的时域分析

一、实验目的：

加深对离散系统的差分方程、单位抽样响应和卷积分析方法的理解。

二、实验原理：

离散系统

其输入、输出关系可用以下差分方程描述：

输入信号分解为冲激信号，

系统单位抽样序列h（n），

则系统响应为如下的卷积计算式：

当

时，h（n）是有限长度的（n：

[0，M]），称系统为FIR系统；反之，称系统为IIR系统。

三、实验内容

编制程序求解下列两个系统的单位抽样响应，并绘出其图形。

（1）

（2）

（1）源程序：

N=21;

b=[1-1];

a=[10.750.125];

x=[1zeros（1,N-1）];

n=0:

1:

N-1;

y=filter（b,a,x）;

stem（n,y）;

xlabel（'n'）;ylabel（'幅度'）;title（'单位抽样响应'）;

图形：

（2）源程序：

N=20;

d=[00.250.250.250.25];

c=[1];

x=[1zeros（1,N-1）];

n=0:

1:

N-1;

y=filter（d,c,x）;

stem（n,y）;

xlabel（'n'）;ylabel（'幅度'）;title（'单位抽样响应'）;

图形如下：

实验2离散系统的频率响应分析和零、极点分布

一、实验目的：

加深对离散系统的频率响应分析和零、极点分布的概念理解。

二、实验原理：

离散系统的时域方程为

其变换域分析方法如下：

系统函数为

分解因式

，

其中

和

称为零、极点。

在MATLAB中，可以用函数[z，p，K]=tf2zp（num，den）求得有理分式形式的系统函数的零、极点，用函数zplane（z，p）绘出零、极点分布图；也可以用函数zplane（num，den）直接绘出有理分式形式的系统函数的零、极点分布图。

使h=freqz（num,den,w）函数可求系统的频率响应，w是频率的计算点，如w=0:

pi/255:

pi,h是复数，abs（h）为幅度响应，angle（h）为相位响应。

另外，在MATLAB中，可以用函数[r，p，k]=residuez（num，den）完成部分分式展开计算；可以用函数sos=zp2sos（z，p，K）完成将高阶系统分解为2阶系统的串联。

三、实验内容

练习1求下列直接型系统函数的零、极点，并将它转换成二阶节形式

解用MATLAB计算程序如下：

num=[1-0.1-0.3-0.3-0.2];

den=[10.10.20.20.5];

[z,p,k]=tf2zp（num,den）;

disp（'零点'）;disp（z）;

disp（'极点'）;disp（p）;

disp（'增益系数'）;disp（k）;

sos=zp2sos（z,p,k）;

disp（'二阶节'）;disp（real（sos））;

zplane（num,den）

输入到“num”和“den”的分别为分子和分母多项式的系数。

计算求得零、极点增益系数和二阶节的系数：

零点

0.9615

-0.5730

-0.1443+0.5850i

-0.1443-0.5850i

极点

0.5276+0.6997i

0.5276-0.6997i

-0.5776+0.5635i

-0.5776-0.5635i 

增益系数

1 

二阶节

1.0000-0.3885-0.55091.00001.15520.6511

1.00000.28850.36301.0000-1.05520.7679 

系统函数的二阶节形式为：

极点图如右图。

2差分方程

所对应的系统的频率响应。

解：

差分方程所对应的系统函数为

用MATLAB计算的程序如下：

k=256;

num=[0.8-0.440.360.02];

den=[10.7-0.45-0.6];

w=0:

pi/k:

pi;

h=freqz（num,den,w）;

subplot（2,2,1）;

plot（w/pi,real（h））;grid

title（'实部'）

xlabel（'\omega/\pi'）;ylabel（'幅度'）

subplot（2,2,2）;

plot（w/pi,imag（h））;grid

title（'虚部'）

xlabel（'\omega/\pi'）;ylabel（'Amplitude'）

subplot（2,2,3）;

plot（w/pi,abs（h））;grid

title（'幅度谱'）

xlabel（'\omega/\pi'）;ylabel（'幅值'）

subplot（2,2,4）;

plot（w/pi,angle（h））;grid

title（'相位谱'）

xlabel（'\omega/\pi'）;ylabel（'弧度'）

3求系统

的零、极点和幅度频率响应和相位响应。

要求：

绘出零、极点分布图，幅度频率响应和相位响应曲线。

源程序：

num=[0.05280.07970.12950.12950.07970.0528];

den=[1-1.81072.4947-1.88010.9537-0.2336];

subplot（3,1,1）;

[z,p,k]=tf2zp（num,den）;

disp（'零点'）;disp（z）;

disp（'极点'）;disp（p）;

disp（'增益系数'）;disp（k）;

sos=zp2sos（z,p,k）;

disp（'二阶节'）;disp（real（sos））;

zplane（num,den）

title（'零极点分布'）

k=256;

m=[0.05280.07970.12950.12950.07970.0528];

n=[1-1.81072.4947-1.88010.9537-0.2336];

w=0:

pi/k:

pi;

h=freqz（m,n,w）;

subplot（3,2,3）;

plot（w/pi,real（h））;grid

title（'实部'）

xlabel（'\omega/\pi'）;ylabel（'幅度'）

subplot（3,2,4）;

plot（w/pi,imag（h））;grid

title（'虚部'）

xlabel（'\omega/\pi'）;ylabel（'Amplitude'）

subplot（3,2,5）;

plot（w/pi,abs（h））;grid

title（'幅度谱'）

xlabel（'\omega/\pi'）;ylabel（'幅值'）

subplot（3,2,6）;

plot（w/pi,angle（h））;grid

title（'相位谱'）

xlabel（'\omega/\pi'）;ylabel（'弧度'）

实验结果：

零点

-1.0000

-0.3018+0.9534i

-0.3018-0.9534i

0.0471+0.9989i

0.0471-0.9989i

极点

0.2788+0.8973i

0.2788-0.8973i

0.3811+0.6274i

0.3811-0.6274i

0.4910

增益系数

0.0528

二阶节

0.05280.052801.0000-0.49100

1.00000.60361.00001.0000-0.76220.5389

1.0000-0.09411.00001.0000-0.55750.8829

图形：

实验3FFT算法的应用

一、实验目的：

加深对离散信号的DFT的理解及其FFT算法的运用。

二、实验原理：

N点序列的DFT和IDFT变换定义式如下：

，

利用旋转因子

具有周期性，可以得到快速算法（FFT）。

在MATLAB中，可以用函数X=fft（x，N）和x=ifft（X，N）计算N点序列的DFT正、反变换。

三、实验内容

1.对连续的单一频率周期信号按采样频率

采样，截取长度N分别选N=20和N=16，观察其DFT结果的幅度谱。

解此时离散序列：

即N=8。

用MATLAB计算并作图，函数fft用于计算离散傅里叶变换DFT，程序如下：

q=8;

k1=[0:

1:

19];

xa1=sin（2*pi*k1/q）;

subplot（2,2,1）

plot（k1,xa1）

xlabel（'k'）;ylabel（'x1（k）'）;

xm1=fft（xa1）;xm1=abs（xm1）;

subplot（2,2,2）

stem（k1,xm1）

xlabel（'m'）;ylabel（'X1（m）'）;

k2=[0:

1:

15];

xa2=sin（2*pi*k2/q）;

subplot（2,2,3）

plot（k2,xa2）

xlabel（'k'）;ylabel（'x2（k）'）;

xm2=fft（xa2）;xm2=abs（xm2）;

subplot（2,2,4）

stem（k2,xm2）

xlabel（'m'）;ylabel（'X2（m）'）;

（d）

（c）

（b）

（a）

        图2.1

计算结果示于图2.1，（c）和（d）分别是N=16时的截取信号和DFT结果，由于截取了两个整周期，得到单一谱线的频谱；而（a）和（b）分别是N=20时的截取信号和DFT结果，由于截取了两个半周期，不再是单一谱线的频谱，而是含有丰富的谐波分量，与N=16时的差别很大；上述频谱的误差主要是由于有限长序列的离散傅里叶变换是把有限长序列作为周期序列的一个周期来表示的，隐含有周期性意义。

2.分别计算16点序列

的16点和32点DFT,绘出幅度谱图形。

源程序：

q=16;

n1=[0:

1:

15];

xn1=sin（5*pi*n1/q）;

subplot（2,2,1）

plot（n1,xn1）

xlabel（'n'）;ylabel（'x1（n）'）;

xm1=fft（xn1）;xm1=abs（xm1）;

subplot（2,2,2）

stem（n1,xm1）

xlabel（'m'）;ylabel（'X1（m）'）;

n2=[0:

1:

31];

xn2=sin（5*pi*n2/q）;

subplot（2,2,3）

plot（n2,xn2）

xlabel（'n'）;ylabel（'x2（n）'）;

xm2=fft（xn2）;xm2=abs（xm2）;

subplot（2,2,4）

stem（n2,xm2）

xlabel（'m'）;ylabel（'X2（m）'）;

实验结果：

实验4无限冲激响应数字滤波器设计

一、实验目的

掌握双线性变换法及冲激响应不变法设计IIR数字滤波器的具体设计方法及其原理；

熟悉用双线性变换法及冲激响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。

二、实验编程函数

在MATLAB中，可以用下列函数辅助设计IIR数字滤波器：

1）利用buttord和cheblord可以确定低通原型巴特沃斯和切比雪夫滤波器的阶数和截止频率；

2）[num,den]=butter（N,Wn）（巴特沃斯）；

3）[num,den]=cheby1（N,Wn）,[num,den]=cheby2（N,Wn）（切比雪夫1型和2型）可以进行滤波器的设计；

4）lp2hp，lp2bp，lp2bs可以完成低通滤波器到高通、带通、带阻滤波器的转换；

5）使用bilinear可以对模拟滤波器进行双线性变换，求得数字滤波器的传输函数系数；

6）利用impinvar可以完成冲激响应不变法的模拟滤波器到数字滤波器的转换。

注：

双线性变换法通过将数字频率

的取值范围从0到

对应到模拟频率

，也就对应于模拟域中所有可能的频率值。

双线性变换法不会出现频率混叠，但非线性关系却导致数字滤波器的频率响应不能逼真地模仿模拟滤波器的频率响应。

冲激响应不变法通过选择满足设计要求的模拟滤波器单位冲激响应h（t）的采样值的h（n），得到的被采样的冲激响应将给出与原模拟滤波器非常相近的滤波器形状。

由于该方法不可避免的要发生频率混叠现象，所以只适合设计低通和带通滤波器。

三、实验内容

1.设采样周期T=250μs，用冲激响应不变法和双线性变换法设计一个三阶巴特沃斯滤波器，其3dB边界频率为fc=1kHz。

[B,A]=butter（3,2*pi*1000,'s'）;

[num1,den1]=impinvar（B,A,4000）;

[h1,w]=freqz（num1,den1）;

[B,A]=butter（3,2/0.00025,'s'）;

[num2,den2]=bilinear（B,A,4000）;

[h2,w]=freqz（num2,den2）;

f=w/pi*2000;

plot（f,abs（h1）,'-.',f,abs（h2）,'-'）;

grid;

xlabel（'频率/Hz'）

ylabel（'幅值/dB'）

程序中第一个butter的边界频率2π×1000，为冲激响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给出了这两种设计方法所得到的频响，虚线为冲激响应不变法的结果；实线为双线性变换法的结果。

冲激响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。

图1

2.设计一巴特沃斯带通滤波器，其通带边界频率分别为f2=110kHz和f1=90kHz处的最大衰减小于3dB，在阻带边界频率分别为f3=80kHz和f4=120kHz处的最小衰减大于10dB，采样频率fs=400Hz。

w1=2*400*tan（2*pi*90/（2*400））;

w2=2*400*tan（2*pi*110/（2*400））;

wr1=2*400*tan（2*pi*80/（2*400））;

wr2=2*400*tan（2*pi*120/（2*400））;

[N,wn]=buttord（[w1w2],[wr1wr2],3,10,'s'）;

[B,A]=butter（N,wn,'s'）;

[num,den]=bilinear（B,A,400）;

[h,w]=freqz（num,den）;f=w/pi*200;

plot（f,20*log10（abs（h）））;axis（[40,160,-30,10]）;grid;

xlabel（'频率/kHz'）;ylabel（'幅度/dB'）

图2

3.利用MATLAB编程设计一个数字带通滤波器，指标要求如下：

通带边缘频率：

，

，通带衰减

。

阻带边缘频率：

，

，最小阻带衰减

。

给出IIR数字滤波器参数，绘出幅度和相位频响曲线，讨论实现形式和特点。

源程序：

wp=[0.450.65];

ws=[0.30.75];

Ap=1;

As=40;

[n1,wn1]=buttord（wp,ws,Ap,As）;

[num,den]=butter（n1,wn1）;

w=0:

pi/255:

pi;

h=freqz（num,den,w）;

g=20*log10（abs（h））;

pha=angle（h）;

subplot（1,2,1）

plot（w/pi,g）;

grid

xlabel（'频率/kHz'）;

ylabel（'幅度'）;

title（'数字带通滤波器幅频曲线'）

subplot（1,2,2）;

plot（w/pi,pha）;

grid

xlabel（'频率/kHz'）;

ylabel（'相位'）;

title（'数字带通滤波器相位曲线'）;

实验结果：

实现形式和特点分析：

实验5有限冲激响应数字滤波器设计

一、实验目的：

（1）掌握用窗函数法设计FIR数字滤波器的原理和方法。

（2）了解各种窗函数对滤波特性的影响。

二、实验原理：

如果所希望的滤波器的理想频率响应函数为

则其对应的单位抽样响应为

窗函数设计法的基本原理是用有限长单位抽样响应序列h（n）逼近hd（n），由于hd（n）往往是无限长序列，且是非因果的，所以用窗函数w（n）将hd（n）截断，并进行加权处理，得到:

h（n）就作为实际设计的FIR数字滤波器的单位抽样响应序列，其频率响应函数

为：

式中，N为所选窗函数的长度

在MATLAB中，可以用b=fir1（M,Wc,’ftype’,taper）等函数辅助设计FIR数字滤波器。

M代表滤波器阶数（单位抽样响应h（n）的长度N=M+1）；Wc代表滤波器的截止频率（对

归一化频率），当设计带通和带阻滤波器时，Wc为双元素相量；ftype代表滤波器类型，如’high’高通，’stop’带阻等；taper为窗函数类型，默认为海明窗；窗函数用blackman,hamming,hanning，kaiser产生。

三、实验内容

1.设计一FIR低通滤波器，通带截止频率

，通带衰减不大于0.25dB，阻带截止频率

，阻带衰减不小于50dB。

选择一个合适的窗函数，绘出所设计的滤波器的频率响应图。

wp=0.2*pi;ws=0.3*pi;

tr_width=ws-wp;

M=ceil（6.6*pi/tr_width）;

n=[0:

1:

M-1];

wc=（ws+wp）/2;

hd=ideal_lp（wc,M）;

w_ham=（hamming（M））';

h=hd.*w_ham;

[H,w]=freqz（h,1,1000,'whole'）;

H=（H（1:

1:

501））';

w=（w（1:

1:

501））';

mag=abs（H）;

db=20*log10（（mag+eps）/max（mag））;

delta_w=2*pi/1000;

Rp=-（min（db（1:

1:

wp/delta_w+1）））;

disp（['实际通带波动为',num2str（Rp）]）;

As=-round（max（db（ws/delta_w+1:

1:

501）））;

disp（['最小阻带衰减为',num2str（As）]）;

subplot（2,2,1）;stem（n,hd）;title（'理想抽样响应'）;

axis（[0M-1-0.10.3]）;xlabel（'n'）;ylabel（'hd（n）'）;

subplot（2,2,2）;stem（n,w_ham）;title（'海明窗'）;

axis（[0M-101.1]）;xlabel（'n'）;ylabel（'w（n）'）;

subplot（2,2,3）;stem（n,h）;title（'实际脉冲响应'）;

axis（[0M-1-0.10.3]）;xlabel（'n'）;ylabel（'h（n）'）;

subplot（2,2,4）;plot（w/pi,db）;title（'衰减幅度'）;

axis（[01-10010]）;xlabel（'以pi为单位的频率'）;ylabel（'分贝数'）;grid;

%子程序ideal_lp

functionhd=ideal_lp（wc,M）

%理想低通滤波器；hd为0到M-1之间的理想脉冲响应；wc为截止频率；M为理想滤波器的长度

alpha=（M-1）/2;

n=[0:

1:

M-1];

m=n-alpha+eps;%加一个小数以避免零做除数

hd=sin（wc*m）./（pi*m）;

FIR滤波器的阶次由过渡带宽度决定，窗函数的选择取决于阻带衰减，海明窗和布莱克曼窗均可提供大于50dB的衰减，由于海明窗提供了较窄的过渡带，因此具有较小的阶数。

在设计中没有使用通带波动值Rp=0.25dB，但必须检查设计的实际波动，验证它是否在给定容限内。

先定义Function函数，在运行主程序调用函数

实验结果：

实际通带波动为0.043502

最小阻带衰减为50

2.（选做） 设计一个数字带通滤波器，指标要求如下：

通带边缘频率：

，

，通带衰减

。

阻带边缘频率：

，

，最小阻带衰减

。

要求：

给出FIR数字滤波器的冲激响应，绘出幅度和相位频响曲线，讨论实现形式和特点。

源程序：

f=[0.30.450.650.8];

a=[010];

dev=[0.010.10870.01];

fs=2;

[n,wn,bta,ftype]=kaiserord（f,a,dev,fs）;

h1=fir1（n,wn,ftype,kaiser（n+1,bta）,'noscale'）;

[hh1,w1]=freqz（h1,1,256）;

figure

（1）

subplot（2,1,1）

plot（w1/pi,20*log10（abs（hh1）））

grid

xlabel（'归一化频率'）;ylabel（'幅度'）;

subplot（2,1,2）

plot（w1/pi,angle（hh1））

grid

xlabel（'归一化频率'）,ylabel（'相位'）;

实现形式和特点的分析：