《谱分析与谱估计课件》由会员分享，可在线阅读，更多相关《谱分析与谱估计课件（43页珍藏版）》请在金锄头文库上搜索。
1、谱分析与谱估计谱分析与谱估计目录 引言引言 功率谱估计功率谱估计 功率谱密度功率谱密度 Matlab应用应用1.引言l l声纳系统搜寻水面舰艇或潜艇l l机械主件及部件故障诊断l l语音识别语音频带压缩2 22. 功率谱估计l l相关函数法l l周期图法l l目标l l获得实随机过程获得实随机过程 的功率谱密度的功率谱密度 的的近似近似估计估计 3 32.1 相关函数法l lBlackman-Tukey算法，BT算法l l基本思路：从时域上先求信号自相关函数，再做Fourier变换，求得功率谱估计值。l l自相关序列4 4l lWiener-Khinchin公式l l弱平稳随机过程的功率谱密度
2、是其相应自相关弱平稳随机过程的功率谱密度是其相应自相关函数的函数的FourierFourier变换变换l l估计方法分为两种l l直接估计法（非参数方法）直接估计法（非参数方法）l l依赖于信号产生模型的方法（参数化方法）依赖于信号产生模型的方法（参数化方法）简单快速、不够精确计算相对复杂、但更精确5 52.2 直接估计法的局限性l l噪声的不利影响；l l对于随机过程，仅有一组可用的信号（实现） ；l l信号的长度有限，如 。6 62.3 假设 对于随机过程的实现l l各态历经(Ergodic)：统计平均可以用算术平均来代替；l l平稳性(Stationary)：信号的无限平均可以用有限平均
3、来代替；l l以上两种平均均可从 求出。7 72.4 窗函数截断l l根据以上假设，可以对数据进行截断，利用窗函数得到以下的新信号 ：l l 为窗函数，从信号中截取一部分信号 .8 83. 周期图（ Periodogram ）l l功率谱密度计算公式为l l很明显当l l可以利用FFT有效计算。9 9l l对其做FFT，有l l其估计值为1010由于序列由于序列x(n)的离散傅里叶变换的离散傅里叶变换X(k)具有周期函数性质，具有周期函数性质，因此称为长度为因此称为长度为N的实平稳随机信号序列的周期图的实平稳随机信号序列的周期图1111l l有限自相关l l也称为短信号的自相关序列，这是将在功
4、率谱密度估计中要用到的参数。12123.1 PSD估计l l周期图法计算式可改写为l l周期图法的有限自相关表达方式；存在问题：l l统计变异性统计变异性l l估计的偏度误差估计的偏度误差13133.2 偏差与方差l l偏差（Bias）l l当采样个数当采样个数N N趋于无穷时，估计值是否收敛于趋于无穷时，估计值是否收敛于真实值？真实值？l l是是无偏：无偏：Unbiased estimate;Unbiased estimate;l l否否有偏：有偏：Biased estimate;Biased estimate;1414偏差分析l lPSD估计式中截尾信号的自相关序列为：1515l l根据卷
5、积定理，有l l偏差定义为谱密度真实均值与期望均值之差1616例：3.3l l选用矩形窗函数l l与真实功率谱密度进行卷积运算时，得到的是平均周期图(平滑PSD)。l l矩形窗的主体宽度为 ，因此当 时，有1717l l因此 是真实功率谱密度的渐近无偏估计。对该结论进行推广可以得到对窗函数的一些具体要求：l l标准化条件：标准化条件：l l窗口的主体部分必须随窗口的主体部分必须随1/N1/N递减。递减。1818方差分析l l1. 谱估计的变异性l l变异系数：谱估计的均值和方差之比变异系数：谱估计的均值和方差之比l l利用周期图法做谱估计时，周期图利用周期图法做谱估计时，周期图X X( (k
6、k) )为复数，为复数，因此有因此有1919l l考虑高斯变量情形，则XR(k)和XI(k)也是高斯变量；l l卡方分布公式：2020l l因此有l l因此变异系数为l l根据周期图公式，自由度n=2，因此变异系数为1，相对误差达到100%，估计极其不准确。2121l l2. 平均化处理l l对周期图进行平均化处理可以减小变异系数，对周期图进行平均化处理可以减小变异系数，得到较高的谱估计精度。得到较高的谱估计精度。l l将序列将序列x(n)x(n)分段，求各段周期图，再做平均，分段，求各段周期图，再做平均，此时有：此时有：l l各段频率分量的实部与虚部互相独立各段频率分量的实部与虚部互相独立2
7、222l l此时卡方分布自由度为n=2q. l l变异系数为 ，表明随着平均谱 分段数q增大，变异系数减小，从而可以通 过增加分段数来减小变异性。l l对连续随机过程，样本总体长度为T(T=Nt), t为采样间隔，分段长度Te(T/q)，那么分析带宽Be=1/Te，因此有 23234. Matlab应用举例l l求卷积 x = randn(1,100);w = 10;y = conv(ones(1,w)/w,x);avgs = y(10:99);plot(avgs)24242525l lEnsemble averagew = 10;for i = 1:w;X(i,:) = randn(1,10
8、0);endAVGS = mean(X);plot(AVGS)26262727l l交叉相关xcorr,交叉协方差xcorrl l两组信号的交叉相关等价于两组信号的卷积(其中一组逆序fliplr)l lxcov的作用是在计算交叉相关之前去掉输入的均值；l l信号与自身的交叉相关称为自相关。2828x = -1 0 1;y = 0 1 2;xcorr(x,y)conv(x,fliplr(y)xcov(x,y)xcov(x,x)xcov(x,y-1)2929例1.交叉相关l l考虑简单的目标测距系统，输出信号x与返回信号y之间的关系为：l l其中alpha为衰减因子，d为时滞，beta为信道噪声，
9、若T为信号返回时间，则x，y将在n=T时相关。目标则可以定位在vT的距离，其中v为信号的传播速率。3030x = 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)31313232l l进行进行DFTDFTdelta1 = 1 zeros(1,11);delta1 = 1 zeros(1,11);fftgui(delta1)fftgui(
10、delta1)delta2 = 0 1 zeros(1,10);delta2 = 0 1 zeros(1,10);fftgui(delta2)fftgui(delta2)deltaNyq = zeros(1,6),1,zeros(1,5);deltaNyq = zeros(1,6),1,zeros(1,5);fftgui(deltaNyq)fftgui(deltaNyq)square = zeros(1,4),ones(1,4),zeros(1,4);square = zeros(1,4),ones(1,4),zeros(1,4);fftgui(square)fftgui(square)t =
11、 linspace(0,1,50);t = linspace(0,1,50);periodic = sin(2*pi*t);periodic = sin(2*pi*t);fftgui(periodic)fftgui(periodic)33333434FFT Demol l太阳黑子sigdemo1playshow fftdemophoneplayshow sunspots3535功率谱分析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)/n
12、fft;Power = Y.*conj(Y)/nfft;plot(f,Power)title(Periodogram)36363737figureryy = xcorr(y,y);Ryy = fft(ryy,512);plot(f, abs(Ryy)title(DFT of Autocorrelation)38383939t = 0:1/100:10-1/100;x = sin(2*pi*15*t) + sin(2*pi*30*t);periodogram(x,512,100);figurepwelch(x,512,100);figurepmtm(x,512,100);404041414242l lSPTool4343
