自己编写算法的功率谱密度的三种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)序列。

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

Nn?0?? (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

nfft=1024

nfft=10000

相关法：

clear;

Fs=1000; %采样频率

n=0:Fs;%产生含有噪声的序列 nfft=1024;

xn=sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n)); cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); %求出功率谱密度 Pxx=abs(CXk);

index=0:round(nfft/2-1); f=index/nfft;

plot_Pxx=10*log10(Pxx(index+1)); plot(f,plot_Pxx); xlabel('频率'); ylabel('功率/DB'); grid on;

nfft=256

nfft=1024

Burg法：

clear

Fs=1000 %设置关键变量，可通过调节这些变量观察不同效果 f1=0.2; f2=0.213;

nfft=128; %取样点数

p=50; %阶数p应该选择在N/3

f=0:1/1000:0.5; n=1:Fs;

xn=sin(2*pi*f1*n)+sqrt(2)*sin(2*pi*f2*n)+randn(size(n)); figure;

plot(n,xn); title('burg时域'); xn= xn(:); N=length(xn); ef = xn; eb = xn; a = 1; for l=1:p

efp = ef(2:end);%m-1阶前向预测误差 ebp = eb(1:end-1);%m-1阶后向预测误差

num = -2.*ebp'*efp;%1km分子多项式

den = efp'*efp+ebp'*ebp;%1km的分母多项式

k(l) = num ./ den;%计算反射系数

% 更新前向和后向预测误差

ef = efp + k(l)*ebp;%各阶前向预测误差 eb = ebp + k(l)*efp;%各阶后向预测误差

% 计算模型参数

a=[a;0] + k(l)*[0;conj(flipud(a))];%AR模型参数a end

a1=a(2:p+1);

for i=1:length(f) %循环递推 sum=0; for k=1:p

sum=sum+a1(k)*exp(-m*2*pi*f(i)*k); end

Pbrg(i)=delta/(abs(1+sum))^2;

Pbrg_f(i)=10*log10(Pbrg(i));%求出功率谱 end figure

plot(f,Pbrg_f); title('burg频域');

nfft=128

nfft=256

五：结果比较分析

(1) 在采样点相同的时，周期图法的特点是离散性大，曲线粗糙，方差较大，但是分辨率较高；采用周期突发估计得出的功率谱很不平滑，相应的估计协方差比较大。而且采用增加采样点的办法也不能吃周期图变得更加平滑，这是周期图法的缺点。周期图法得出的估计谱方差特性不好：当数据长度N太大时，扑线的起伏加剧；N太小时谱的分辨率又不好。对其改进的主要方法有二种，即平均和平滑，平均就是将截取的数据段xN(n)再分成L个小段，分别计算功率谱后取功率谱的平均，这种方法使估计的方差减少，但偏差加大，分辨率下降。平滑是用一个适当的窗函数W(e)与算出的功率谱SX(e)进行卷

jwjw?积，使谱线平滑。这种方法得出的谱估计是无偏的，方差也小，但分辨率下降。

(2)相关法当延迟与数据长度之比很小时，可以有良好的估计精度，相关法的收敛性较好，曲线平滑，方差较小，但是功率谱主瓣较宽，分辨率低。

（3）用Burg算法进行功率谱估计时令前后向预测误差功率之和最小,即对前向序列误差和后向序列误差前后都不加窗。Burg算法是建立在数据基础之上的，避免了先计算自相关函数从而提高计算速度。是较为通用的方法，计算不太复杂 并且分辨率优于自相关法。但对于白噪声加正弦信号有时会出现谱线分裂现象，并且从上两个图中可以看出burg法产生的功率谱曲线比较平滑即方差小，分辨率高，可以明显的观察到两个谱峰，在降低模型阶次后谱的分辨率降低(两个谱峰几乎变成一个谱峰)，但是曲线的平滑性更好。

并且采样点数越大，谱图的分辨率就越高。对比nfft=128和nfft=256即可发现。除此之外还发现对于上面三种情况采样点数越大，其功率谱密度也越大。还有就是阶数p应该选择在N/3

第一次作业在上课的时候被老师点到，当时老师问到burg法产生的图像是否正确，当时我觉得应该没错误啊，不过因为自己对整个过程的理解有限，没听出老师要表达什么意思，所以就只能沉默了。。不过这篇实验报告确实不是复制粘贴得到的，是经过查询很多资料自己写

出来的。

后来下课后自己看了下程序，发现老师说的应该是burg法产生的图像没有两个信号的谱峰，这是因为X轴没有进行单位的调节导致谱图在最左侧位置所以看不到。并且自己上次是直接调用matlab中的Pburg函数，这虽然简单省事但是对burg算法无法进行较为深入的认识，并且后来问了同学才知道老师是不允许直接调用burg函数的。虽然上课时被老师点到，出了点小丑，不过我觉得这倒是对自己的一个激励，否则的话我到以为自己做的还行，实际上自己做的根本没有达到老师的要求，在下课后自己为了用burg算法进行求解功率谱自己又在网上查找了很多资料并结合课本上的介绍，网上很多资料要么是直接调用burg函数的方法，要么就是描述有些问题，所以花费了不少时间，我看了很多次，对有些不懂得语句和公式来回看课本和程序，最后结合课本和网上的介绍，自己搞懂了burg算法的设计过程，现在看来虽然这个过程比着直接调用matlab中的pburg函数耗时很多，但是在这个过程中我不仅对matlab的编程有了更深入的了解，而且也明白了burg算法的设计原理，自己收获很大，这也应了有付出就有回报的那句话。

本文来源：https://www.bwwdw.com/article/8hs3.html