功率谱密度估计方法计划MATLAB实现

功率谱密度估计方法计划MATLAB实现功率谱密度估计方法计划MATLAB实现功率谱密度估计方法计划MATLAB实现

功率谱密度估计方法的MATLAB实现

在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用看法，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号平时是波的形式，比方电磁波、随机振动也许声波。当波的频谱密度乘以一个合适的系数后将获取每单位频率

波携带的功率，这被称为信号的功率谱密度（powerspectraldensity,PSD）也许谱功率分布（spectralpowerdistribution,SPD）。功率谱密度的单位平时用每赫兹的瓦特数（W/Hz）表示，也许使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。若是信号不是平稳过程，那么自相关函数必然是两个变量的函数，这样就不存在功率谱密度，但是可以使用近似的技术估计时变谱密度。信号功率谱的看法和应用是电子工程的基础，特别是在电子通信系统中，比方无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要，借助于Matlab工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大熵法、

Levinson递推法和Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。

以下程序运行平台：MatlabR2015a（）

一、周期图法谱估计程序

1、源程序

Fs=100000;%采样频率100kHz

N=1024;%数据长度N=1024

n=0:N-1;

t=n/Fs;

xn=sin(2000*2*pi*t);%正弦波,f=2000Hz

Y=awgn(xn,10);%加入信噪比为10db的高斯白噪声

subplot(2,1,1);

plot(n,Y)

title('信号')

xlabel('时间');ylabel('幅度');1/12gridon;

window=boxcar(length(xn));%矩形窗

nfft=N/4;%采样点数

[Pxxf]=periodogram(Y,window,nfft,Fs);%直接法

subplot(2,1,2);

plot(f,10*log10(Pxx));

gridon;

title(['周期图法谱估计,',int2str(N),'点']);

xlabel('频率（Hz）');ylabel('功率谱密度');

2、仿真结果

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

1、源程序

Fs=100000;%采样频率100kHz

N=512;%数据长度

M=32;%汉明窗宽度

n=0:N-1;2/12t=n/Fs;

xn=sin(2000*2*pi*t);%正弦波,f=2000Hz

Y=awgn(xn,10);%加入信噪比为10db的高斯白噪声

subplot(2,1,1);subplot(2,1,1);

plot(n,Y)

title('信号')

xlabel('时间');ylabel('幅度');

gridon;

window=hamming(M);%汉明窗

[Pxxf]=pwelch(Y,window,10,256,Fs);

subplot(2,1,2);

plot(f,10*log10(Pxx));

gridon;

title(['修正周期图法谱估计N=',int2str(N),'M=',int2str(M)]);

xlabel('频率（Hz）');ylabel('功率谱密度');

2、仿真结果3/12三、最大熵法谱估计程序

1、源程序

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);%画出原始信号和观察信号

subplot(2,1,1);

plot(s,'b'),xlabel('时间'),ylabel('幅度'),title('原始信号s');

grid;

subplot(2,1,2);

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)]);4/122、仿真结果

5/12四、Levinson递推法谱估计程序

1、源程序

fs=1;%设采样频率为1

N=1000;%数据长度改变数据长度会以致分辨率的变化；

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

f2=0.3*fs;%第二个sin信号的频率，f1/fs=0.2也许0.3

M=16;%滤波器阶数的最大取值,高出则认为代价太大而放弃

L=2*N;%有限长序列进行失散傅里叶变换前，序列补零的长度

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);%画出原始信号和观察信号

subplot(2,1,1);

plot(s,'b'),axis([0100-33]),xlabel('时间'),ylabel('幅度'),title('原始信号s');

grid;

subplot(2,1,2);

plot(x,'r'),axis([0100-33]),xlabel('时间'),ylabel('幅度'),title('观察信号x');

grid;

%计算自相关函数

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;%一阶方差6/12FPE(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);

fork=1:p-1%求a(p,k)

a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);

end

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;

end

end

%功率谱估计

W=0.01:0.01:pi;%功率谱以2*pi为周期，又信号为实信号，只需输出0到PI即可；

he=ones(1,length(W));%length()求向量的长度

fork=1:p7/12he=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))

xlabel('频率/Hz'),ylabel('功率谱P'),title(['AR模型的最正确阶数p='int2str(p)]);

grid;

2、仿真结果

8/12

五、Burg法谱估计程序

1、源程序

fs=1;%设采样频率为1

N=900;%数据长度改变数据长度会以致分辨率的变化；

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

f2=0.3*fs;%第二个sin信号的频率，f1/fs=0.2也许0.3

M=512;%滤波器阶数的最大取值,高出则认为代价太大而放弃

n=1:N;

s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%s为原始信号

x=awgn(s,10);%x为观察信号，即对原始信号加入白噪声，信噪比10dB

fori=1:N

ef(1,i)=x(i);

eb(1,i)=x(i);

end

sum=0;

fori=1:N

sum=sum+x(i)*x(i);

9/12end

r(1)=sum/N;

%Burg递推

forp=2:M

求解第p个反射系数

sum1=0;

forn=p:N

sum1=sum1+ef(p-1,n)*eb(p-1,n-1);

end

sum1=-2*sum1;

sum2=0;

forn=p:N

sum2=sum2+ef(p-1,n)*ef(p-1,n)+eb(p-1,n-1)*eb(p-1,n-1);

end

k(p-1)=sum1/sum2;

求解展望误差平均功率

r(p)=(1-k(p-1)*k(p-1))*r(p-1);

求解p阶白噪声方差

q(p)=r(p);

系数aifp>2

fori=1:p-2a(p-1,i)=a(p-2,i)+k(p-1)*a(p-2,p-1-i);

end

end

a(p-1,p-1)=k(p-1);

求解前向展望误差forn=p+1:N

ef(p,n)=ef(p-1,n)+k(p-1)*eb(p-1,n-1);

10/12end

%求解后向展望误差

forn=p:N-1

eb(p,n)=eb(p-1,n-1)+k(p-1)*ef(p-1,n);

end

end

计算功率谱

forj=1:N

sum3=0;

sum4=0;

fori=1:p-1

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

end

sum3=1+sum3;

fori=1:p-1

sum4=sum4+a(p-1