,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,谱分析与谱估计,目录,引言,功率谱估计,功率谱密度,Matlab,应用,谱分析与谱估计目录 引言,1.,引言,声纳系统,搜寻水面舰艇或潜艇,机械主件及部件,故障诊断,语音识别,语音频带压缩,2,1.引言声纳系统搜寻水面舰艇或潜艇2,2.,功率谱估计,相关函数法,周期图法,目标,获得实随机过程 的功率谱密度 的,近似,估计,3,2.功率谱估计相关函数法3,2.1,相关函数法,Blackman-Tukey,算法，,BT,算法,基本思路：从时域上先求信号自相关函数，再做,Fourier,变换，求得功率谱估计值。,自相关序列,4,2.1 相关函数法Blackman-Tukey算法，BT算法,Wiener-Khinchin,公式,弱平稳随机过程的功率谱密度是其相应自相关函数的,Fourier,变换,估计方法分为两种,直接估计法（非参数方法）,依赖于信号产生模型的方法（参数化方法）,简单快速、不够精确,计算相对复杂、但更精确,5,Wiener-Khinchin公式简单快速、不够精确计算相对,2.2,直接估计法的局限性,噪声的不利影响；,对于随机过程，仅有一组可用的信号（实现）；,信号的长度有限，如 。,6,2.2 直接估计法的局限性噪声的不利影响；6,2.3,假设,对于随机过程的实现,各态历经,(Ergodic),：统计平均可以用算术平均来代替；,平稳性,(Stationary),：信号的无限平均可以用有限平均来代替；,以上两种平均均可从 求出。,7,2.3 假设 对于随机过程的实现7,2.4,窗函数截断,根据以上假设，可以对数据进行截断，利用窗函数得到以下的新信号 ：,为窗函数，从信号中截取一部分信号,.,8,2.4 窗函数截断根据以上假设，可以对数据进行截断，利用窗函,3.,周期图（,Periodogram,）,功率谱密度计算公式为,很明显当,可以利用,FFT,有效计算。,9,3.周期图（Periodogram）功率谱密度计算公式,对其做,FFT,，有,其估计值为,10,对其做FFT，有10,由于序列,x,(,n,),的离散傅里叶变换,X,(,k,),具有周期函数性质，,因此称为长度为,N,的实平稳随机信号序列的周期图,11,由于序列x(n)的离散傅里叶变换X(k)具有周期函数性质，1,有限自相关,也称为短信号的自相关序列，这是将在功率谱密度估计中要用到的参数。,12,有限自相关12,3.1 PSD,估计,周期图法计算式可改写为,周期图法的有限自相关表达方式；存在问题：,统计变异性,估计的偏度误差,13,3.1 PSD估计周期图法计算式可改写为13,3.2,偏差与方差,偏差（,Bias,）,当采样个数,N,趋于无穷时，估计值是否收敛于真实值？,是,无偏：,Unbiased estimate;,否,有偏：,Biased estimate;,14,3.2 偏差与方差偏差（Bias）14,偏差分析,PSD,估计式中截尾信号的自相关序列为：,15,偏差分析PSD估计式中截尾信号的自相关序列为：15,根据卷积定理，有,偏差定义为谱密度真实均值与期望均值之差,16,根据卷积定理，有16,例：,3.3,选用矩形窗函数,与真实功率谱密度进行卷积运算时，得到的是平均周期图,(,平滑,PSD),。,矩形窗的主体宽度为 ，因此当,时，有,17,例：3.3选用矩形窗函数17,因此 是真实功率谱密度的渐近无偏估计。对该结论进行推广可以得到对窗函数的一些具体要求：,标准化条件：,窗口的主体部分必须随,1/N,递减。,18,因此 是真实功率谱密度的渐近无偏估计,方差分析,1.,谱估计的变异性,变异系数：谱估计的均值和方差之比,利用周期图法做谱估计时，周期图,X,(,k,),为复数，因此有,19,方差分析1.谱估计的变异性19,考虑高斯变量情形，则,X,R,(,k,),和,X,I,(,k,),也是高斯变量；,卡方分布公式：,20,考虑高斯变量情形，则XR(k)和XI(k)也是高斯变量；20,因此有,因此变异系数为,根据周期图公式，自由度,n=2,，因此变异系数为,1,，相对误差达到,100%,，估计极其不准确。,21,因此有21,2.,平均化处理,对周期图进行平均化处理可以减小变异系数，得到较高的谱估计精度。,将序列,x(n),分段，求各段周期图，再做平均，此时有：,各段频率分量的实部与虚部互相独立,22,2.平均化处理22,此时卡方分布自由度为,n=2q.,变异系数为 ，表明随着平均谱,分段数,q,增大，变异系数减小，从而可以通,过增加分段数来减小变异性。,对连续随机过程，样本总体长度为,T,(,T,=,N,t,),t,为采样间隔，分段长度,T,e,(T/q),，那么分析带宽,B,e,=1/,Te,，因此有,23,此时卡方分布自由度为n=2q.23,4.Matlab,应用举例,求卷积,x=randn(1,100);,w=10;,y=conv(ones(1,w)/w,x);,avgs=y(10:99);,plot(avgs),24,4.Matlab应用举例求卷积24,25,25,Ensemble average,w=10;,for i=1:w;,X(i,:)=randn(1,100);,end,AVGS=mean(X);,plot(AVGS),26,Ensemble average26,27,27,交叉相关,xcorr,交叉协方差,xcorr,两组信号的交叉相关等价于两组信号的卷积,(,其中一组逆序,fliplr),xcov,的作用是在计算交叉相关之前去掉输入的均值；,信号与自身的交叉相关称为自相关。,28,交叉相关xcorr,交叉协方差xcorr28,x=-1 0 1;,y=0 1 2;,xcorr(x,y),conv(x,fliplr(y),xcov(x,y),xcov(x,x),xcov(x,y-1),29,x=-1 0 1;29,例,1.,交叉相关,考虑简单的目标测距系统，输出信号,x,与返回信号,y,之间的关系为：,其中,alpha,为衰减因子，,d,为时滞，,beta,为信道噪声，若,T,为信号返回时间，则,x,，,y,将在,n=T,时相关。目标则可以定位在,vT,的距离，其中,v,为信号的传播速率。,30,例1.交叉相关考虑简单的目标测距系统，输出信号x与返回信号y,x=zeros(1,25),1,zeros(1,25);,subplot(311),stem(x),y=0.75*zeros(1,20),x+0.1*randn(1,71);,subplot(312),stem(y),c lags=xcorr(x,y);,subplot(313),stem(lags,c),31,x=zeros(1,25),1,zeros(1,25,32,32,进行,DFT,delta1=1 zeros(1,11);,fftgui(delta1),delta2=0 1 zeros(1,10);,fftgui(delta2),deltaNyq=zeros(1,6),1,zeros(1,5);,fftgui(deltaNyq),square=zeros(1,4),ones(1,4),zeros(1,4);,fftgui(square),t=linspace(0,1,50);,periodic=sin(2*pi*t);,fftgui(periodic),33,进行DFT33,34,34,FFT Demo,太阳黑子,sigdemo1,playshow fftdemo,phone,playshow sunspots,35,FFT Demo太阳黑子35,功率谱分析,Fs=100;,t=0:1/Fs:10;,y=sin(2*pi*15*t)+sin(2*pi*30*t);,nfft=512;,Y=fft(y,nfft);,f=Fs*(0:nfft-1)/nfft;,Power=Y.*conj(Y)/nfft;,plot(f,Power),title(Periodogram),36,功率谱分析Fs=100;36,37,37,figure,ryy=xcorr(y,y);,Ryy=fft(ryy,512);,plot(f,abs(Ryy),title(DFT of Autocorrelation),38,figure38,39,39,t=0:1/100:10-1/100;,x=sin(2*pi*15*t)+sin(2*pi*30*t);,periodogram(x,512,100);,figure,pwelch(x,512,100);,figure,pmtm(x,512,100);,40,t=0:1/100:10-1/100;40,41,41,42,42,SPTool,43,SPTool43,