实验二 信号分析

实验二信号分析1

1.熟悉离散信号和系统的时域特性；2.熟悉线性卷积和相关的计算编程方法；3.掌握序列付氏变换的计算机实现方法，利用序列付氏变换对离散信号、系统和系统的响应进行频域分析；4.利用MATLAB对wav文件进行频谱分析。一、实验目的

2二、实验原理1.离散时间与系统(1)离散时间信号信号可以粗略地分模拟信号和数字信号两大类，模拟信号用x(t)表示，其中变量t可以表示任何物理量，一般假定它代表以秒为单位的时间；离散信号用x(n)表示，其中变量n为整数并代表时间的离散时刻。在工程实际中，一个序列一般如下表示：3二、实验原理x(n)={x(n)}={…,x(-1),x(0),x(1),…}在MATLAB中，一般用一个列向量来表示有限长度的序列。然而这样一个向量并没有包含采样位置的信息，因此序列x(n)一般用两个向量来表示，一个表示采样时刻，一个表示采样值，例如x(n)={2,6,1,2,0,3,4,5,6}就应该用如下MATLAB语句来表示：4二、实验原理n=[-4-3-2-101234];x=[261203456];

当不需要采样位置信息或这个信息是多余的时候(例如序列从n=0开始)，用户可以只使用x向量来表示序列。由于有限的内存，MATLAB无法表示一个任意无限序列。5二、实验原理[例1]编写产生单位采样序列、单位阶跃序列的函数。

IMPSEQ.Mfunction[x,n]=impseq(n0,n1,n2)%产生x(n)=delta(n-n0);n1<=n<=n2;n=[n1:n2];x=[(n-n0)==0];

STEPAEQ.Mfunction[x,n]=stepseq(n0,n1,n2)%产生x(n)=u(n-n0);n1<=n<=n2;n=[n1:n2];x=[(n-n0)>=0];6二、实验原理[例2]编写序列相加、相乘、位移和反褶的函数。SIGADD.Mfunction[y,n]=sigadd(x1,n1,x2,n2)%实现序列相加，即y(n)=x1(n)+x2(n)%n1可以不等于n2n=min(min(n1),min(n2)):max(max(n1),max(n2));%初始化y(n)的长度y1=zeros(1,length(n));y2=y1;y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;%具有y(n)长度的x1y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;%具有y(n)长度的x2y=y1+y2;7二、实验原理SIGMULT.Mfunction[y,n]=sigmult(x1,n1,x2,n2)%实现序列相加，即y(n)=x1(n)*x2(n)%n1可以不等于n2n=min(min(n1),min(n2)):max(max(n1),max(n2));%初始化y(n)的长度y1=zeros(1,length(n));y2=y1;y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;%具有y(n)长度的x1y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;%具有y(n)长度的x2y=y1.*y2;8二、实验原理SIGSHIFT.Mfunction[y,n]=sigshift(x,m,n0)%实现序列移位，即y(n)=x(n-n0)n=m+n0;y=x;9二、实验原理SIGFOLD.Mfunction[y,n]=sigfold(x,n)%实现序列反褶，即y(n)=x(-n)y=fliplr(x);n=-fliplr(n);10二、实验原理[例3]产生如下复数信号：

x(n)=e(-0.3+j0.1)n-15<=n<=15

要求绘出该信号的幅度、相位、实部和虚部的图形。n=-15:15alpha=-0.3+0.1j;x=1.5*exp(alpha*n);subplot(221)stem(n,real(x));title(‘实部’);xlabel(‘n’);11二、实验原理subplot(222)stem(n,imag(x));title(‘虚部’);xlabel(‘n’);subplot(223)stem(n,abs(x));title(‘幅度’);xlabel(‘n’);subplot(224)stem(n,(180/pi)*angle(x));title(‘相位’);xlabel(‘n’);12二、实验原理13二、实验原理(2)离散系统的卷积和相关一个离散时间系统可以用一个算子来描述，即T[.]，它取得一个序列x(n)并将它变换为另一个序列y(n)，即有下式：

y(n)=T[x(n)]

在数字信号处理中，一般认为系统将输入信号处理为输出信号。离散系统可以大致分为线性系统和非线性系统，这里只讨论线性系统。14二、实验原理1)

卷积一个线性系统一般应满足下式：y(n)=T[x(n)]=x(k)h(n-k)由数字信号处理理论可知，线性系统的脉冲响应即为h(n)，上式一般称为线性卷积，一般可以表示为：y(n)=x(n)*h(n)

如果任意序列是无限长度的，就不能用MATLAB来直接计算卷积。MATLAB提供了一个内部函数conv()来计算两个有限长度序列的卷积。conv()函数假定两个序列都从n=0开始，调用格式如下所示：y=conv(x,h)

15二、实验原理[例4]已知某系统的输入和冲激响应如下所示：

x(n)=n(u(n)-u(n-15))h(n)=(0.4)n(u(n)-u(n-10))

求系统的响应。执行下面的m文件：n=1:20;x=zeros(1,20);h=zeros(1,20);h1=zeros(1,20);x=n.*(stepseq(0,0,19)-stepseq(15,0,19));h1=stepseq(0,0,19)-stepseq(10,0,19);h=0.4.^n.*h1;y=conv(x,h)16二、实验原理n=length(y);x=[x,zeros(1,n-length(x))];h=[h,zeros(1,n-length(h))];n=1:n;subplot(311)stem(n,x);title(‘输入’);subplot(312)stem(n,h);title(‘冲激响应’);subplot(313)stem(n,y);title(‘输出’);17二、实验原理18二、实验原理2)相关在工程实际当中，随机信号是DSP处理的主要内容，为了计算随机信号的功率谱和谱密度，DSP理论引入了相关的概念。工程实际当中，定义互相关的两个信号应该满足下式：rxy(t)=y(t)*x(t)

相应地，自相关定义为下式：rxx(t)=x(t)*x(t)

如果序列具有有限的长度的话，计算序列的相关可以直按通过conv()函数来实现。另外，MATLAB还提供了xcorr()函数来计算相关，但它也只能计算从0时刻开始的两个序列的互相关，因此一般地，并不建议使用此函数。

19二、实验原理[例5]已知原型序列如下所示：

x(n)=[2,3,6,5,7,6,4]令y(n)为x(n)加入噪声干扰并位移后的序列，如下所示：y(n)=x(n-l)+v(n)

其中v(n)为具有零均值和单位方差的高斯序列。要求计算x(n)与y(n)之间的互相关。从y(n)的结构可以看出它相似于x(n-l)，因而它们之间的相关在n=l处应该显示出最强。为了测试这一点，用两个不同的噪声来验证这一点，执行下面的M文件：

20二、实验原理%噪声序列1x=[2365764];nx=[-3:3];[y,ny]=sigshift(x,nx,1);w=randn(1,length(y));nw=ny;[y,ny]=sigadd(y,ny,w,nw);[x,nx]=sigfold(x,nx);[rxy,nrxy]=conv_m(y,ny,x,nx);subplot(211)stem(nrxy,rxy);axis([-510–50250]);xlabel(‘延迟量1’);ylabel(‘rxy’);title(‘噪声序列1的互相关’);21二、实验原理%噪声序列2x=[2365764];nx=-3:3;[y,ny]=sigshift(x,nx,1);w1=randn(1,length(y));nw=ny;[y,ny]=sigadd(y,ny,w1,nw);[x,nx]=sigfold(x,nx);[rxy,nrxy]=conv_m(y,ny,x,nx);subplot(212)stem(nrxy,rxy);axis([-510–50250]);xlabel(‘延迟量2’);ylabel(‘rxy’);title(‘噪声序列2的互相关’);22二、实验原理23二、实验原理2.离散傅立叶变换由傅立叶分析可知，一个周期函数可以由其各个谐波分量的线性组合得到，这就是离散傅级数(DFS)。将其推广至有限持续时间序列，并产生一个新的变换，这就是离散傅立叶变换(DFT)。(1)离散傅立叶级数在DSP理论中定义周期序列x(n)满足下式：

x(n)=x(n+kN),n,k

其中N为序列的基本周期。由傅立叶分析可知，周期函数可由复指数的线性组合叠加得到。其频率为基本频率(在这里为2/N)的倍数(或谐波)。24二、实验原理因此周期序列x(n)可以表示为下式：

N-1x(n)=1/NX(k)ej2/Nkn,n=0,1,…

k=0

其中，{X(k),k=0,1,…}被称为离散傅立叶级数系数，由下式给出：

N-1X(k)=x(n)e-j2/Nkn,k=0,1,…

n=0

注意X(k)本身就是一个基本周期为N的(复值)周期序列，即满足下式：

X(k+N)=X(k)25二、实验原理[例7]编写实现离散傅立叶级数和逆离散傅立叶级数的函数。DFS.Mfunction[Xk]=dfs(xn,N)%计算离散傅立叶级数(DFS)的系数%Xk为在0<=k<=N-l间的DFS系数数组%xn为周期信号在0<=n<=N-l之间的一个单周期信号%N为xn的基本周期n=0:N-l;k=0:N-l;WN=exp(-j*2*pi/N);nk=n’*k;

WNnk=WN.^nk;

Xk=xn*WNnk;26二、实验原理IDFS.Mfunction[xn]=idfs(Xk,N)%计算逆离散傅立叶级数(IDFS)%xn为周期信号在0<=n<=N-l之间的一个单周期信号%Xk为在0<=k<=N-l之间的DFS系数数组%N为Xk的基本周期n=0:N-l;k=0:N-l;WN=exp(-j*2*pi/N);nk=n’*k;

WNnk=WN.^(-nk);

xn=(Xk*WNnk)/N;27二、实验原理(2)离散傅立叶变换离散傅立叶级数提供了一种对离散时间傅立叶变换作数值计算的技巧，同时它也提出了时域出现混叠的潜在问题。数学分析表明对离散时间傅立叶变换的采样，可得到一个周期信号x(n)。但实际中，大多数信号不具有周期性，它们很可能具有有限持续时间。对这些信号，在DSP理论中，可以通过定义一个周期信号，其基本形状为有限持续时间信号，然后采用此周期信号的DFS。实际上，在DSP理论中定义了一个新的变换，叫做离散傅立叶变换(DFT)，它是DFS的主周期。DFT为对任意有限持续时间序列可数值计算的傅立叶变换。

28二、实验原理离散傅立叶变换的定义如下所示;

N-1X(k)=DFT[x(n)]=x(n)WNkn,0<=k<=N-1n=0N-1x(n)=IDFT[(X(k))=X(k)WN-kn,0<=n<=N-1k=029二、实验原理[例8]编制实现DFT和IDFT的函数。DFT.Mfunction[Xk]=dft(xn,N)%计算离散傅立叶变换%Xk为在0<=n<=N-l间的DFT系数数组%xn为N点有限持续时间序列%N为DFT的长度n=0:N-l;k=0:N-l;WN=exp(-j*2*pi/N);nk=n’*k;

WNnk=WN.^nk;

Xk=xn*WNnk;30二、实验原理IDFT.Mfunction[xn]=idft(Xk,N)%计算逆离散傅立叶变换%xn为在0<=n<=N-l间的N点有限持续时间序列%Xk为在0<=k<=N-l间的DFT系数数组%N为DFT的长度n=0:N-l;k=0:N-l;WN=exp(-j*2*pi/N);nk=n’*k;

WNnk=WN.^(-nk);

xn=Xk*WNnk/N;31三、实验内容本实验将通过分析计算机中的wav文件信号来进一步讨论数字信号处理中的信号分析方法，只要用户的电脑有麦克和声卡都可以完成本实验。1.实验准备

在进行实验之前，首先介绍几个将要用到的函数：[x,fs,bits]=wavread(‘filename’)这是一个MATLAB中读取wav文件的数据的函数。其中的x表示一长串的数据，一般是两列(立体声)；fs是该wav文件在采集时用的采样频率；bits是指在进行A/D转化时用的量化位长(一般是8bits或I6bits)。32三、实验内容[d]=FFT(w,l)这是MATLAB中快速傅立叶变换(FFT)函数的一种输人输出形式。w是一列波形数据；l是指示用多少点的FFT，我们应该选择2的乘方的数(如16，128，1024等)，因为这样就可以使用优化的蝶形算法；d是频域的输出。由于FFT(DFT)的对称性，又输入的是实数，FFT的结果的复数序列是共轭反对称的，所以它们的模的大小对称，一般来说只用取一半的数据就可以了。

sound(w,fs,bits)和前面的wavread一样的参数表示，它将数列的数据通过声卡转化为声音。

33三、实验内容2.wav文件的一次性傅立叶变换下面，就开始用MATLAB分析wav文件。(l)步骤一：选择一个wav文件作为分析的对象，当然也可以用麦克录一段用户自己的声音。这里选择的是每个Windows系统都有的ding.wav(在windows目录中)一个比较单纯的声音“叮……”。(2)步骤二：

首先执行下面的语句；

34三、实验内容[x,fs,bits]=wavread(‘ding.wav’);sound(x,fs,bits);%%因为x有两列，这里用其中的一列来进行傅氏变化

y=x(:,1);%%看一下要处理的数据的大小

size(y)

运行结果如下所示：

ans=20191135三、实验内容由上面的结果可以看出将要处理的数据量很大，接着执行下面的语句：

plot(y)

运行结果如图所示。

36三、实验内容接着执行下面的语句：[k]=fft(y,32768);%%对其进行FFT，点数为32768(2^15)plot(abs(k));%%显示频域的幅值

37三、实验内容从图可以看出，它是以x=32768/2对称的。接着执行下面的语句：[m,i]=max(abs(k))

运行结果如下所示：

m=490.3636i=1171

最大值在1171的位置上，它代表的频率和时城的采样时间有关，相邻的两点之间的距离为f=l/(NTs)。N是离散傅氏变换用的点数，Ts是采样的时间，前面在读取wav文件时得到了采样频率fs=l/Ts。

38三、实验内容接着执行下面的语句：

f=1171/32768*fs

运行结果如下所示:

f=787.9807%%主频为787Hz

由图可知，除了1171一点外，在坐标4000~5000之间还有一个峰值。按着执行下面的语句：[m,i]=max(abs(k(4000:5000)))

运行结果如下所示：

m=22.6522i=1739三、实验内容从而知道该波的频率为：(4000+717)/32768*fs

即被处理文件的信号的主频率为3174Hz。这样，就完成了对该声波主要频谱的分析。(3)步骤三：既然知道了该声波的频谱，接下来就可以反演它的时域值。当然可以使用傅氏反变换，但在本例中将利用刚才分析的两个主要的冲击峰值来重构，执行下面的语句：40三、实验内容t=[0:0.0001:1];%%在一秒钟内采样一万次p=[sin(2*pi*787*t)*1+sin(2*pi*3174*t)*(22/490)]*0.18;%%包括两个正弦波的叠加(22/490)是两个波幅度的比由于没有考虑相位和其它的频谱分量，所以波形和原来的波形相去甚远，但大体的频率是没有错的。听一下重构的声波，执行下面的语句：

sound(p,10000);sound(x,fs,bits);41三、实验内容(4)结论这个实验通过剖析wav文件中的波形，使用户了解了声音波形的大致形状，知道了利用FFT来进行频谱分析，还掌握了MATLAB中处理多媒体的几个函数。本实验还有一些欠缺，一个DSP不可能将无限长的序列读入后再加工(本次试验用的FFT的点数也太大了)，它必须学会如何分段处理数据；并会连接每个段的频城结果。

42三、实验内容3.wav文件的分段傅立叶分析上一实验是把wav文件中的信号做一次性的傅立叶变换，其中有不少缺点。最重要的是语音是分音节的，应把它分段分析，而且实际运用中的DSP的FFT的点数是有限的，一般只能达到千点，这个实验是对前者的补充。首先执行下面的语句：[w,fs,bits]=wavread(‘ding.wav’);43三、实验内容(l)显示双声道波形：接下来显示双声道的波形，执行下面的语句：subplot(2,1,1);plot(w(:,1));subplot(2,1,2);plot(w(:,2));

运行后得到如图所示的结果。

44三、实验内容45三、实验内容(3)作1024点的分段FFT；

接下来作1024点的分段FFT,执行下面的语句：u=w(:,1);length(u)/1024z=zeros(20,1024);fori=1:19;z(i,:)=(fft(u(1024*(i-1)+1:1024*i),1024))’;%%用1024点的FFTend%%剩下的一部分，在后面补0，凑成1024个数据z(20,:)=fft(u(1024*19+1:length(u)),1024)’;[x,y]=meshgrid(1:20,1:200);%%显示频率幅值变化的细节plot3(x,y,abs(z(:,1:200)’));46三、实验内容47三、实验内容应该注意到第一条线不是很光滑，线上有许多的锯齿状，这是什么原因呢?仔细的观察比较时域图和频域曲线，可以得到的结论是这是发声的第一个阶段，声音从无到有，因此难免有广泛的频谱。就像人们要发声时，刚开始时总不是很准，由于气流从无到有，声音很粗糙，但人的耳朵可以将声音传回大脑，产生负反馈时，才能很好的把握声音。

48三、实验内容继续进行分析，执行下面的语句：u=w(:,2);z2=zeros(20,1024);fori=1:19;z2(i,:)=(fft(u(1024*(i-1)+1:1024*i),1024))’;end%%剩下的一部分，在后面补0，凑成1024个数据z2(20,:)=f