功率谱密度估计方法的MATLAB实现文档格式.docx

功率谱密度估计方法的MATLAB实现文档格式.docx

《功率谱密度估计方法的MATLAB实现文档格式.docx》由会员分享，可在线阅读，更多相关《功率谱密度估计方法的MATLAB实现文档格式.docx（12页珍藏版）》请在冰豆网上搜索。

%采样频率100kHzN=1024;

%数据长度N=1024n=0:

N-1;

t=n/Fs;

xn=sin（2000*2*pi*t）;

%正弦波,f=2000HzY=awgn（xn,10）;

%加入信噪比为10db的高斯白噪声subplot（2,1,1）;

plot（n,Y）title（信号）xlabel（时间）;

ylabel（幅度）;

gridon;

window=boxcar（length（xn）;

%矩形窗nfft=N/4;

%采样点数Pxxf=periodogram（Y,window,nfft,Fs）;

%直接法subplot（2,1,2）;

plot（f,10*log10（Pxx）;

title（周期图法谱估计,int2str（N）,点）;

xlabel（频率（Hz））;

ylabel（功率谱密度）;

2、仿真结果二、修正周期图法（加窗）谱估计程序1、源程序Fs=100000;

%采样频率100kHzN=512;

%数据长度M=32;

%汉明窗宽度n=0:

subplot（2,1,1）;

window=hamming（M）;

%汉明窗Pxxf=pwelch（Y,window,10,256,Fs）;

subplot（2,1,2）;

title（修正周期图法谱估计N=,int2str（N）,M=,int2str（M）;

2、仿真结果三、最大熵法谱估计程序1、源程序fs=1;

%设采样频率N=128;

%数据长度改变数据长度会导致分辨率的变化；

f1=0.2*fs;

%第一个sin信号的频率,f1/fs=0.2f2=0.3*fs;

%第二个sin信号的频率，f2/fs=0.2或者0.3P=10;

%滤波器阶数n=1:

N;

s=sin（2*pi*f1*n/fs）+sin（2*pi*f2*n/fs）;

%s为原始信号x=awgn（s,10）;

%x为观测信号，即对原始信号加入白噪声，信噪比10dBfigure

（1）;

%画出原始信号和观测信号subplot（2,1,1）;

plot（s,b）,xlabel（时间）,ylabel（幅度）,title（原始信号s）;

grid;

plot（x,r）,xlabel（时间）,ylabel（幅度）,title（观测信号x）;

Pxx1,f=pmem（x,P,N,fs）;

%最大熵谱估计figure

（2）;

plot（f,10*log10（Pxx1）;

xlabel（频率（Hz）;

ylabel（功率谱（dB）;

title（最大熵法谱估计模型阶数P=,int2str（P）,数据长度N=,int2str（N）;

2、仿真结果四、Levinson递推法谱估计程序1、源程序fs=1;

%设采样频率为1N=1000;

%第二个sin信号的频率，f1/fs=0.2或者0.3M=16;

%滤波器阶数的最大取值,超过则认为代价太大而放弃L=2*N;

%有限长序列进行离散傅里叶变换前，序列补零的长度n=1:

plot（s,b）,axis（0100-33）,xlabel（时间）,ylabel（幅度）,title（原始信号s）;

plot（x,r）,axis（0100-33）,xlabel（时间）,ylabel（幅度）,title（观测信号x）;

%计算自相关函数rxx=xcorr（x,x,M,biased）;

%计算有偏估计自相关函数，长度为-M到M，%共2M+1r0=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:

Msum=0;

fork=1:

p-1%求a（p,p）sum=sum+a（p-1,k）*R（p-k）;

enda（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）;

endvar（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:

MifFPE（k）2fori=1:

p-2a（p-1,i）=a（p-2,i）+k（p-1）*a（p-2,p-1-i）;

endenda（p-1,p-1）=k（p-1）;

%求解前向预测误差forn=p+1:

Nef（p,n）=ef（p-1,n）+k（p-1）*eb（p-1,n-1）;

end%求解后向预测误差forn=p:

N-1eb（p,n）=eb（p-1,n-1）+k（p-1）*ef（p-1,n）;

endend%计算功率谱forj=1:

Nsum3=0;

sum4=0;

fori=1:

p-1sum3=sum3+a（p-1,i）*cos（2*pi*i*j/N）;

endsum3=1+sum3;

p-1sum4=sum4+a（p-1,i）*sin（2*pi*i*j/N）;

endpxx=sqrt（sum3*sum3+sum4*sum4）;

pxx=q（M）/pxx;

pxx=10*log10（pxx）;

pp（j）=pxx;

end%画出功率谱ff=1:

ff=ff/N;

figure;

plot（ff,pp）,axis（00.5-2010）,xlabel（频率）,ylabel（幅度（dB））,title（功率谱P）;

2、仿真结果