功率谱密度估计方法的MATLAB实现Word文件下载.docx
- 文档编号:18227807
- 上传时间:2022-12-14
- 格式:DOCX
- 页数:13
- 大小:119.25KB
功率谱密度估计方法的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)'
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 功率 密度 估计 方法 MATLAB 实现
![提示](https://static.bdocx.com/images/bang_tan.gif)