功率谱密度估计方法的MATLAB实现Word文件下载.docx
《功率谱密度估计方法的MATLAB实现Word文件下载.docx》由会员分享，可在线阅读，更多相关《功率谱密度估计方法的MATLAB实现Word文件下载.docx（13页珍藏版）》请在冰豆网上搜索。
Y=awgn（xn,10）;
%加入信噪比为10db的高斯白噪声
subplot（2,1,1）;
plot（n,Y）
二、修正周期图法（加窗）谱估计程序
N=512;
%数据长度
M=32;
%汉明窗宽度
title（'
信号'
）
xlabel（'
时间'
）;
ylabel（'
幅度'
gridon;
window=hamming（M）;
%汉明窗
[Pxxf]=pwelch（Y,window,10,256,Fs）;
subplot（2,1,2）;
plot（f,10*log10（Pxx））;
title（['
修正周期图法谱估计N='
int2str（N）,'
M='
int2str（M）]）;
频率（Hz）'
功率谱密度'
2、仿真结果
三、最大熵法谱估计程序
fs=1;
%设采样频率
N=128;
%数据长度改变数据长度会导致分辨率的变化；
f1=0.2*fs;
%第一个sin信号的频率,f1/fs=0.2
f2=0.3*fs;
%第二个sin信号的频率，f2/fs=0.2或者0.3
P=10;
%滤波器阶数
n=1:
N;
s=sin（2*pi*f1*n/fs）+sin（2*pi*f2*n/fs）;
%s为原始信号
x=awgn（s,10）;
%x为观测信号，即对原始信号加入白噪声，信噪比10dB
figure
（1）;
%画出原始信号和观测信号
plot（s,'
b'
）,xlabel（'
）,ylabel（'
）,title（'
原始信号s'
grid;
plot（x,'
r'
观测信号x'
[Pxx1,f]=pmem（x,P,N,fs）;
%最大熵谱估计
figure
（2）;
plot（f,10*log10（Pxx1））;
频率（Hz）'
功率谱（dB）'
最大熵法谱估计模型阶数P='
int2str（P）,'
数据长度N='
int2str（N）]）;
四、Levinson递推法谱估计程序
%设采样频率为1
N=1000;
%第二个sin信号的频率，f1/fs=0.2或者0.3
M=16;
%滤波器阶数的最大取值,超过则认为代价太大而放弃
L=2*N;
%有限长序列进行离散傅里叶变换前，序列补零的长度
%s为原始信号
%x为观测信号，即对原始信号加入白噪声，信噪比10dB
）,axis（[0100-33]）,xlabel（'
%计算自相关函数
rxx=xcorr（x,x,M,'
biased'
%计算有偏估计自相关函数，长度为-M到M，
%共2M+1
r0=rxx（M+1）;
%r0为零点上的自相关函数，相对于-M，第M+1个点为零点
R=rxx（M+2:
2*M+1）;
%R为从1到第M个点的自相关函数矩阵
%确定矩阵大小
a=zeros（M,M）;
FPE=zeros（1,M）;
%FPE：
最终预测误差，用来估计模型的阶次
var=zeros（1,M）;
%求初值
a（1,1）=-R
（1）/r0;
%一阶模型参数
var
（1）=（1-（abs（a（1,1）））^2）*r0;
%一阶方差
FPE
（1）=var
（1）*（M+2）/（M）;
%递推
forp=2:
M
sum=0;
fork=1:
p-1%求a（p,p）
sum=sum+a（p-1,k）*R（p-k）;
end
a（p,p）=-（R（p）+sum）/var（p-1）;
p-1%求a（p,k）
a（p,k）=a（p-1,k）+a（p,p）*a（p-1,p-k）;
var（p）=（1-a（p,p）^2）*var（p-1）;
%求方差
FPE（p）=var（p）*（M+1+p）/（M+1-p）;
%求最终预测误差
end
%确定AR模型的最佳阶数
min=FPE
（1）;
%求出FPE最小时对应的阶数
p=1;
fork=2:
M
ifFPE（k）<
min
min=FPE（k）;
p=k;
%功率谱估计
W=0.01:
0.01:
pi;
%功率谱以2*pi为周期，又信号为实信号，只需输出0到PI即可；
he=ones（1,length（W））;
%length（）求向量的长度
fork=1:
p
he=he+（a（p,k）.*exp（-j*k*W））;
end
Pxx=var（p）./（（abs（he））.^2）;
%功率谱函数；
F=W*fs/（pi*2）;
%将角频率坐标换算成HZ坐标，便于观察；
重要！
figure;
plot（F,abs（Pxx））
频率/Hz'
功率谱P'
）,title（['
AR模型的最佳阶数p='
int2str（p）]）;
五、Burg法谱估计程序
%设采样频率为1
N=900;
%数据长度改变数据长度会导致分辨率的变化；
%第一个sin信号的频率,f1/fs=0.2
%第二个sin信号的频率，f1/fs=0.2或者0.3
M=512;
%滤波器阶数的最大取值,超过则认为代价太大而放弃
s=sin（2*pi*f1*n/fs）+sin（2*pi*f2*n/fs）;
x=awgn（s,10）;
fori=1:
N
ef（1,i）=x（i）;
eb（1,i）=x（i）;
sum=0;
sum=sum+x（i）*x（i）;
r
（1）=sum/N;
%Burg递推
%求解第p个反射系数
sum1=0;
forn=p:
sum1=sum1+ef（p-1,n）*eb（p-1,n-1）;
sum1=-2*sum1;
sum2=0;
sum2=sum2+ef（p-1,n）*ef（p-1,n）+eb（p-1,n-1）*eb（p-1,n-1）;
k（p-1）=sum1/sum2;
%求解预测误差平均功率
r（p）=（1-k（p-1）*k（p-1））*r（p-1）;
%求解p阶白噪声方差
q（p）=r（p）;
%系数a
ifp>
2
fori=1:
p-2
a（p-1,i）=a（p-2,i）+k（p-1）*a（p-2,p-1-i）;
a（p-1,p-1）=k（p-1）;
%求解前向预测误差
forn=p+1:
ef（p,n）=ef（p-1,n）+k（p-1）*eb（p-1,n-1）;
%求解后向预测误差
N-1
eb（p,n）=eb（p-1,n-1）+k（p-1）*ef（p-1,n）;
%计算功率谱
forj=1:
sum3=0;
sum4=0;
p-1
sum3=sum3+a（p-1,i）*cos（2*pi*i*j/N）;
sum3=1+sum3;
sum4=sum4+a（p-1,i）*sin（2*pi*i*j/N）;
end
pxx=sqrt（sum3*sum3+sum4*sum4）;
pxx=q（M）/pxx;
pxx=10*log10（pxx）;
pp（j）=pxx;
%画出功率谱
ff=1:
ff=ff/N;
plot（ff,pp）,axis（[00.5-2010]）,xlabel（'
频率'
幅度（dB）'
升级会员 