《2022年自己编写算法的功率谱密度的三种matlab实现方法》由会员分享，可在线阅读，更多相关《2022年自己编写算法的功率谱密度的三种matlab实现方法（13页珍藏版）》请在金锄头文库上搜索。
1、功率谱密度的三种matlab 实现方法一：实验目的：(1)掌握三种算法的概念、应用及特点；(2)了解谱估计在信号分析中的作用；(3) 能够利用 burg 法对信号作谱估计，对信号的特点加以分析。二; 实验内容：（1）简单说明三种方法的原理。（2）用三种方法编写程序，在matlab 中实现。（3）将计算结果表示成图形的形式，给出三种情况的功率谱图。（4）比较三种方法的特性。（5）写出自己的心得体会。三：实验原理 ：1. 周期图法：周期图法又称直接法。它是从随机信号x(n)中截取 N 长的一段，把它视为能量有限x(n)真实功率谱)(jwxeS的估计)(jwxeS的抽样 . 认为随机序列是广义平稳且
2、各态遍历的，可以用其一个样本x(n)中的一段)(nxN来估计该随机序列的功率谱。这当然必然带来误差。由于对)(nxN采用 DFT，就默认)(nxN在时域是周期的，以及)(kxN在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(nxN的周期延拓，这也就是周期图法这个名字的来历。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 13 页2.相关法（间接法 ): 这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。 这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列)(nxN第二步：由 N 长序列
3、)(nxN求 （2M-1） 点的自相关函数)(mRx序列。)()(1)(10mnxnxNmRNnNNx(2-1) 这里，m=-(M-1),-1,0,1,M-1，MN，)(mRx是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。 。 。,M-1 的傅里叶变换，另一半也就知道了。第 三 步 ： 由 相 关 函 数 的 傅 式 变 换 求 功 率 谱 。 即jwmMMmXjwxemReS)()(1)1(以上过程中经历了两次截断， 一次是将 x(n)截成 N 长，称为加数据窗，一次是将 x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jwxeS代表
4、估值。一般取 MN，因为只有当M 较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。因此，在FFT问世之前，相关法是最常用的谱估计方法。三：Burg 法：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 13 页AR 模型功率谱估计又称为自回归模型，它是一个全极点的模型，要利用 AR模型进行功率谱估计须通过levinson_dubin递推算法由Yule-Walker 方程求得 AR的参数: 2，12p。计算中，预测系数必须满足Lenvinson-Durbin递推关系，并且可直接计算而无需首先计算自
5、相关系数。这种方法的优点就是对未知数据不需要做任何假设， 估计精度较高。 其缺点是在分析噪声中的正弦信号时，会引起谱线分裂， 且谱峰的位置和正弦信号的相位有很大的关系。Burg 算法是使前向预测误差和后向预测误差均方误差之和最小来求取 Km的，它不对已知数据段之外的数据做认为假设。计算m阶预测误差的递推表示公式如下：x(n)(n)(n)(n)1)-(n(n)1)-(n(n)(n)0f0f1-mm1-bm1-mf1-mmeeeeekeekebbmbmf求取反射系数的公式如下：1)-(n(n)1)-(n(n)2-2b1 -m2f1 -mb1-mf1-mmeeeeEEk对于平稳随机过程，可以用时间平
6、均代替集合平均， 因此上式可写成：p,2, 1,1)-(n(n)1)-(n(n)2-1 -21 -21-1-mn1 -1-，mNmnbmfmNbmfmmeeeek这样便可求得 AR模型的反射系数。将 m 阶 AR 模型的反射系数和m-1 阶 AR 模型的系数代入到精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 13 页Levinson 关系式中，可以求得AR模型其他的 p-1 个参数。Levinson 关系式如下：1-m1,2,ii),-(m(i)(i)1-m1-mm，akaamm阶 AR模型的第 m+1个参数 G ，2mG其中m是预测
7、误差功率， 可由递推公式m21(1)mmK求得。易知为进行该式的递推，必须知道0 阶 AR模型误差功率0，20( )(0)xE XnR可知该式由给定序列易于求得。完成上述过程， 即最终求得了表征该随机信号的 AR模型的 p+1 个参数 。然后根据22()()jwjwxwSeH e即可求得该随机信号的功率谱密度。四实验内容：实验程序及实验图像周期法：Fs=1000;nfft=10000; %2nn=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);
8、 % 求解 PSDt=0:round(nfft/2-1);f=t/nfft;P=10*log10(Pxx(t+1); % 纵坐标的单位为dBplot(f,P);grid on精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 13 页nfft=200 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 13 页 nfft=1024 nfft=10000 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 13 页相关法：clear;F
9、s=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 ;
10、nfft=256 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 13 页nfft=1024 Burg 法：clearFs=1000 % 设置关键变量，可通过调节这些变量观察不同效果f1=0.2;f2=0.213;nfft=128; % 取样点数p=50; % 阶数 p应该选择在 N/3pN/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(bur
11、g时域 );xn= xn(:);N=length(xn);ef = xn;eb = xn;a = 1;for l=1:p精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 13 页 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;% 各阶
12、前向预测误差 eb = ebp + k(l)*efp;% 各阶后向预测误差% 计算模型参数 a=a;0 + k(l)*0;conj(flipud(a);%AR 模型参数 aenda1=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);% 求出功率谱endfigureplot(f,Pbrg_f);title(burg频域 ); 精选学习资料 - - - - - -
13、 - - - 名师归纳总结 - - - - - - -第 9 页，共 13 页nfft=128 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 13 页 nfft=256 五：结果比较分析(1) 在采样点相同的时，周期图法的特点是离散性大，曲线粗糙，方差较大，但是分辨率较高； 采用周期突发估计得出的功率谱很不平滑，相应的估计协方差比较大。 而且采用增加采样点的办法也不能吃周期图变得更加平滑， 这是周期图法的缺点。 周期图法得出的估计谱方差特性不好：当数据长度N 太大时，扑线的起伏加剧；N 太小时谱的分辨率又不好。对其改进的主要方法有二
14、种，即平均和平滑，平均就是将截取的数据段)(nxN再分成 L 个小段，分别计算功率谱后取功率谱的平均，这种方法使估计的方差减少，但偏差加大，分辨率下降。平滑是用一个适当的窗函数)(jweW与算出的功率谱)(jweSX进行卷精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 13 页积，使谱线平滑。这种方法得出的谱估计是无偏的，方差也小，但分辨率下降。(2)相关法当延迟与数据长度之比很小时，可以有良好的估计精度，相关法的收敛性较好，曲线平滑，方差较小，但是功率谱主瓣较宽，分辨率低。（3） 用 Burg 算法进行功率谱估计时令前后向预测误差功
15、率之和最小,即对前向序列误差和后向序列误差前后都不加窗。Burg 算法是建立在数据基础之上的， 避免了先计算自相关函数从而提高计算速度。是 较 为 通 用 的方法，计算不太复杂 并且分 辨 率 优 于 自 相 关 法 。但对于白噪声加正弦信号有时会出现谱线分裂现象， 并且从上两个图中可以看出 burg 法产生的功率谱曲线比较平滑即方差小，分辨率高，可以明显的观察到两个谱峰，在降低模型阶次后谱的分辨率降低(两个谱峰几乎变成一个谱峰 )，但是曲线的平滑性更好。并且采样点数越大，谱图的分辨率就越高。对比 nfft=128和 nfft=256即可发现。除此之外还发现对于上面三种情况采样点数越大， 其功
16、率谱密度也越大。还有就是阶数p 应该选择在 N/3pN/2比较合适，这个通过测试可以得到验证。六：心得体会第一次作业在上课的时候被老师点到，当时老师问到burg 法产生的图像是否正确，当时我觉得应该没错误啊， 不过因为自己对整个过程的理解有限，没听出老师要表达什么意思，所以就只能沉默了。 。不过这篇实验报告确实不是复制粘贴得到的， 是经过查询很多资料自己写精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 13 页出来的。后来下课后自己看了下程序，发现老师说的应该是burg 法产生的图像没有两个信号的谱峰，这是因为X 轴没有进行单位的调节
17、导致谱图在最左侧位置所以看不到。并且自己上次是直接调用matlab中的Pburg函数， 这虽然简单省事但是对burg算法无法进行较为深入的认识， 并且后来问了同学才知道老师是不允许直接调用burg 函数的。虽然上课时被老师点到， 出了点小丑， 不过我觉得这倒是对自己的一个激励，否则的话我到以为自己做的还行，实际上自己做的根本没有达到老师的要求， 在下课后自己为了用burg算法进行求解功率谱自己又在网上查找了很多资料并结合课本上的介绍，网上很多资料要么是直接调用 burg函数的方法，要么就是描述有些问题，所以花费了不少时间，我看了很多次， 对有些不懂得语句和公式来回看课本和程序， 最后结合课本和网上的介绍， 自己搞懂了 burg算法的设计过程，现在看来虽然这个过程比着直接调用matlab中的 pburg函数耗时很多，但是在这个过程中我不仅对matlab的编程有了更深入的了解，而且也明白了 burg算法的设计原理，自己收获很大，这也应了有付出就有回报的那句话。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 13 页
