功率谱密度的三种matlab 实现方法一：实验目的：(1)掌握三种算法的概念、应用及特点； (2)了解谱估计在信号分析中的作用；(3) 能够利用burg 法对信号作谱估计，对信号的特点加以分析。
二;实验内容：（1）简单说明三种方法的原理。
（2）用三种方法编写程序，在matlab 中实现。
（3）将计算结果表示成图形的形式，给出三种情况的功率谱图。
（4）比较三种方法的特性。
（5）写出自己的心得体会。
三：实验原理： 1.周期图法：周期图法又称直接法。
它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样.认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(n x N 来估计该随机序列的功率谱。
这当然必然带来误差。
由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。
这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。
2.相关法（间接法):这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。
这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列)(n x N第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x ∧序列。
)()(1)(1m n x n xNm R N n N Nx +=∑-=∧(2-1) 这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。
,M-1的傅里叶变换，另一半也就知道了。
第三步：由相关函数的傅式变换求功率谱。
即jwm M M m Xjwx e m Re S ----=∧∧∑=)()(1)1(以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。
因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。
一般取M<<N ，因为只有当M 较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。
因此，在FFT 问世之前，相关法是最常用的谱估计方法。
三：Burg 法：AR 模型功率谱估计又称为自回归模型，它是一个全极点的模型，要利用AR 模型进行功率谱估计须通过levinson_dubin 递推算法由 Yule-Walker 方程求得AR 的参数:σ2，α1α2…αp 。
计算中，预测系数必须满足Lenvinson-Durbin 递推关系，并且可直接计算而无需首先计算自相关系数。
这种方法的优点就是对未知数据不需要做任何假设，估计精度较高。
其缺点是在分析噪声中的正弦信号时，会引起谱线分裂，且谱峰的位置和正弦信号的相位有很大的关系。
Burg 算法是使前向预测误差和后向预测误差均方误差之和最小来求取Km 的，它不对已知数据段之外的数据做认为假设。
计算m 阶预测误差的递推表示公式如下：x(n)(n)(n)(n)1)-(n (n)1)-(n (n)(n)0f 0f1-m m 1-b m 1-m f 1-m m e e e e ==+=+=e k e e k e bb m bm f求取反射系数的公式如下：}1)]-(n [(n)]{[1)]-(n (n)[2-2b1-m 2f1-m b1-m f 1-m me e e e +=E E k对于平稳随机过程，可以用时间平均代替集合平均，因此上式可写成：[][][][]{}p ,2,1,1)-(n(n)1)-(n (n)2-1-21-21-1-m n 1-1-，⋯=+=∑∑==m N mn bm fm N bm f m me e e e k这样便可求得AR 模型的反射系数。
将m 阶AR 模型的反射系数和m-1阶AR 模型的系数代入到Levinson 关系式中，可以求得AR 模型其他的p-1个参数。
Levinson 关系式如下：1-m 1,2,i i),-(m (i)(i)1-m 1-m m，⋯=+=a k a a mm 阶AR 模型的第m+1个参数G ，2m G ρ=其中m ρ是预测误差功率，可由递推公式m ρ=21(1)m mK ρ--求得。
易知为进行该式的递推，必须知道0阶AR 模型误差功率0ρ，20()(0)x E X n R ρ⎡⎤==⎣⎦可知该式由给定序列易于求得。
完成上述过程，即最终求得了表征该随机信号的AR 模型的p+1个参数 。
然后根据22()()jw jw x w S e H e σ=即可求得该随机信号的功率谱密度。
四．实验内容： 实验程序及实验图像 周期法：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 onnfft=200nfft=1024nfft=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=256nfft=1024Burg法：clearFs=1000 %设置关键变量，可通过调节这些变量观察不同效果f1=0.2;f2=0.213;nfft=128; %取样点数p=50; %阶数p应该选择在N/3<p<N/2delta=1;m=sqrt(-1);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:pefp = 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 enda1=a(2:p+1);for i=1:length(f) %循环递推sum=0;for k=1:psum=sum+a1(k)*exp(-m*2*pi*f(i)*k);endPbrg(i)=delta/(abs(1+sum))^2;Pbrg_f(i)=10*log10(Pbrg(i));%求出功率谱endfigureplot(f,Pbrg_f);title('burg频域');nfft=128nfft=256五：结果比较分析(1) 在采样点相同的时，周期图法的特点是离散性大，曲线粗糙，方差较大，但是分辨率较高；采用周期突发估计得出的功率谱很不平滑，相应的估计协方差比较大。
而且采用增加采样点的办法也不能吃周期图变得更加平滑，这是周期图法的缺点。
周期图法得出的估计谱方差特性不好：当数据长度N 太大时，扑线的起伏加剧；N 太小时谱的分辨率又不好。
对其改进的主要方法有二种，即平均和平滑，平均就是将截取的数据段)(n x N 再分成L 个小段，分别计算功率谱后取功率谱的平均，这种方法使估计的方差减少，但偏差加大，分辨率下降。
平滑是用一个适当的窗函数)(jw e W 与算出的功率谱)(jw e S X进行卷积，使谱线平滑。
这种方法得出的谱估计是无偏的，方差也小，但分辨率下降。
(2)相关法当延迟与数据长度之比很小时，可以有良好的估计精度，相关法的收敛性较好，曲线平滑，方差较小，但是功率谱主瓣较宽，分辨率低。
（3）用Burg算法进行功率谱估计时令前后向预测误差功率之和最小,即对前向序列误差和后向序列误差前后都不加窗。
Burg算法是建立在数据基础之上的，避免了先计算自相关函数从而提高计算速度。
是较为通用的方法，计算不太复杂并且分辨率优于自相关法。
但对于白噪声加正弦信号有时会出现谱线分裂现象，并且从上两个图中可以看出burg法产生的功率谱曲线比较平滑即方差小，分辨率高，可以明显的观察到两个谱峰，在降低模型阶次后谱的分辨率降低(两个谱峰几乎变成一个谱峰)，但是曲线的平滑性更好。
并且采样点数越大，谱图的分辨率就越高。
对比nfft=128和nfft=256即可发现。
除此之外还发现对于上面三种情况采样点数越大，其功率谱密度也越大。
还有就是阶数p应该选择在N/3<p<N/2比较合适，这个通过测试可以得到验证。
六：心得体会第一次作业在上课的时候被老师点到，当时老师问到burg法产生的图像是否正确，当时我觉得应该没错误啊，不过因为自己对整个过程的理解有限，没听出老师要表达什么意思，所以就只能沉默了。
不过这篇实验报告确实不是复制粘贴得到的，是经过查询很多资料自己写出来的。
后来下课后自己看了下程序，发现老师说的应该是burg法产生的图像没有两个信号的谱峰，这是因为X轴没有进行单位的调节导致谱图在最左侧位置所以看不到。
并且自己上次是直接调用matlab中的Pburg函数，这虽然简单省事但是对burg算法无法进行较为深入的认识，并且后来问了同学才知道老师是不允许直接调用burg函数的。
