好文档 - 专业文书写作范文服务资料分享网站

(完整word版)自己编写算法的功率谱密度的三种matlab实现方法

天下 分享 时间: 加入收藏 我要投稿 点赞

功率谱密度的三种matlab实现方法

一:实验目的:

(1)掌握三种算法的概念、应用及特点; (2)了解谱估计在信号分析中的作用;

(3) 能够利用burg法对信号作谱估计,对信号的特点加以分析。 二;实验内容:

(1)简单说明三种方法的原理。

(2)用三种方法编写程序,在matlab中实现。

(3)将计算结果表示成图形的形式,给出三种情况的功率谱图。 (4)比较三种方法的特性。 (5)写出自己的心得体会。 三:实验原理: 1.周期图法:

周期图法又称直接法。它是从随机信号x(n)中截取N长的一段,把它视为能量有限x(n)真实功率谱Sx(ejw)的估计Sx(ejw)的抽样.

认为随机序列是广义平稳且各态遍历的,可以用其一个样本x(n)中的一段xN(n)来估计该随机序列的功率谱。这当然必然带来误差。由于对xN(n)采用DFT,就默认xN(n)在时域是周期的,以及xN(k)在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段xN(n)的周期延拓,这也就是周期图法这个名字的来历。

2.相关法(间接法):

这种方法以相关函数为媒介来计算功率谱,所以又叫间接法。这种方法的具体步骤是:

第一步:从无限长随机序列x(n)中截取长度N的有限长序列列

xN(n)

第二步:由N长序列xN(n)求(2M-1)点的自相关函数Rx(m)序列。

1Rx(m)?N???xn?0N?1N(n)xN(n?m)

(2-1) 这里,m=-(M-1)…,-1,0,1…,M-1,MN,Rx(m)是双边序列,但是由自相关函数的偶对称性式,只要求出m=0,。。。,M-1的傅里叶变换,另一半也就知道了。

第三步:由相关函数的傅式变换求功率谱。即

Sx(e)??jwM?1?Xm??(M?1)?R(m)e?jwm

以上过程中经历了两次截断,一次是将x(n)截成N长,称为加数据窗,一次是将x(n)截成(2M-1)长,称为加延迟窗。因此所得的功率谱仅是近似值,也叫谱估计,式中的Sx(ejw)代表估值。一般取M<

AR模型功率谱估计又称为自回归模型,它是一个全极点的模型,要利用AR模型进行功率谱估计须通过levinson_dubin递推算法由 Yule-Walker方程求得AR的参数:σ

2,

α1α2…αp。

计算中,预测系数必须满足Lenvinson-Durbin递推关系,并且可直接计算而无需首先计算自相关系数。这种方法的优点就是对未知数据不需要做任何假设,估计精度较高。其缺点是在分析噪声中的正弦信号时,会引起谱线分裂,且谱峰的位置和正弦信号的相位有很大的关系。

Burg算法是使前向预测误差和后向预测误差均方误差之和最小来求取Km的,它不对已知数据段之外的数据做认为假设。计算m阶预测误差的递推表示公式如下:

em(n)?em-1(n)?kmem-1(n-1)em(n)?em-1(n-1)?kmem-1(n)e0(n)?e0(n)?x(n)

求取反射系数的公式如下:

-k?E{[?[-1)]}e(n)]e(n

mf2b2m-1m-1ffbbbffb

2E[(n)(n-1)]eem-1m-1fb对于平稳随机过程,可以用时间平均代替集合平均,因此上式可写成:

???-,m?1,2,?,pk????????-1)?e(n)e(nmn?mN-1f2(n)(n-1)?eem-1m-12bm-1m-1N-1?fb2n?m

这样便可求得AR模型的反射系数。

将m阶AR模型的反射系数和m-1阶AR模型的系数代入到

Levinson关系式中,可以求得AR模型其他的p-1个参数。 Levinson关系式如下:

(i)?(i)?(m-i),i?1,2,?,m-1aaka

mm-1mm-1m阶AR模型的第m+1个参数G,G2?ρm其中ρm是预测误差功率,可由

2)求得。 递推公式ρm?ρm?1(1?Km易知为进行该式的递推,必须知道0阶AR模型误差功率ρ0,

?2?ρ0?E?X(n)??Rx(0)

可知该式由给定序列易于求得。完成上述过程,即最终求得了表征该随机信号的AR模型的p+1个参数 。然后根据

2Sx(ejw)?σwH(ejw)

2即可求得该随机信号的功率谱密度。

四.实验内容: 实验程序及实验图像 周期法:

Fs=1000;

nfft=10000; %2^n n=0:Fs;

x=sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n)); X=fft(x,nfft);

Pxx=abs(X).^2/length(n); %求解PSD t=0:round(nfft/2-1); f=t/nfft;

P=10*log10(Pxx(t+1)); %纵坐标的单位为dB plot(f,P); grid on

nfft=200

44t9x7umfg79ew80o94h77xpo5846y00r0s
领取福利

微信扫码领取福利

微信扫码分享