MATLAB实验二(修改)

实验二 信号的表示及其基本运算

一、实验目的

1、掌握连续信号及其MATLAB 实现方法；

2、掌握离散信号及其MA TLAB 实现方法

3、掌握离散信号的基本运算方法，以及MA TLAB 实现

4 熟悉应用MATLAB 实现求解系统响应的方法 4、了解离散傅里叶变换的MA TLAB 实现 5、了解IIR 数字滤波器设计 6、了解FIR 数字滤波器设计1

二、实验设备

计算机，Matlab 软件 三、实验内容

（一）、 连续信号及其MATLAB 实现 1、 单位冲激信号

()0,0()1,0

t t t dt ε

ε

δδε-?=≠??=?>???

例1.1：单位冲击信号的MATLAB 实现程序如下：

t1=-4; t2=4; t0=0; dt=0.01; t=t1:dt:t2; n=length(t); x=zeros(1,n);

x(1,(-t0-t1)/dt+1)=1/dt; stairs(t,x);

axis([t1,t2,0,1.2/dt]);

2、 任意函数

()()()f t f t d τδττ+∞

-∞

=-?

例1.2：用MA TLAB 画出如下表达式的脉冲序列

()0.4(2)0.8(1) 1.2() 1.5(1) 1.0(2)0.7(3)f n n n n n n n δδδδδδ=-+-+++++++

3 单位阶跃函数

1,0()0,

t u t t ?≥?=?

例1.3：用MA TLAB 实现单位阶跃函数

clear all;

t=-0.5:0.001:1; t0=0;

u=stepfun(t,t0); plot(t,u)

axis([-0.5 1 -0.2 1.2])

4 斜坡函数

0()()g t B t t =-

例1.4：用MA TLAB 实现g(t)=3(t-1) clear all; t=0:0.01:3; B=3; t0=1;

u=stepfun(t,t0); n=length(t); for i=1:n

u(i)=B*u(i)*(t(i)-t0); end plot(t,u)

axis([-0.2 3.1 -0.2 6.2])

5 抽样信号 抽样信号Sa(t)=sin(t)/t 在MATLAB 中用 sinc 函数表示。 定义为 )/(sin )(πt c t Sa =

t=-3*pi:pi/100:3*pi; ft=sinc(t/pi); plot(t,ft); grid on;

axis([-10,10,-0.5,1.2]); %定义画图范围，横轴，纵轴 title('抽样信号') %定义图的标题名字

6 指数函数

()at f t Ae =

例1.5：用MA TLAB 实现0.5()3t

f t e =

7 正弦函数

2()cos(

)t

f t A T π?=+ 例1.6：用MA TLAB 实现正弦函数f(t)=3cos(10πt+1) 8 虚指数信号

例 虚指数信号 调用格式是f=exp((j*w)*t) t=0:0.01:15;

w=pi/4; X=exp(j*w*t);

Xr=real(X); %取实部 Xi=imag(X); %取虚部 Xa=abs(X); %取模 Xn=angle(X); %取相位

subplot(2,2,1),plot(t,Xr),axis([0,15,-(max(Xa)+0.5),max(Xa)+0.5]), title('实部');

subplot(2,2,3),plot(t,Xi),axis([0,15,-(max(Xa)+0.5),max(Xa)+0.5]), title('虚部');

subplot(2,2,2), plot(t,Xa),axis([0,15,0,max(Xa)+1]),title('模');

subplot(2,2,4),plot(t,Xn),axis([0,15,-(max(Xn)+1),max(Xn)+1]),title('相角'); %subplot(m,n,i) 命令是建立m 行n 列画图窗口，并指定画图位置i 9 复指数信号

例 复指数信号 调用格式是f=exp((a+j*b)*t) t=0:0.01:3;

a=-1;b=10;

f=exp((a+j*b)*t);

subplot(2,2,1),plot(t,real(f)),title('实部') subplot(2,2,3),plot(t,imag(f)),title('虚部') subplot(2,2,2),plot(t,abs(f)),title('模') subplot(2,2,4),plot(t,angle(f)),title('相角')

（二）、离散信号及其MATLAB 实现 1、 单位冲激序列

1,0()0,

n n n δ?=?=?

≠??

例2.1：用MA TLAB 产生64点的单位冲激序列

clear all; N=64;

x=zeros(1,N); x(1)=1; xn=0:N-1; stem(xn,x)

axis([-1 65 0 1.1])

2、 任意序列

()()()m f n f m n m δ∞

=-∞

=

-∑

例2.2：用MA TLAB 画出如下表达式的脉冲序列

()8.0() 3.4(1) 1.8() 5.6(3) 2.9(4)0.7(5)f n n n n n n n δδδδδδ=+-++-+-+-

3、 单位阶跃序列

1,0()0,

n u n n ?≥?=?

例2.3：用MA TLAB 实现单位阶跃函数

4、 斜坡序列

0()()g n B n n =-

例2.4：用MA TLAB 实现g(n)=3(n-4)点数为32的斜坡序列 clear all; N=32; k=4 B=3; t0=1;

x=[zeros(1,k) ones(1,N-k)]; for i=1:N

x(i)=B*x(i)*(i-k); end

xn=0:N-1; stem(xn,x)

axis([-1 32 0 90])

5、 正弦序列

()sin(2)x n A fn π?=+

例2.5：用MA TLAB 实现幅度A=3，频率f=100，初始相位Φ=1.2，点数为32的正弦信号 6、 实指数序列

()n x n Aa =

例2.6：用MA TLAB 实现0.7

()3x n e =，点数为32的实指数序列

clear all; N=32; A=3; a=0.7; xn=0:N-1; x=A*a.^xn; stem(xn,x)

7、 复指数序列

()(),a j n x n Ae n ω+=?

例2.7：用MA TLAB 实现幅度A=3，a=0.7，角频率ω=314，点数为32的实指数序列 clear all;

N=32;

A=3;

a=0.7;

w=314;

xn=0:N-1;

x=A*exp((a+j*w)*xn);

stem(xn,x)

8、随机序列

利用MATLAB产生两种随机信号：

rand(1,N)在区间上产生N点均匀分布的随机序列

randn(1,N)产生均值为0，方差为1的高斯随机序列，即白噪声序列

例2.8：用MA TLAB产生点数为32的均匀分布的随机序列与高斯随机序列

clear all;

N=32;

x_rand=rand(1,N);

x_randn=randn(1,N);

xn=0:N-1;

figure(1);

stem(xn,x_rand)

figure(2);

stem(xn,x_randn)

（三）、离散信号的基本运算

1、信号的延迟

给定离散信号x(n)，若信号y(n)定义为：y(n)=x(n-k)，那么y(n)是信号x(n)在时间轴上右移k个抽样周期得到的新序列。

例3.1：正弦序列y(n)=sin(100n)右移3个抽样周期后所得的序列，MATLAB程序如下：clear all;

N=32;

w=100;

k=3;

x1=zeros(1,k);

xn=0:N-1;

x2=sin(100*xn);

figure(1)

stem(xn,x2)

x=[x1 x2];

axis([-1 N -1.1 1.1])

N=N+k;

xn=0:N-1;

figure(2)

stem(xn,x)

axis([-1 N -1.1 1.1])

利用for 循环语句实现周期延迟. 2、 信号相加

若信号12()()()x n x n x n =+，值得注意的是当序列1()x n 和2()x n 的长度不相等或者位置不对应时，首先应该使两者的位置对齐，然后通过zeros 函数左右补零使其长度相等后再相加 例3.2：用MA TLAB 实现两序列相加 clear all; n1=0:3

x1=[2 0.5 0.9 1]; figure(1) stem(n1,x1)

axis([-1 8 0 2.1] ) n2=0:7

x2=[ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7]; figure(2) stem(n2,x2)

axis([-1 8 0 0.8] ) n=0:7;

x1=[x1 zeros(1,8-length(n1))]; x2=[ zeros(1,8-length(n2)),x2]; x=x1+x2; figure(3) stem(n,x)

axis([-1 8 0 2.1])

已知f1(t)=sinwt , f2(t)=sin8wt , w=2pi , 求f1(t)+f2(t)和f1(t)f2(t) 的波形图

3、 信号相乘

信号序列1()x n 和2()x n 相乘所得信号()x n 的表达式为：12()()()x n x n x n =

这是样本与样本之间的点乘运算，在MA TLAB 中可采用“.*”来实现，但是在信号序列相乘之前，应对其做与相加运算一样的操作。 例3.3：用MA TLAB 实现上例中两序列相乘 clear all; n1=0:3

x1=[2 0.5 0.9 1]; figure(1) stem(n1,x1)

axis([-1 8 0 2.1] ) n2=0:7

x2=[ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7]; figure(2) stem(n2,x2)

axis([-1 8 0 0.8] ) n=0:7;

x1=[x1 zeros(1,8-length(n1))]; x2=[ zeros(1,8-length(n2)),x2]; x=x1.*x2; figure(3) stem(n,x)

axis([-1 8 0 0.35])

4、 信号翻转

信号翻转的表达式为：y(n)=x(-n)，在MATLAB 中可以用fliplr 函数实现此操作 例3.4：用MA TLAB 实现“信号相加”中的1()x n 序列翻转 clear all; n=0:3

x1=[2 0.5 0.9 1]; x=fliplr(x1); stem(n,x)

axis([-1 4 0 2.1] )

5、 信号和

对于N 点信号()x n ，其和的定义为：1

()N

n y x n ==

∑

例3.5：用MA TLAB 实现“信号相加”中的1()x n 序列和 clear all; n=0:3

x1=[2 0.5 0.9 1]; x=sum(x1)

6、 信号积

对于N 点信号()x n ，其积的定义为：1

()N

n y x n ==

∏

例3.5：用MA TLAB 实现“信号相加”中的1()x n 序列积 clear all; n=0:3

x1=[2 0.5 0.9 1]; x=prod(x1) 7 卷积

1、 完成)(1t f 与)(2t f 两函数的卷积运算 其中：)4()()(),

()(221--==-t u t u t f t u e

t f t

在一个图形窗口中，画出)(1t f 、

)(2t f 以及卷积结果。要求每个坐标系有标题、坐标轴名称。 p=0.1; t=0:p:10;

f1=exp(-2*t).*u(t); f2=u(t)-u(t-4); f=conv(f1,f2); subplot(1,3,1); plot(t,f1,'r');

title('f1(t)=e^-2*t*u(t)'); xlabel('t(sec)'); ylabel('f1(t)'); subplot(1,3,2); plot(t,f2,'g');

title('f2(t)=u(t)-u(t-4)'); xlabel('t(sec)'); ylabel('f2(t)'); subplot(1,3,3); plot(f);

title('f(t)=f1(t)*f2(t)'); xlabel('t(sec)'); ylabel('f(t)');

四 MATLAB 在信号与系统中的应用

该实验用MATLAB 中库函数，如tf2zp(b,a)，ss2zp(A,B,C,D)，zplane(z,p)，freqz(b,a)等。例如：

1．传递函数为

14.02

5.0)(2

2+++-=s s s s s H ，求其零极点图。 程序如下：

num=[1 0.5 2]； 分子系数，按降幂顺序排列 den=[1 0.4 1]； 分母系数，按降幂顺序排列 [z,p]=tf2zp(num,den)； 用tf2zp 函数求出其零点z 和极点p zplane(z,p) 作出零极点图

2．若给出的是滤波器的输入与输出的状态方程，如：

u x x ??????+??????-='013101，[]u x y 0141+???

???-=，求其零极点图。 程序如下： A=[1，0；1，-3]；

B=[1；0]；

C=[-41

，1]；

D=0；

[z,p]=ss2zp(A,B,C,D)；求出其零极点z,p zplane(z,p)

在连续时间系统中，当极点在虚轴的右半平面时，系统不稳定，在虚轴上， 单阶极点系统稳定；若零点均处于左半平面内，则系统为最小相位系统。在离散系统中，极点在单位圆外系统不稳定，在单位圆上，单阶极点系统稳定；零点在单位圆内，系统为最小相位系统。

对一滤波器，我们不仅要知道它的零点和极点，还要了解它的频率特性，本实验可求其频率特性。

对模拟滤波器，可用freqs 函数求得其频率特性，对数字滤波器，则用freqz 函数求得。

3 已知模拟滤波器的传递函数为

14.01

3.02.0)(2

2++++=s s s s s H ，求其频率特性。 程序如下：

num=[0.2 0.3 1]; den=[1 0.4 1];

w=logspace(-1,1); 频率范围 freqs(num,den,w)

例4．数字滤波器212

14.013.02.0)(----++++=

z z z z z H ，取样点数为128点，求其频率特性。

程序如下：

num=[0.2 0.3 1]; den=[1 0.4 1]; freqz(num,den,128) 实验内容:

1．已知下列传递函数H(s)或H(z)，求其零极点，并画出零极点图。

a. )2)(1()2)(1(3)(++--=s s s s s H

b.

s s H 1)(=

c.

521

)(2

2+++=s s s s H d.

1235

2236.0)(23

23++++++?=z z z z z z z H

2．已知下列H(s)或H(z)，求其频响。 a. 111)(-+=

z z H

b.

122)(2++=

s s s s H

c.

22

161.01)1()(--++=

z z z H d. )2)(1()2)(1(3)(++--=

s s s s s H

1： 若某连续系统的输入为e (t )，输出为r (t )，系统的微分方程为：

''()5'()6()3'()2()y t y t y t f t f t ++=+

①求该系统的单位冲激响应h (t )及其单位阶跃响应g (t )。 ②若2()()t f t e t ε-= 求出系统的零状态响应y(t )

分析： ① 求冲激响应及阶跃响应的MATLAB 程序：

a=[1 5 6];b=[3 2];

subplot(2,1,1), impulse(b,a,4) subplot(2,1,2), step(b,a,4)

运行结果如右：

② 求零状态响应的MATLAB 程序：

a=[1 5 6];b=[3 2];

p1=0.01; %定义取样时间间隔为0.01 t1=0:p1:5; %定义时间范围 x1=exp(-2*t1); %定义输入信号

lsim(b,a,x1,t1), %对取样间隔为0.01时系统响应进行仿真

hold on; %保持图形窗口以便能在同一窗口中绘制多条曲线 p2=0.5; %定义取样间隔为0.5 t2=0:p2:5; %定义时间范围 x2=exp(-2*t2); %定义输入信号

lsim(b,a,x2,t2), hold off %对取样间隔为0.5时系统响应进行仿真并解除保持

运行结果如下：

2 已知一个过阻尼二阶系统的状态方程和输出方程分别为：

010'()()()232x t X t f t ????=+????--????

， r (t )=[0 1]X (t ) 。

若系统初始状态为X (0)=[4 -5]T ， 求系统在4()3()t f t e t ε-=作用下的全响应。 求全响应程序如下：

A=[0 1 ; -2 -3] ;B=[0 2]';C=[0 1];D=[0];

X0=[4 -5]'; %定义系统初始状态 t=0: 0.01:10;

E =[3*exp(-4*t).*ones(size(t))]'; %定义系统激励信号

[r , x]=lsim(A,B,C,D,E,t,X0); %求出系统全响应的数值解 plot(t,r) %绘制系统全响应波形 运行结果如右。

3 已知描述离散系统的差分方程为：()0.25(1)0.5(2)()(1)y k y k y k f k f k --+-=+-，

且已知系统输入序列为12()()()

k

f k k ε=， ① 求出系统的单位函数响应h (k )在-3 ~10离散时间范围内响应波形。

② 求出系统零状态响应在0 ~15区间上的样值；并画出输入序列的时域波形以及系统零状态响应的波形

分析：①求系统的单位函数响应的MATLAB 程序：

a=[1,-0.25,0.5]; b=[1,1,0];

impz(b,a,-3:10), title('单位响应') %绘出单位函数响应在-3 ~10区间上的波形

运行结果如图a 。

②求零状态响应的MATLAB 程序：

a=[1,-0.25,0.5];b=[1,1,0]

k=0:15; %定义输入序列取值范围 x=(1/2).^k; %定义输入序列表达式

y=filter(b,a,x) %求解零状态响应样值 subplot(2,1,1),stem(k,x) %绘制输入序列的波形 title('输入序列')

subplot(2,1,2),stem(k,y) %绘制零状态响应的波形 title('输出序列')

运行结果如下：

y =

Columns 1 through 10

1.0000 1.7500 0.6875 -0.3281 -0.2383 0.1982 0.2156 -0.0218 -0.1015 -0.0086 Columns 11 through 16

图a. ①运行结果 图b. ②运行结果

实验内容

1. 已知描述系统的微分方程和激励信号e (t ) 分别如下，并用MATLAB 绘出系统单位冲激响应和系统零状态响应的波形.

①''()4'()4()'()3()y t y t y t f t f t ++=+；()()t f t e t ε-= ②''()2'()26()'()y t y t y t f t ++=；()()f t t ε= ③''()4'()3()()y t y t y t f t ++=；2()()t f t e t ε-=

④如下图所示的电路中，已知1234()R R R ===Ω，121()L L H ==，且两电感上初始电流分别为12(0)2(),(0)0()i A i A ==，如果以电阻3R 上电压()y t 作为系统输出，请求出系统在激励()12()f t t ε=（v ）作用下的全响应。

2. 请用MATLAB 分别求出下列差分方程所描述的离散系统，在0~20时间范围内的单位函数响应、阶跃响应和系统零状态响应的数值解，并绘出其波形。另外，请将理论值与MATLAB 仿真结果在对应点上的值作比较，并说出两者的区别和产生误差的原因。

① ()2(1)(2)()y k y k y k f k +-+-=；14

()()f k k ε=

② (2)0.7(1)0.1()7(2)2(1)y k y k y k f k f k +-++=+-+；()()f k k ε=

③ 51

66()(1)(2)()(2)

y k y k y k f k f k --+-=--；()()f k k ε= ④一带通滤波器可由下列差分方程描述：()0.81(2)()(2)y k y k f k f k +-=--， 其中()f k 为系统输入， ()y k 为系统输出。请求出当激励为

[]()1010cos(/2)10cos()()f k kn kn k ε=++（选取适当的n 值）时滤波器的稳态输出。

（四）、傅里叶变换的MATLAB 实现

①在MATLAB 中实现傅里叶变换的函数为：

F=fourier( f ) 对f(t)进行傅里叶变换，其结果为F(w) F =fourier(f,v) 对f(t)进行傅里叶变换，其结果为F(v) F=fourier( f,u,v ) 对f(u)进行傅里叶变换，其结果为F(v) ②傅里叶反变换

f=ifourier( F ) 对F(w)进行傅里叶反变换，其结果为f(x) f=ifourier(F,U) 对F(w)进行傅里叶反变换，其结果为f(u) f=ifourier( F,v,u ) 对F(v)进行傅里叶反变换，其结果为f(u) 例① 求门函数()(1)(1)f t t t εε=+--的傅里叶变换，并画出幅度频谱图

MATLAB 程序如下：

syms t w %定义两个符号变量t,w Gt=sym('Heaviside(t+1)-Heaviside(t-1)'); %产生门宽为2的门函数

Fw=fourier(Gt,t,w); %对门函数作傅氏变换求F(jw)

FFw=maple('convert',Fw,'piecewise'); %数据类型转换，转为分段函数，此处可以去掉

FFP=abs(FFw); %求振幅频谱| F(jw)|

ezplot(FFP,[-10*pi 10*pi]);grid; %绘制函数图形，并加网格 axis([-10*pi 10*pi 0 2.2]) %限定坐标轴范围

运行结果：Fw= exp(i*w)*(pi*Dirac(w)-i/w)-exp(-i*w)*(pi*Dirac(w)-i/w)

% Dirac(w)为δ(ω)，即傅立叶变换结果中含有奇异函数，故绘图前需作函数类型

转换

FFw= -i*exp(i*w)/w+i*exp(-i*w)/w % FFw 为复数

FFP= abs(-i*exp(i*w)/w+i*exp(-i*w)/w) %求FFw 的模值

例② 求函数2

1

()1F j ωω

=

+的傅里叶反变换f (t ) MATLAB 程序如下：

syms t w %定义两个符号变量t,w Fw=sym('1/(1+w^2)'); %定义频谱函数F(jw)

ft=ifourier(Fw,w,t); %对频谱函数F(jw)进行傅氏反变换

运行结果： ft =

1/2*exp(-t)*Heaviside(t)+1/2*exp(t)*Heaviside(-t)

2、傅里叶变换的数值计算实现法

严格说来，如果不使用symbolic 工具箱，是不能分析连续时间信号的。采用数值计算方法实现连续时间信号的傅里叶变换，实质上只是借助于MATLAB 的强大数值计算功能，特别是其强大的矩阵运算能力而进行的一种近似计算。傅里叶变换的数值计算实现法的原理如下：

对于连续时间信号f(t)，其傅里叶变换为：

()F j ω0

()lim

()j t

j n n f t e

dt f n e ωωττττ∞

∞

---∞

→=-∞

==∑

?

其中τ为取样间隔，如果f(t)是时限信号，或者当|t|大于某个给定值时，f(t)的值已经衰减得很厉害，可以近似地看成是时限信号，则上式中的n 取值就是有限的，假定为N ，有： ()F j ω1

()N j n n f n e

ωτ

τ

τ--==∑

若对频率变量ω进行取样，得：

()()k F k F j ω=1

()0k N j n n f n e k M ωτ

ττ--==<<∑

通常取：

2k k k M

M ωπ

ωτ

==

，其中0ω是要取的频率范围，或信号的频带宽度。采用MATLAB 实现上式时，其要点是要生成f(t)的N 个样本值()f n τ的向量，以及向量k j n e

ωτ

-，两向量的

内积（即两矩阵的乘积），结果即完成上式的傅里叶变换的数值计算。

注意：时间取样间隔τ的确定，其依据是τ必须小于奈奎斯特（Nyquist ）取样间隔。如果f(t)不是严格的带限信号，则可以根据实际计算的精度要求来确定一个适当的频率0ω为信号的带宽。

例③ 用数值计算法实现上面门函数()(1)(1)f t t t εε=+--的傅里叶变换，并画出幅度频谱图.

分析： 该信号的频谱为()2()F j Sa ωω=，其第一个过零点频率为π，一般将此频率认为是信号的带宽。但考虑到()F j ω的形状（为抽样函数），假如将精度提高到该值的50倍，即取05050B ωωπ==，则据此确定的Nyquist 取样间隔为：0011

0.02222F τωπ

≤

==?。

MATLAB 程序如下：

R=0.02; %取样间隔τ=0.02 t=-2:R:2; % t 为从-2到2，间隔为0.02的行向量,有201个样本点 ft=[zeros(1,50),ones(1,101),zeros(1,50)]; % 产生f(t)的样值矩阵（即f(t)的样本值组

成的行向量）

W1=10*pi; %取要计算的频率范围

M=500; k=0:M; w=k*W1/M; %频域采样数为M, w 为频率正半轴的采样点 Fw=ft*exp(-j*t'*w)*R; %求傅氏变换F(jw) FRw=abs(Fw); %取振幅

W=[-fliplr(w),w(2:501)] ; %由信号双边频谱的偶对称性，利用fliplr(w)

形成负半轴的点，% w(2:501)为正半轴的点，函数fliplr(w)对矩阵w 行向量作180度反转

FW=[fliplr(FRw),FRw(2:501)]; %形成对应于2M+1个频率点的值

Subplot(2,1,1) ; plot(t,ft) ;grid; %画出原时间函数f(t)的波形，并加网格 xlabel('t') ; ylabel('f(t)'); %坐标轴标注 title('f(t)=u(t+1)-u(t-1)'); %文本标注

subplot(2,1,2) ; plot(W,FW) ;grid on; %画出振幅频谱的波形,并加网格 xlabel ('W') ; ylabel ('F(W)'); %坐标轴标注 title('f(t)的振幅频谱图'); %文本标注

运行结果如下：

三、 实验内容

1.编程实现求下列信号的幅度频谱

(1) 求出1()(21)(21)f t t t εε=+--的频谱函数F 1(j ω)，请将它与上面门宽为2的门

函数()(1)(1)f t t t εε=+--的频谱进行比较，观察两者的特点，说明两者的关系。

(2) 三角脉冲 21||||1

()0

||1

t t f t t -≤?=?

>? (3) 单边指数信号3()()t f t e t ε-= (4) 高斯信号2

3()t f t e -=

2.利用ifourier( ) 函数求下列频谱函数的傅氏反变换

(1)22()16F j j ωωω=-+ (2) 22

()58

()()65

j j F j j j ωωωωω+-=++

例４：若()s i n ()4

n x n π

=是一个N=32的有限序列，利用MATLAB 计算它的DFT 并画出图形。

N=32; n=0:N-1;

xn=cos(pi*n/6); k=0:N-1;

WN=exp(-j*2*pi/N); nk=n ’*k;

WNnk=WN.^nk; Xk=xn*WNnk; figure(1)

stem(n,xn) figure(2)

stem(k,abs(Xk))

在MATLAB 中，可以直接利用内部函数fft 来实现FFT 算法，该函数是机器语言，而不是MATLAB 指令写成的，执行速度很快。常用格式为： y=fft(x) y=fft(x,N)

4．2 )4/sin()8/sin()(ππn n n x +=是一个N ＝16的有限序列，用MATLAB 求其DFT 的结果，并画出其结果图，如图3－1所示。

图 3－1 有限长序列的DFT 结果图

程序 N=16;

n=0:1:N-1; %时域采样 xn=sin(n*pi/8)+sin(n*pi/4); k=0:1:N-1; %频域采样 WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^nk;

Xk=xn*WNnk;

subplot(2,1,1)

stem(n,xn);

subplot(2,1,2)

stem(k,abs(Xk));

n

j

(3/π

=，n=[0,10]

)

n

9.0(

e

x)

运算结果

Xk =

Columns 1 through 5

0.0000 -0.0000 - 8.0000i -0.0000 - 8.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i

Columns 6 through 10

-0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i

Columns 11 through 15

0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 8.0000i

Column 16

0.0000 + 8.0000i

用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。

解 MATLAB程序如下：

a=[-2 0 1 -1 3];

b=[1 2 0 -1];

c=conv(a,b);

M=length(c)-1;

n=0:1:M;

stem(n,c);

xlabel('n'); ylabel('幅度');

图1.1给出了卷积结果的图形，求得的结果存放在数组c中为：{-2 -4 1 3 1 5 1 -3}。

1-2 用MATLAB计算差分方程

当输入序列为时的输出结果。

解 MATLAB程序如下：

N=41;

a=[0.8 -0.44 0.36 0.22];

b=[1 0.7 -0.45 -0.6];

x=[1 zeros(1,N-1)];

k=0:1:N-1;

y=filter(a,b,x);

stem(k,y)

xlabel('n');ylabel('幅度')

图 1.2 给出了该差分方程的前41个样点的输出，即该系统的单位脉冲响应。

用MATLAB 计算例1-2差分方程

所对应的系统函数的DTFT 。

解 例1-2差分方程所对应的系统函数为：

123

123

0.80.440.360.02()10.70.450.6z z z H z z z z -------++=

+--

其DTFT 为

23230.80.440.360.02()10.70.450.6j j j j j j j e e e H e

e e e ωωωω

ωωω--------++=

+--

用MATLAB 计算的程序如下： k=256;

num=[0.8 -0.44 0.36 0.02]; den=[1 0.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);