自己编写算法的功率谱密度的三种matlab实现方式.docx
《自己编写算法的功率谱密度的三种matlab实现方式.docx》由会员分享，可在线阅读，更多相关《自己编写算法的功率谱密度的三种matlab实现方式.docx（14页珍藏版）》请在冰豆网上搜索。
自己编写算法的功率谱密度的三种matlab实现方式
功率谱密度的三种matlab实现方式
一：
实验目的：
（1）把握三种算法的概念、应用及特点；
（2）了解谱估量在信号分析中的作用；
（3）能够利用burg法对信号作谱估量，对信号的特点加以分析。
二;实验内容：
（1）简单说明三种方式的原理。
（2）用三种方式编写程序，在matlab中实现。
（3）将计算结果表示成图形的形式，给出三种情形的功率谱图。
（4）比较三种方式的特性。
（5）写出自己的心得体会。
三：
实验原理：
1.周期图法：
周期图法又称直接法。
它是从随机信号x（n）中截取N长的一段，把它视为能量有限x（n）真实功率谱
的估量
的抽样.
以为随机序列是广义平稳且各态遍历的，能够用其一个样本x（n）中的一段
来估量该随机序列的功率谱。
这固然必然带来误差。
由于对
采纳DFT，就默许
在时域是周期的，和
在频域是周期的。
这种方式把随机序列样本x（n）看成是截得一段
的周期延拓，这也确实是周期图法那个名字的来历。
2.相关法（间接法）:
这种方式以相关函数为媒介来计算功率谱，因此又叫间接法。
这种方式的具体步骤是：
第一步：
从无穷长随机序列x（n）中截取长度N的有限长序列列
第二步：
由N长序列
求（2M-1）点的自相关函数
序列。
（2-1）那个地址，m=-（M-1）…,-1,0,1…,M-1，M
N，
是双边序列，可是由自相关函数的偶对称性式，只要求出m=0,。
。
。
M-1的傅里叶变换，另一半也就明白了。
第三步：
由相关函数的傅式变换求功率谱。
即
以上进程中经历了两次截断，一次是将x（n）截成N长，称为加数据窗，一次是将x（n）截成（2M-1）长，称为加延迟窗。
因此所得的功率谱仅是近似值，也叫谱估量，式中的
代表估值。
一样取M<因此，在FFT问世之前，相关法是最经常使用的谱估量方式。
三：
Burg法：
AR模型功率谱估量又称为自回归模型，它是一个全极点的模型，要利用AR模型进行功率谱估量须通过levinson_dubin递推算法由
Yule-Walker方程求得AR的参数:
σ2，α1α2…αp。
计算中，预测系数必需知足Lenvinson-Durbin递推关系，而且可直接计算而无需第一计算自相关系数。
这种方式的优势确实是对未知数据不需要做任何假设，估量精度较高。
其缺点是在分析噪声中的正弦信号时，会引发谱线割裂，且谱峰的位置和正弦信号的相位有专门大的关系。
Burg算法是使前向预测误差和后向预测误差均方误差之和最小来求取Km的，它不对已知数据段之外的数据做以为假设。
计算m阶预测误差的递推表示公式如下：
求取反射系数的公式如下：
关于平稳随机进程，能够历时刻平均代替集合平均，因此上式可写成：
如此即可求得AR模型的反射系数。
将m阶AR模型的反射系数和m-1阶AR模型的系数代入到Levinson关系式中，能够求得AR模型其他的p-1个参数。
Levinson关系式如下：
m阶AR模型的第m+1个参数G，
其中
是预测误差功率，可由递推公式
求得。
易知为进行该式的递推，必需明白0阶AR模型误差功率
，
可知该式由给定序列易于求得。
完成上述进程，即最终求得了表征该随机信号的AR模型的p+1个参数。
然后依照
即可求得该随机信号的功率谱密度。
四．实验内容：
实验程序及实验图像
周期法：
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）;
gridon
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'）;
gridon;
nfft=256
nfft=1024
Burg法：
clear
Fs=1000%设置关键变量，可通过调剂这些变量观看不同成效
f1=0.2;
f2=0.213;
nfft=128;%取样点数
p=50;%阶数p应该选择在N/3delta=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;
forl=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）;
fori=1:
length（f）%循环递推
sum=0;
fork=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过小时谱的分辨率又不行。
对其改良的要紧方式有二种，即平均和滑腻，平均确实是将截取的数据段
再分成L个小段，别离计算功率谱后取功率谱的平均，这种方式使估量的方差减少，但误差加大，分辨率下降。
滑腻是用一个适当的窗函数
与算出的功率谱
进行卷积，使谱线滑腻。
这种方式得出的谱估量是无偏的，方差也小，但分辨率下降。
（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算法的设计原理，自己收成专门大，这也应了有付出就有回报的那句话。
升级会员 