实验4 离散信号的时频描述与分析

一、实验目的

1、熟悉常见离散信号的意义、特性及波形。 2、学会用MATLAB表示常用离散信号的方法。 3、学会运用Matlab进行离散时间信号的基本运算。 二、实验环境

计算机，Matlab软件 三、实验原理

离散信号的绘制一般用stem函数，MATLAB只能表示一定时间范围内有限长度的序列。而对于无限长序列，只能在一定范围内表示出来。对于任意离散序列f[k]，需用两个向量来表示：一个表示k的取值范围，另一个表示序列的值。

例如序列f[k]={2,1,1,-1,3,0,2}（这里的“1”表示k=0所对应的位置，下同）可用MATLAB表示为： k=-2:4; f=[2,1,1,-1,3,0,2];

若序列是从k=0开始的，则只用一个向量f就可以表示该序列了。由于计算机内存的限制，MATLAB无法表示一个任意的无穷序列。

在MATLAB中，绘制离散序列波形图的专用命令为stem()，其格式有： stem(k,f)：在图形窗口中，绘制出样值顶部为空心圆的序列波形图。

stem(k,f,’fill’)：在图形窗口中，绘制出样值顶部为实心圆的序列波形图。 1、常用信号的MATLAB表示 1）指数序列

离散指数序列的一般形式为ak，可以用MATLAB中的数组幂运算a.^k来实现。 指数衰减序列的MATLAB源程序如下（取A=1，a=-0.6）: k=0:10;A=1;a=-0.6; fk=A*a.^k; stem(k, fk);

2）正弦序列

离散正弦序列的MATLAB表示与连续信号类似，只不过是用stem函数而不是用plot函数来画出序列的波形。正弦序列sin()k的MATLAB源程序如下：

k=0:39;

fk=sin(pi/6*k); stem(k,fk);

3）单位脉冲序列 单位冲击序列定义为：

(k) k01

k00 一种简单的表示方法是借助MATLAB中的全零矩阵函数zeros。全零矩阵zeros（1，N）

产生一个由N个零组成的行向量，对于有限区间的可以在MATLAB中表示为：

k=-50:50;

delta=[zeros(1,50),1,zeros(1,50)];%用数组的形式构成。 stem(k,delta);

另一种更有效的表示方法是将单位冲激序列写成MATLAB函数的形式，利用关系运算符“= =”来实现。单位冲激序列在范围内，用MATLAB函数形式可写为：

function[f,k]=impseq(k0,k1,k2) %产生f[k]=delta(k-k0);k1<=k<=k2 k=[k1:k2]; f=[(k-k0)==0]; 程序中关系运算“（k-k0）==0”的结果是一个由“0”和“1”组成的向量，即k=k0时返回True值1，k≠k0时False返回值0。

4）单位阶跃序列 k=-50:50;

uk=[zeros(1,50),ones(1,51)]; stem(k,uk);

5）随机序列

x=rand(1,20); stem(x);

xlabel('Time indexn');ylabel('Amplitude'); title('Unit Sample Sequence');

2、离散信号的基本运算 1）信号的相加、相乘

信号相加： x=x1+x2; 信号相乘： x=x1.*x2;

注意：x1和x2序列应该具有相同的长度，位置对应，才能相加。

w0=pi/20;a=1.05;

n1=[-20:20];n2=[-20:20];

subplot(2,2,1);stem(n1,sin(w0.*n1),'.k');title('x1(n)');axis([-20,20,-1,3]); subplot(2,2,2);stem(n2,a.^n2,'.k');title('x2(n)');axis([-20,20,-1,3]);

n=[min(min(n1),min(n2));max(max(n1),max(n2))];%计算两个信号长度的区间范围。

x1=zeros(1,length(n));%保证两个信号的长度相等。 x2=zeros(1,length(n)); x1=sin(w0.*n1); x2=a.^n2; y1=x1+x2; y2=x1.*x2;

subplot(2,2,3);stem(n1,y1,'.k');title('y1(n)');axis([-20,20,-1,3]); subplot(2,2,4);stem(n2,y2,'.k');title('y2(n)');axis([-20,20,-1,3]);

2）信号的翻转：调用格式：y=fliplr(x)

x1=[1,2,3,4,2,3,4,5]; y1=x1;

y2=fliplr(x1) ;

subplot(2,1,2);stem(x1,y1,'.r');title('x(-n)'); subplot(2,1,1);stem(x1,y2,'.r');title('x(n)');

3）平移：（同连续信号）

4）求和：调用格式：y=sum(x(k1:k2)) 5）差分：调用格式：y=diff(x)

6）卷积和:两个有限长序列f1，f2卷积可调用MATLAB函数conv，调用格式是f=conv(f1,f2)。f的长度等于f1和f2长度之和减一， f的起点是f1和f2的起点之和，f的终点是f1和f2的终点之和。

3、DFT

在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，它的DFT为：

有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等

距采样，因此可以用于序列的谱分析。 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;

补充：圆周移位的源程序： function yn=cirshift(x,m,N); if length(x)>N

error('N must be > tht length of x'); end

x=[x zeros(1,N-length(x))]; n=[0:1:N-1];ss

n=mod(n-m,N)' yn=x(n+1);

4、快速FFT

在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，它的DFT和反变换定义分别为：

MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有fft、Ifft、fft2 、ifft2, fftn、ifftn和fftshift、Ifftshift等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。

fft和ifft函数 一维快速正傅里叶变换和逆傅里叶变换 Y=fft(X)

参数说明:·如果X是向量，则采用傅立叶变换来求解X的离散傅立叶变换； ·如果X是矩阵，则计算该矩阵每一列的离散傅立叶变换；

·如果X是（ND)维数组，则是对第一个非单元素的维进行离散傅立叶变换； Y＝fft(X,N)

参数说明:N是进行离散傅立叶变换的X的数据长度，可以通过对X进行补零或截取来实现。

fft2和ifft2函数 调用方式 Y=fft2(X)

参数说明:·如果X是向量，则此傅立叶变换即变成一维傅立叶变换fft； ·如果X是矩阵，则是计算该矩阵的二维快速傅立叶变换；

·数据二维傅立叶变换fft2（X）相当于fft(fft(X))，即先对X的列做一维傅立叶变换，然后再对变换结果的行做一维傅立叶变换。 Y＝fft2(X,M,N)

通过对X进行补零或截断，使得成为（M*N）的矩阵。 函数Ifft2的参数应用与函数Fft2完全相同。

''5、利用FFT求解离散傅立叶变换（IFFT） 逆变换IFFT为

1x(n)NX(n)Wk0N1nkN r0,1,N1 2求X（k）的逆FFT可以分为以下3个步骤： （1）取X（k）的共扼，得到X*(k)； （2）求X*(k) 的FFT，得Nx*(k)；

（3）取x*(n) 的共扼，并除以N，即得到x(n)。

采用这种方法可以直接用FFT程序来计算逆FFT。有关IFFT的具体应用，与FFT一一对应，在此不再赘述。

四、实验内容

1、上机前，认真阅读实验原理，掌握信号表示和信号运算的方法，并对原理中的程序进行验证。

2、利用MATLAB命令画出下列离散信号的波形图。 1）()ku[k]

12

k2）(0.8)cos(0.9πk)

3）u[k2]u[k5]

3、绘出信号x(k)1.5sin(2*0.1k)的频率是多少？周期是多少？产生一个数字频率为0.9的正弦序列，并显示该信号，说明其周期?

根据图可以看出：x(k)=1.5sin(2*pi*0.1*k)是周期信号，周期为10，而数字频率为0.9的正弦序列不是周期信号，对于x(k)，N=2*pi/(0.1*pi)=10，第二个0.9不是pi的倍数，所以不是周期的。

4、对连续的单一频率周期信号x(n)=sin⁡(2πn𝑓𝑎)按采样频率𝑓𝑠=8𝑓𝑎 采样，截取长度N分别选N =20和N =16，用函数fft用于计算该信号的离散傅里叶变换DFT，并观察其DFT结果的幅度谱。

5、分别对下列序列进行频谱分析：编制DFT程序及FFT程序，并比较DFT程序与FFT程序的运行时间。

(1) 实指数序列(1.08)n

在程序两端加上tic和toc判断时间

时间对比：DFT时间已过 0.046848 秒。 FFT时间已过 0.009625 秒。

(2) 周期为N的正弦序列sin(

2n),0nN1N

时间对比：DFT时间已过 0.046568 秒。FFT时间已过 0.009999 秒。

6、序列x(n)2sin(0.48n)cos(0.52n),0n100，试绘制x(n)及它的离散傅里叶变换X(k)图。

五、思考题

1、对于周期序列，如果周期不知道，如何用FFT进行谱分析？

答：周期信号的周期预先不知道时，可先截取M点进行DFT,再将截取长度扩大1倍截取，比较结果，如果二者的差别满足分析误差要求，则可以近似表示该信号的频谱，如果不满足误差，要求就继续将截取长度加倍，重复比较，直到结果满足要求。