完整word版功率谱密度估计方法的MATLAB实现Word格式文档下载.docx

完整word版功率谱密度估计方法的MATLAB实现Word格式文档下载.docx

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

Y=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）,'

点'

]）;

频率（Hz）'

功率谱密度'

2、仿真结果

二、修正周期图法（加窗）谱估计程序

N=512;

%数据长度

M=32;

%汉明窗宽度

window=hamming（M）;

%汉明窗

[Pxxf]=pwelch（Y,window,10,256,Fs）;

修正周期图法谱估计N='

M='

int2str（M）]）;

三、最大熵法谱估计程序

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）'