谱分析与谱估计.ppt

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