数字信号处理实验报告3

实验3 离散时间系统的特性离散时间系统的输出

一实验内容

给定离散时间系统

y(n)=x(n)+x(n-1)+0.9y(n-1)-0.81y(n-2)

1）求系统函数，调用zplane，画出系统的零极图；

2）调用freqz，画出系统的幅度响应和相位响应。

3）系统输入信号为

x(n)=sin(n*pi/3)+5cos(n*pi) (0≤n ≤ 49），

用四种方法计算系统的输出y(n)

a. 用差分方程递推

b. 调用conv函数

c. 调用filter函数

d. 利用DFT计算

观察结果，判断四种方法得到的结果是否一致？如果使得上述方法结果一致，编程时有关参数应该怎样选取？

二编程原理、思路和公式

1．对离散时间系统y(n)=x(n)+x(n-1)+0.9y(n-1)-0.81y(n-2)，用B表示x(n)，即系统函数的分母部分，用A表示y(n)，即系统函数的分子部分。调用zplane(B,A

函数，即可绘制由行向量B和A构成的系统函数的零极点分布图。

H时，调用[h, w] = freqz(b, a, n)可2．在求解离散时间系统的频率响应函数)

( j e

得到数字滤波器的n点幅频响应值，这n个点均匀地分布在[0, p]上，并将这n个

频点的频率记录在w中，相应的频响值记录在h中。即可得到系统函数的幅度响应和相位响应。

3. 差分方程递推法：首先定义数组x和y x=zeros(1,50);y=zeros(1,50);

然后用for循环定义输入信号x(n)=sin(n*pi/3)+5cos(n*pi) (0≤n ≤49），即i=0:49 x(i+1)=sin(i*pi/3)+5*cos(i*pi);由递推法写出不能循环的前三项。再应用for循环求y i=4:49；y(i)=x(i)+x(i-1)+0.9*y(i-1)-0.81*y(i-2);

4. 调用filter函数：先分别写出x y的定义x=sin(n*pi/3)+5*cos(n*pi);

B=[1,1];A=[1,-0.9,0.81];调用filter函数实现对信号的滤波处理

y = filter(B, A x); %对于由矢量B A定义的数字系统，当输入信号为x时，对x中的数据进行滤波，结果放于y中。

5. 调用conv函数：先分别写出x y的定义x=sin(n*pi/3)+5*cos(n*pi);

B=[1,1];A=[1,-0.9,0.81] 然后求解系统的冲击响应hn [hn,n]=impz(B,A,14); 再调用conv函数用线性卷积求系统响应y=conv(hn,x)

6.利用DFT计算：先编写两个子函数DFT 和IDFT 然后定义函数x和y，求系统的冲击响应hn [hn,n]=impz(B,A,14);然后直接调用两个子函数求输出y

三、程序脚本，并注释

1 B=[1,1];A=[1,-0.9,0.81];%给定义函数

subplot(1,3,1);

zplane(B,A); %绘制由行向量AB构成的系统函数确定的零极点分布图xlabel('实部Re');ylabel(‘虚部Im');title('系统函数零极点分布');grid;

[H,w]=freqz(B,A,'whole'); %求解离散事件系统的频率响应

subplot(1,3,2);plot(w/pi,abs(H),'linewidth',2);grid;

xlabel('w/pi');ylabel('|X(e^j^w)|');title('幅频响应特性');

subplot(1,3,3);plot(w/pi,angle(H),'linewidth',2);grid;

xlabel('w/pi');ylabel('Angle[X(e^j^w)]');title('相频响应特性');

2 .N=50;n=0:N-1;x=sin(n*pi/3)+5*cos(n*pi);

B=[1,1];A=[1,-0.9,0.81];

y=filter(B,A,x); %对于由AB组成的数字系统，当输入信号为x时，对x中的数据进行滤波，结果放在y中

stem(n,y);

3 . N=50;n=0:N-1;x=sin(n*pi/3)+5*cos(n*pi);

B=[1,1];A=[1,-0.9,0.81];

[hn,n]=impz(B,A,14); %求解数字响应的冲击h,取样点数由n确定

y=conv(hn,x); %用线性卷积求系统响应

stem(y);

4.dft子函数

function Xk=DFT(xn,N)

n=[0:1:N-1];

k=[0:1:N-1];

WN=exp(-j*2*pi/N);

nk=n'*k;

WNnk=WN.^nk;

Xk=xn*WNnk;

5.idft子函数

function xn=IDFT(Xk,N)

n=[0:1:N-1];

k=[0:1:N-1];

WN=exp(-j*2*pi/N);

nk=-n'*k;

WNnk=WN.^nk;

xn=Xk*WNnk/N;

6. A=[1,-0.9,0.81];

B=[1,1];

N=100;

n=0:N-1;

xn1=sin(n*pi/3)+5*cos(n*pi);

hn=impz(B,A,50); %求解数字响应的冲击h,取样点数50

hn1=[hn' zeros(1,50)];

Hk=DFT(hn1,N);

Xk=DFT(xn,N);

Yk=Xk.*Hk;

y=IDFT(Yk,N);

figure;

stem(n,y);title('DFT ');

7. n=0:49;

x=zeros(1,50);

y=zeros(1,50); %定义数组x y

for i=0:49

x(i+1)=sin(i*pi/3)+5*cos(i*pi);

end

y(1)=x(1);

y(2)=x(2)+x(1)+0.9*y(1);

y(3)=x(3)+x(2)+0.9*y(2)-0.81*y(1); %先给出前三个不能循环的值for i=4:49

y(i)=x(i)+x(i-1)+0.9*y(i-1)-0.81*y(i-2); %递推法求出函数y end

n=0:49;

stem(n,y); %x y长度相同

四仿真结果、图形