数字信号处理实验报告

实验一：离散时间信号的表示与运算

一：实验内容、原理描述及实验结果 1. 离散时间信号的表示

离散时间信号定义为一时间函数，它只在某些离散的瞬时给出函数值，而在其他处无定义。因此，它是时间上不连续按一定先后次序排列的一组数的集合，故称为时间序列，简称序列，通常表示为

{x （n ）} -∞

例1.1 用matlab 编写生成单位脉冲序列函数的程序，n ∈(-5,5)。

程序代码如下：n0=0;n1=-5;n2=5;n=[n1:n2];nc=length(n);x=zeros(1,nc);

for i=1:nc

if n(i)==n0 x(i)=1; end

end

stem(n,x)

xlabel('n');ylabel('x(n)');title('delta sequence'); grid

或者：n0=0;n1=-5;n2=5;n=[n1:n2];x=[(n-n0)==0];stem(n,x)；

xlabel('n');ylabel('x(n)');title('delta sequence');grid 图形如下：

n

x (n )

（2）单位阶跃序列

例1.2 用matlab 编写生成单位阶跃序列函数的程序，n ∈(-5,5)。

程序代码如下：n0=0;n1=-5;n2=5;n=[n1:n2];x=[(n-n0)>=0];stem(n,x)

xlabel('n');ylabel('x(n)');title('step sequence');grid

图形如下：

n

x (n )

（3）单位斜坡序列

例1.3 用matlab 编写生成单位斜坡序列函数的程序，n ∈(0,5)。 程序代码如下：n1=0;n2=5;n=[n1:1:n2];x=n;stem(n,x)

xlabel('n');ylabel('x(n)');title('ramp sequence');grid

图形如下：

n

x (n )

（4）正余弦序列

例1.4 用matlab 编写正弦序列x(n)=5sin(0.1πn+π/3)函数的程序。

程序代码如下：n=[0:50];x=5*sin(0.1*pi*n+pi/3);stem(n,x)

xlabel('n');ylabel('x(n)');title('sine sequence');grid

图形如下：

n

x (n )

（5）实指数序列

例1.5 用matlab 编写生成实指数序列的程序。 程序代码如下：n=[0:20];x=(0.78).^n;stem(n,x)

xlabel('n');ylabel('x(n)');title('real power sequence');grid

图形如下：

n

x (n )

（6）复指数序列

例1.6 用matlab 编写生成复指数序列x(n)=2e (-0.08+j0.26)n 的程序, n ∈(-20,20)。 程序代码如下：clf

n=[-20:20];a=-0.08+0.26*j;x=2*exp(a*n);b=real(x);c=imag(x);d=abs(x);e=(180/pi)*angle(x); subplot(221)

stem(n,b);xlabel('n');title('real part');grid subplot(222)

stem(n,c);xlabel('n');title('imaginary part');grid

stem(n,d);xlabel('n');title('magnitude');grid subplot(224)

stem(n,e);xlabel('n');title('phase');grid 图形如下：

-20

-10

01020

n

real part

-20

-10

01020

n imaginary part

n

magnitude

n

phase

（7）随机序列

MATLAB 提供两个产生随机序列的函数：

rand(1,n)用于产生[0,1]上均匀分布随机序列，长度为n 。

randn(1,n)用于产生均值为零、方差为一的高斯随机序列，即白噪声序列，长度为n 。 例1.7 用matlab 编写生成长度为14的随机序列。 程序代码如下：n=[1:14];x=rand(1,14);

subplot(211)

stem(n,x);xlabel('n');ylabel('x(n)');title('rand sequence');grid x=randn(1,14); subplot(212)

stem(n,x);xlabel('n');ylabel('x(n)');title('randn sequence');grid

图形如下：

n

x (n )

rand sequence

n

x (n )

若x(n)=x(n+T),则x(n)为周期序列，周期为T。

用matlab产生一个信号x(n)的P个周期的周期信号xp(n)，用matlab语言可直接产生

。

2.离散时间信号的运算

离散时间信号一般可用数的序列表示，和连续时间信号类似，序列也可进行运算。（1）信号加

序列的加性条件是x1（n）和x2（n）具有相等的长度，且在相同的采样位置上相加，否则需进行转换。

例2.1 用matlab实现信号加的函数sigadd（）的程序。

程序代码如下：function[y,n]=sigadd(x1,n1,x2,n2)

n=min(min(n1),min(n2)):max(max(n1),max(n2));y1=zero(1,length(n));y2=y1;

y1(find((n>=min(n1))& (n<=max(n1))==1))=x1;

y2(find((n>=min(n2))& (n<=max(n2))==1))=x2;

y=y1+y2;

（2）信号乘

序列的乘性条件是x1（n）和x2（n）具有相等的长度，且在相同的采样位置上相乘，否则需进行转换。

例2.2 用matlab实现信号乘的函数sigmult（）的程序。

程序代码如下：function[y,n]=sigmult(x1,n1,x2,n2)

n=min(min(n1),min(n2)):max(max(n1),max(n2));y1=zeros(1,length(n));y2=y1;

y1(find((n>=min(n1))& (n<=max(n1))==1))=x1;

y2(find((n>=min(n2))& (n<=max(n2))==1))=x2;

y=y1.*y2;

（3）信号移位

y(n)=x(n+n0)

序列y(n)相对于序列x(n)左移n0个采样周期。

例2.3 用matlab实现序列移位的函数sigshift（）的程序。

程序代码如下：function[y,n]=sigshift(x,m,n0)

n=m+n0;y=x;

（4）序列折叠

y(n)=x(-n)

例2.4 用matlab实现序列折叠函数sigfold（）的程序。

程序代码如下：function[y,n]=sigfold(x,m,n0)

Y=fliplr(x);n=-fliplr(n);

（5） 序列的奇偶性

任何一个序列x(n)都可以分解为偶分量x e (n)和奇分量x o (n)两部分之和。 x(n)= x e (n)+ x o (n)

例2.5 用matlab 编写函数sigevenodd （）用于将序列分成奇序列和偶序列两部分。 程序代码如下：function[xeven,xodd,m]=sigevenodd(x,n) If(imag(x)~=0)

error(‘x is not a real sequence ’); end

m=fliplr(n);m1=min([m,n]);m2=max([m,n]);m=m1:m2;nm=n(1)-m(1); n1=1:length(n);x1(n1+nm)=x;x=x1;xeven=0.5*(x+fliplr(x)); xodd=0.5*(x-fliplr(x)); 例2.6 设单位阶跃序列

将其分解为偶分量和奇分量。

程序代码如下：clf

n0=0;n1=-10;n2=10;n=[n1:n2];x=[(n-n0)>=0]; subplot(311)

stem(n,x);xlabel('n');ylabel('x(n)');title('step sequence');grid [xeven,xodd,m]=sigevenodd(x,n); subplot(312)

stem(m,xeven);xlabel('m');ylabel('x even(n)');title('even part'); grid

subplot(313)

stem(m,xodd);xlabel('m');ylabel('x odd(n)');title('odd part');grid

图形如下：

n x (n )

-10

-8

-6

-4

-2

02

4

6

8

10

m x e v e n (n )

m

x o d d (n )

实验二：离散时间系统分析

例1.1使用MATLAB 的residuez 函数，求出

1218153332

5644162)(234234-+-+++++=

z z z z z z z z z X 的部分分式展开和。

代码：B=[2 16 44 56 32];A=[3 3 -15 18 -12];[R,P,K]=residuez(B,A) 结果：R = -0.0177 ，9.4914 ，-3.0702 + 2.3398i ，-3.0702 - 2.3398i

P =-3.2361，1.2361 ，0.5000 + 0.8660i ，0.5000 - 0.8660i

K =-2.6667

例1.2使用MATLAB 画出下列因果系统的系统函数的零极点图，并判断系统的稳定性。

48.096.15.29

.06.12)(23

2-+---=z z z z z z H

代码：B=[2,-1.6,-0.9];A=[1,-2.5,1.96,-0.48];zplane(B,A),grid on legend('零点','极点')；title('零极点分布图') 图形：

Real Part

I m a g i n a r y P a r t

零极点分布图

由图可以看出有极点在单位圆外，所以系统不稳定。 例1.3已知IIR 滤波传递函数为试确定其零极点表

达式。

代码： b=[2 3 4];a=[1 3 3 1];f=filt(b,a)；q=roots(b)

p=roots(a) k=b(1)/a(1) bb=k*poly(q) aa=poly(p)

结果： Transfer function:

2 +

3 z^-1 +

4 z^-2 --------------------------

1 + 3 z^-1 + 3 z^-

2 + z^-3

Sampling time: unspecified

q =-0.7500 + 1.1990i

-0.7500 - 1.1990i

p =-1.0000

-1.0000 + 0.0000i

-1.0000 - 0.0000i

k = 2

bb =2.0000 3.0000 4.0000

aa =1.0000 3.0000 3.0000 1.0000

例1.4已知滤波器的传递函数试求该滤波器的零极点和状态空间表达式

代码：num=[2 3];den=[1 0.4 1];f=filt(num,den)

disp('Convert to zero-pole-gain')

[z,p,k]=tf2zp(num,den)

disp('Convert to state space')

[A,B,C,D]=tf2ss(num,den)

结果：Sampling time: unspecified

Convert to zero-pole-gain

z =-1.5000

p =-0.2000 + 0.9798i

-0.2000 - 0.9798i

k =2

Convert to state space

A =-0.4000 -1.0000

1.0000 0

B =1

C =2 3

D =0

例1.5已知系统传递函数求其部分分式展开式。

代码：b=[-4 8];a=[1 6 8];[r,p,k]=residuez(b,a)

结果：r =-12；8

p =-4

-2

k = [ ]

即

例1.6已知滤波器的二阶形式确定其空间状态表达式

和零极点表达式。

代码：sos=[1 1 1 1 0 -1;-2 3 1 1 10 1];[num,den]=sos2tf(sos);f-filt(num,den)

disp('convert to zero-pole form')

[z,p,k]=sos2zp(sos)

结果：Transfer function:

4 + 22.8 z^-1 + 29.6 z^-2 - 25.8 z^-3 - 36.6 z^-4 - 7.4 z^-

5 - z^-6

-------------------------------------------------------------------

1 + 10.4 z^-1 + 5 z^-

2 - 5 z^-4 - 10.4 z^-5 - z^-6

Sampling time: unspecified

convert to zero-pole form

z =-0.5000 + 0.8660i

-0.5000 - 0.8660i

1.7808

-0.2808

p =1.0000

-1.0000

-9.8990

-0.1010

k =-2

例1.7把下面函数传递形式的IIR滤波器转变为格型、梯型结构

代码：num=[1 2 2 1];den=[1 13/24 5/8 1/3];%Transfer function to lattice filter conversion [k,v]=tf2latc(num,den)

结果：k =0.2500

0.5000

0.3333

v =-0.2695

0.8281

1.4583

1.0000

例1.8已知h=[1 2 3 2 ],x为4*1随机序列，求其卷积矩阵。

程序：h=[1 2 3 2];n=4;randn('seed',0);x=randn(n,1);C=convmtx(h',n)

y1=C*x

y2=conv(h,x)

结果：C =1 0 0 0

2 1 0 0

3 2 1 0

2 3 2 1 0 2 3 2 0 0 2 3 0 0 0 2 y1 =1.1650 2.9567 4.8236 4.7122 2.1821 1.2050 0.7032 y2 =1.1650 2.9567 4.8236 4.7122 2.1821 1.2050 0.7032

例1.9已知滤波器的传递函数1

8.0115

.0)(--=

z z H ，输入信号为)()05.0sin(2)(t w t t x +=π，

W(t)为随机信号，幅值为0.2.试绘出滤波器输出信号。

程序：b=0.15;a=[1 -0.8];n=[0:100];x=2*sin(0.05*pi*n)+0.2*randn(1,101);imp=[1;zeros(100,1)];

h=filter(b,a,imp);yc=conv(h,x);y=yc(1:101);y1=filter(b,a,x);plot(n,x,'r',n,y1,'m'); xlabel('n');ylabel('x y yc');title('time response');grid 图形：

10

20

30

40

5060

70

80

90

100

n

x y y c

time response

实验三：Z 变换

例1.1求z 的逆变换：已知2

12

143)(----+-+=z

z z z z X ，1z >，求其逆变换。 代码：num=[0 1 1];den=[3 -4 1];[r,p,k]=residuez(num,den) 结果：r =1.0000 -2.0000 p =1.0000 0.3333 k =1

则实验结果为：

1

1

311211

1)(----

-+

=z z z x

例1.2利用z 变换求卷积：已知序列x(n)={2,3,4},y(n)={3,4,5,6},求h(n)=x(n)*y(n)。

程序：xn=[2 3 4];Xz=[2 3 4]; yn=[3 4 5 6];Yz=[3 4 5 6]; xnCyn=conv(xn,yn) XzMYz=conv(Xz,Yz) 结果：

xnCyn =6 17 34 43 38 24 XzMYz =6 17 34 43 38 24

则其结果为：h(n)=x(n)*y(n)={ 6 17 34 43 38 24}

例1.3利用变换解差分方程：求解系统差分方程y(n)=x(n)-5x(n-1)+8x(n-3)。

解：其取Z 变换为：H(z)=3

1851--+-z z

程序：b=[1 -5 0 8];N=30;n=0:N-1;x=0.8.^n;y=filter(b,1,x);stem(n,y);grid 图形：

例1.4求解差分方程：y(n)-0.4y(n-1)-0.45y(n-2)=0.45x(n)+0.4x(n-1)-x(n-2),其中

)(8.0)(n u n x n =，初始状态y(-1)=0,y(-2)=1,x(-1)=1,x(-2)=2.

解：将方程两边进行变换得：

2

12

145.04.014.045.0)(-------+=z z z z z H

代码：num=[0.45 0.4 -1];den=[1 -0.4 -0.45];x0=[1 2];y0=[0 1];N=50;

n=[0:N-1];x=0.8.^n ；Zi=filtic(num,den,y0,x0);[y,Zf]=filter(num,den,x,Zi); plot(n,x,'r-',n,y,'b--');title('response')

xlabel('n');ylabel('x(n)-y(n)');legend('Input x','Output y',1);grid

图形：

response

n

x (n )-y (n )

实验四 ：DTF 变换的性质及应用

例1.1周期序列——离散傅里叶级数：求下面周期序列的DFS 表达式

（1）求其

主值序列；（2）求

的主值序列。

（1）程序如下：function[Xk]=dfs(xn,N)

n=[0:1:N-1];k=n;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk; Xk=xn*WNnk;xn=[0 1 2 3];N=4;Xk=dfs(xn,N)'

结果：Xk =6.0000

-2.0000 - 2.0000i -2.0000 + 0.0000i -2.0000 + 2.0000i

（2）程序如下：function[xn]=idfs(Xk,N)

n=[0:1:N-1];k=n;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^(-nk); xn=Xk*WNnk/N;Xk=[6.0000-2.0000+2.0000i-2.0000-2.0000-2.0000i];

k

|X (k )|

k

|X (k )|

N=4;xn=idfs(Xk,N)'

结果：Xk =0

1.0000 - 0.0000i

2.0000 + 0.0000i

3.0000 + 0.0000i

例1.2.周期方波序列：

其中，N 是基本周期，L/N 是占空比，绘出下列情况下的()X k 曲线： (a) L=5,N=20; (b) L=5,N=40; (c) L=5,N=60; (d) L=7,N=20. 代码：L=5;N=20;k=[-N/2:N/2]；xn=[ones(1,L),zeros(1,N-L)]；Xk=dfs(xn,N)；

magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)])； subplot(2,2,1);

stem(k,magXk);axis([-N/2,N/2,0,0.3]);xlabel('k');ylabel('|X(k)|');

L=5;N=40;k=[-N/2:N/2]；xn=[ones(1,L),zeros(1,N-L)]；Xk=dfs(xn,N)； magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]) subplot(2,2,2);

stem(k,magXk);axis([-N/2,N/2,0,0.2]);xlabel('k');ylabel('|X(k)|');

L=5;N=60;k=[-N/2:N/2]；xn=[ones(1,L),zeros(1,N-L)]；Xk=dfs(xn,N)； magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]) subplot(2,2,3);

stem(k,magXk);axis([-N/2,N/2,0,0.1]);xlabel('k');ylabel('|X(k)|');

L=7;N=20;k=[-N/2:N/2]；xn=[ones(1,L),zeros(1,N-L)]；Xk=dfs(xn,N) magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]) subplot(2,2,4);

stem(k,magXk);axis([-N/2,N/2,0,0.5]);xlabel('k');ylabel('|X(k)|'); 图形：

例1.3有限长序列——离散傅里叶变换：

已知序列()cos(0.48)cos(0.52),0100x n n n n ππ=+≤≤,试绘制x(n)及它的傅立叶变换

()X k 图。

解：函数调用：dft.m 文件

function [Xk]=dft(xn,N)

n=[0:1:N-1];k=n;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnk; 代码：N=100;n=0:N-1;xn=cos(0.48*pi*n)+cos(0.52*pi*n);Xk=dft(xn,N);

magXk=abs(Xk);phaXk=angle(Xk); subplot(121);

plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n) N=100'); subplot(122);

k=0:length(magXk)-1;plot(k,magXk);xlabel('k');ylabel('|X(k)|');title('X(k) N=100'); 图形：

-2

-1.5-1

-0.500.51

1.52n

x (n )

x(n) N=100

k

|X (k )|

X (k) N=100

实验五 ：离散时间系统的频域响应

例1.1z域重构—由离散傅里叶变换到z变换：

已知序列的x(n)=(0.9)n u(n),(1)球序列的z变换；（2）当频率采样点数M=5，10，20，25时，序列的傅立叶变换X（k）及其逆变换x(n)；(3)比较M对x（n）的影响。

解：调用函数idft：

function[xn]=idft(Xk,N)

n=[0:1:N-1];k=n;WN=exp(-j*2*pi/N);nk=n'*k;

WNnk=WN.^(-nk);xn=Xk*WNnk/N;

代码：clf

a=0.8；Nsum=40;n=0:Nsum-1;

for i=1:4

switch i

case 1

N=5;

case 2

N=10;

case 3

N=20;

case 4

N=40;

end

k=0:N-1;zk=exp(j*2*pi*k/N);Xk=zk./(zk-a);xnN=real(idft(Xk,N));xn=xnN'*ones(1,Nsum/N);

xn=(xn(:))';gsptext=['22' int2str(i)];

subplot(gsptext)

stem(n,xn);axis([0,40,-0.1,1.5]);xtext=['IDFT N=' int2str(N)];xlabel('n');ylabel('x(n)');title(xtext); hold on

plot(n,zeros(1,length(n)));

hold off

end

图形：

n

x (n )

IDFT N=5

n

x (n )

IDFT N=10

n

x (n )

IDFT N=20

n

x (n )

IDFT N=40

由图可以看出，当频域采样点数N 小于10时，引起较大的时域混叠，恢复的x(n)会产生较大失真。当N 较大，不会引起明显的混叠。

例 2.2 离散时间系统的频域响应：已知系统差分方程为

y(n)=x(n)+2x(n-1)-0.5y(n-1)-0.25y(n-2),绘制系统的幅值响应和相位响应曲线。 解：对系统差分方程两边进行z 变换可得系统传输函数为

代码：方法一：

b=[1 2 1];a=[1 0.5 0.25];m=0:length(b)-1;n=0:length(a)-1;K=512;k=1:K;w=pi*k/K;

num=b*exp(-j*m'*w);den=a*exp(-j*n'*w);H=num./den;magH=20*log(abs(H));phaH=angle(H); subplot(121)

plot(w/pi,magH);xlabel('w(pi)');ylabel('H(jw)')title('Magnitude');grid subplot(122)

plot(w/pi,phaH*180/pi);xlabel('w(pi)');ylabel('H(jw)')title('phase')grid

图形：

-200

-150

-100

-50

50

w(pi)

H (j w )

Magnitude

w(pi)

H (j w )

phase

方法二：b=[1,2,1];a=[1,0.5,0.25];K=512;figure(1);freqz(b,a,K);

[H,w]=freqz(b,a,K);magH=20*log(abs(H));phaH=angle(H);figure(2); subplot(1,2,1) plot(w/pi,magH);

xlabel('w(pi)');ylabel('|H(jw)| dB');title('Magnitude');grid subplot(1,2,2)

plot(w/pi,phaH*180/pi);

xlabel('w(pi)');ylabel('

Normalized Frequency (?π rad/sample)

P h a s e (d e g r e e s

)

Normalized Frequency (?π rad/sample)

M a g n i t u d e (d B )

w(pi)

|H (j w )| d B

Magnitude

w(pi)

Phase

由图形可知，使用freqz 函数所求频域响应与方法一求的结果完全一致。

实验六 ：FFT 与Chirp-z 变换

例1.1已知序列x=[4,3,2,6,7,8,9,0],求X(k)=DFT[x （n ）]。 代码：N=8;n=0:N-1;xn=[4,3,2,6,7,8,9,0]';Xk=fft(xn) 结果：Xk =39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000

4.7782 + 7.7071i

0 + 5.0000i

-10.7782 - 6.2929i

例1.2已知信号由15Hz幅值0.5的正弦信号和40Hz幅值2的正弦信号组成，数据采样频率为100Hz，试分别绘制N=128点DFT的幅频图和N=1024点的DFT幅频图。

由题意的，信号可写为x=0.5sin(2πf1t)+2 sin(2πf2t),其中，f1=15Hz,f2=40Hz。

代码：clf

fs=100;N=128;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);

y=fft(x,N);mag=abs(y);f=(0:length(y)-1)'*fs/length(y);

subplot(221)

plot(f,mag);xlabel('frequence');ylabel('magnitude');title('N=128')；

grid

subplot(222)

plot(f(1:N/2),mag(1:N/2));xlabel('frequence');ylabel('magnitude');title('N=128')；

grid

fs=100;N=1024;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);mag=abs(y); f=(0:length(y)-1)'*fs/length(y);

subplot(223)

plot(f,mag);xlabel('frequence');ylabel('magnitude');title('N=1024')；

grid

subplot(224)

plot(f(1:N/2),mag(1:N/2));xlabel('frequence');ylabel('magnitude');title('N=1024')；

grid

图形：

50100

050

100

150

frequence m a g n i t u d e

N=128

204060

frequence m a g n i t u d e

N=128

0200400600800frequence

m a g n i t u d e

N=1024

frequence

m a g n i t u d e

N=1024

结果分析：用FFT 算法实现离散傅里叶变换，求结果与所选点数N 有关。当N 大于序列点数时，先补零到N 点，最够结果幅值的大小也随N 点增加而增大。当N 小于序列点数时，则要对函数截断。从图中可知，不论N 为何值，其幅频特性，总是对抽样频率的一半呈偶对称，此频率计为奈圭斯特频率。

例1.3对信号()sin(2*40)sin(2*15)x t t t ππ=+进行DFT ，对其结果进行IDFT ，并将IDFT 的结果和原信号进行比较。 代码：clf

fs=100;N=128;n=0:N-1;t=n/fs;x=sin(2*pi*15*t)+sin(2*pi*40*t); subplot(221)

plot(t,x);xlabel('t(s)');ylabel('x');title('Original signal');grid y=fft(x,N);mag=abs(y);f=(0:length(y)-1)'*fs/length(y); subplot(222)

plot(f(1:N/2),mag(1:N/2));xlabel('Frequency(Hz)');ylabel('Magnitude'); title('FFT to Original signal');grid

xifft=ifft(y);magx=real(xifft);ti=[0:length(xifft)-1]/fs; subplot(223)

plot(ti,magx);xlabel('t(s)');ylabel('x');title('signal from IFFT');grid; yif=fft(xifft,N);mag=abs(yif);f=(0:length(y)-1)'*fs/length(y); subplot(224)

plot(f(1:N/2),mag(1:N/2));xlabel('Frequency(Hz)');ylabel('Magnitude'); title('FFT to signal from IFFT');grid