自功率谱估计性能分析及Matlab仿真Word文档格式.docx
《自功率谱估计性能分析及Matlab仿真Word文档格式.docx》由会员分享,可在线阅读,更多相关《自功率谱估计性能分析及Matlab仿真Word文档格式.docx(18页珍藏版)》请在冰豆网上搜索。
周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。
取平稳随机信号的有限个观察值,求出其傅里叶变换
然后进行谱估计
周期图法应用比较广泛,主要是由于它与序列的频谱有直接的对应关系,并且可以采用FFT快速算法来计算。
但是,这种方法需要对无限长的平稳随机序列进行截断,相当于对其加矩形窗,使之成为有限长数据。
同时,这也意味着对自相关函数加三角窗,使功率谱与窗函数卷积,从而产生频谱泄露,容易使弱信号的主瓣被强信号的旁瓣所淹没,造成频谱的模糊和失真,使得谱分辨率较低[1]。
该方法基于Matlab实现的程序:
clearall;
load testx;
N=4096;
Fn=-0.5:
1/N:
0.5-1/N;
px=fft(x,N);
pmax=max(px);
%归一化
px=px/pmax;
px=10*log10(px+0.000001);
plot(Fn,fftshift(px));
gridon;
图1周期图法
图2周期图法
说明:
(1)本报告仿真中所采用的用于功率谱估计的数据文件来自参考文献[3]的test.dat。
该数据为128点复序列(图3),由复数噪声加上四个复正弦组成。
其归一化频率分别是:
。
图3复序列
(2)从仿真图可以清晰看到,和不能完全分开,仅在波形的顶部能看出是两个频率分量;
此外,当数据长度太大时(图1),谱曲线呈现较大的起伏;
当数据长度太小时(图2),谱的分辨率又不好。
据此,周期图法不满足一致性估计条件。
2.2自相关法(BT法)
自相关法的理论基础是维纳—辛钦定理。
1958年Blackman和Tukey给出了这一方法的具体实现。
对于平稳随机信号来说,其自相关函数是确定性函数,故其功率谱也是确定的。
这样可由平稳随机离散信号的有限个离散值求出自相关函数
然后在内对做傅里叶变换,得到功率谱
该方法基于Matlab实现的程序:
clear all;
loadtestx;
N=4096;
Fn=-0.5:
1/N:
0.5-1/N;
Mlag=64;
rx=xcorr(x,Mlag,'
unbiased'
);
px=fft(rx,N);
pmax=max(px);
px=px/pmax;
px=10*log10(px+0.000001);
plot(Fn,fftshift(px));
gridon;
图4自相关法不加窗
图5自相关法不加窗
图6自相关法使用汉明窗( Hamming )
(1)该方法先由序列估计出自相关函数,然后对进行傅里叶变换,便得到的功率谱估计。
当延迟与数据长度之比很小时,可以有良好的估计精度。
(2)图4是用自相关法(BT法)求出的功率谱,没有加窗;
图5也是用自相关法(BT法)求出的功率谱,,没有加窗;
图6同样是采用自相关法求出的功率谱,,使用了汉明窗。
显然,自相关函数的延迟越小,谱变得越平滑。
2.3Welch法
该方法的基本原理是在对随机序列分段时,使每一段有部分重叠,然后对每一段数据用一个合适的窗函数进行平滑处理,最后对各段谱求平均。
这样可得功率谱
其中
这里为窗函数。
该方法基于Matlab实现的程序:
clearall;
loadtest x;
N=4096;
Fn=-0.5:
0.5-1/N;
xpsd=pwelch(x,hamming(33),16,N,'
whole'
);
mmax=max(xpsd);
xpsd=xpsd/mmax;
xpsd=10*log10(xpsd+0.000001);
plot(Fn,fftshift(xpsd));
gridon;
图7Welch法不叠合使用汉明窗(Hamming)
图8Welch法叠合16点使用汉明窗(Hamming)
图9Welch法叠合16点使用矩形窗( Boxcar)
图10Welch法叠合16点使用布莱克曼窗(Blackman )
说明:
(1)因为Welch法允许各段数据交叠,所以数据段数会增加,使方差得到更大的改善,但是数据的交叠减小了每一段数据的不相关性,使方差的减小不会达到理论程度。
另外,采用合适的窗函数可以减少信号的频谱泄露,同时也可以增加谱峰的宽度,从而提高分辨率。
(2) 图7是利用Welch法求出的周期图,共分四段,每段32点,没有叠合,使用了汉明窗;
图8也是利用Welch法求出的周期图,共分四段,每段32点,,使用了汉明窗;
图9是利用Welch法求出的周期图,共分四段,每段32点,交叠数为16,且使用了矩形窗;
图10是利用Welch法求出的周期图,共分四段,每段32点,交叠数为16,使用了布莱克曼窗。
从图中可以看出,由矩形窗处理的谱估计的主瓣宽度最窄,分辨率最好,但是其起伏性较大,所以其方差特性最差。
由汉明窗和布莱克曼窗得到的谱估计的主瓣宽度最宽,因此其分辨率相对较差,但其旁瓣较小,大大改善了由矩形窗处理的谱估计旁瓣较大所产生的谱失真。
因此,选择不同的窗函数其主瓣宽度不一样,造成谱估计的分辨率也不相同。
2.4经典功率谱估计的性能比较
由以上的Matlab仿真图形和相关结果分析,我们得到了经典谱估计算法性能的直观比较:
(1)周期图法得到的功率谱分辨率最高,但是方差性能最差,功率谱起伏剧烈,容易出现虚假谱峰。
(2)自相关法(BT法)由于使用了平滑窗对周期图法估计的功率谱进行了平滑,因此方差性能较好,功率谱比周期图法估计的要平滑,但其分辨率比周期图法低。
(3) Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的,估计的功率谱也最为平滑,但这是以分辨率的下降及偏差的增大为代价的。
综合上述讨论,我们可以对经典谱估计的算法作大致的总结[3]:
(1)功率谱估计,不论是直接法还是间接法都可以用FFT快速计算,且物理概念明确,因而仍是目前较常用的谱估计方法。
(2) 谱的分辨率较低,它正比于,是所使用的数据长度。
(3)方差性能不好,不是真实功率谱的一致估计,且增大时,功率谱起伏加剧。
(4) 周期图的平滑和平均是和窗函数的使用紧密关联的,平滑和平均主要是用来改善周期图的方差性能,但往往又减小了分辨率且增加了偏差,没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善,因此使用窗函数只是改进估计质量的一个技巧问题,并不能从根本上解决问题。
3 现代功率谱估计
由前一章的讨论我们可知,经典功率谱估计方法的方差性能较差,分辨率较低。
而现代谱估计技术的目标都是旨在努力改善谱估计的分辨率。
参数模型法是现代谱估计的主要内容,参数模型主要分为AR模型、MA模型和ARMA模型。
由于AR模型具有一系列好的性能,因此是被研究最多并获得广泛应用的一种模型。
本报告中现代功率谱估计的仿真基于的是AR模型。
3.1自相关法
假定观察到得数据为,而对于无法观察到得区间(即),的样本假定为0,观测数据区间之外的数据为0,在均方误差意义下使得数据的预测误差功率最小。
由于自相关矩阵是Toeplitz矩阵,而且又为正定的,故可利用Levinson-Durbin递归算法高效求解,得到AR模型参数。
该方法基于Matlab实现的程序:
clearall;
loadtestx;
N=4096;
fn=-0.5:
1/N:
0.5-1/N;
xpsd=pyulear(x,20,N);
pmax=max(xpsd);
xpsd=xpsd/pmax;
xpsd=10*log10(xpsd+0.000001);
plot(fn,fftshift(xpsd));
gridon;
图11 自相关法
图12自相关法
图13 自相关法
(1)图11、12和13是用自相关法求出的AR谱曲线,阶次分别等于10,20和30。
可以看出,在阶次较低时(图11),分辨率和检测能力均不好。
当时,和处的两个正弦刚刚可以分开,在和处的两个正弦也可以检出。
因此必须通过提高阶次来达到分辨出间隔较小的频率点的效果。
(2) AR模型的自相关法等效于对前向预测的误差序列前后加窗,加窗的结果是使得自相关法的分辨率降低。
数据越短,分辨率越不好。
3.2协方差法
协方差法与自相关法的区别主要在于预测误差功率求和式的上下限取得不同。
由于协方差法对于观察区间外样本并未假定为0,故预测误差功率表达式中的总是落在观察区间中,为此预测误差功率的求和上下限必须在之间。
但由此得到的自相关矩阵是对此的半正定矩阵,且不具有Toeplitz性质,故不能采用Levinson-Durbin递归算法求解,因此得到的AR模型可能不稳定。
该方法基于Matlab实现的程序:
clearall;
loadtestx;
N=4096;
tn=-0.5:
0.5-1/N;
xpsd=pcov(x,10,N);
pmax=max(xpsd);
xpsd=xpsd/pmax;
xpsd=10*log10(xpsd+0.000001);
plot(tn,fftshift(xpsd));
gridon;
图14协方差法
(1) 可以看到谱图在信号源频率处:
谱线狭窄突出,其他处谱线起伏较为剧烈。
(2)采用协方差法对信号进行建模,能够较好地反映出信号真正的模型。
3.3修正的协方差法
AR谱估计的协方差算法基于的是最小化前向预测误差。
而修正的协方差算法基于的是最小化前向和后向预测误差。
这样使得它的误差功率的计算是在相对于协方差法多一倍的数据点上进行,这在观察数据长度很短的情况下,是非常有利的,但这要求信号在正反两个方向上呈现相同的特性。
此外,由此得到的自相关矩阵不具有Toeplitz性质,故其正则方程不能采用Levinson-Durbin递归算法求解。
该方法基于Matlab实现的程序:
clearall;
loadtestx;
N=4096;
fn=-0.5:
1/N:
0.5-1/N;
xpsd=pmcov(x,10,N);
pmax=max(xpsd);
xpsd=xpsd/pmax;
xpsd=10*log10(xpsd+0.000001);
plot(fn,fftshift(xpsd));
grid on;
图15修正的协方差法