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

1、功率谱密度估计方法的MATLA破现在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度(powerspectraldensity,PSD)或者谱功率分布(spectralpowerdistribution,SPD)。功率谱密度的单位通常用每赫兹的瓦特数(W/Hz)表示,或者使用波长而不是频率，即每纳米的瓦特数(W/nm)来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的

2、时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中，例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要，借助于Matlab工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大嫡法、Levinson递推法和Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。周期图法谱估计程序1、源程序Fs=100000;N=1024;n=0:N-1;t=n/Fs;xn=sin(2000*2*pi*t)

3、;Y=awgn(xn,10);subplot(2,1,1);以下程序运行平台：MatlabR2015a(8.5.0.197613)%采样频率100kHz%数据长度N=1024%正弦波，f=2000Hz%加入信噪比为10db的高斯白噪声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(

4、f,10*log10(Pxx);gridon;title('周期图法谱估计,',int2str(N),'点');xlabel('频率(Hz)');ylabel('功率谱密度,);2、仿真结果,伯号-2020040060080010001200时间20周期图法谓估计,1024点TO-60-8000.511.522.533.544.55频率(Hz)x104修正周期图法(加窗)谱估计程序1、源程序Fs=100000;%采样频率100kHzN=512;%数据长度M=32;%汉明窗宽度n=0:N-1;2/12t=n/Fs;xn=sin(2000*2

5、*pi*t);%正弦波,f=2000HzY=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

6、=',int2str(M);xlabel('频率(Hz)');ylabel('功率谱密度');2、仿真结果_信号1-12010020030。400500600时间“修正周期图法谱估计N=512M=3230MO-50-6000.511.522.533.544.55顺率HzJxio43/12三、最大嫡法谱估计程序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=

7、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;subplot(2,1,2);plot(x,"r"),xlabel("时间"),ylabel("幅度"),

8、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、仿真结果原始信号S21O-1-220406080100120140时间观测信号X42O-2-4«蹙2040608010012。140时间(mp)去能会

9、-5最大燧法谱佶计模型阶数P=10数据氏度N=128O0,050.10.150.20,250.30.350,40,450.5频率（H,5/12四、Levinson递推法谱估计程序1、源程序fs=1;%设采样频率为1N=1000;%数据长度改变数据长度会导致分辨率的变化；f1=0.2*fs;%第一个sin信号的频率,f1/fs=0.2f2=0.3*fs;%第二个sin信号的频率，f1/fs=0.2或者0.3M=16;%滤波器阶数的最大取值，超过则认为代价太大而放弃L=2*N;%有限长序列进行离散傅里叶变换前，序列补零的长度n=1:N;s=sin(2*pi*f1*n/fs)+sin(2*pi*f2

10、*n/fs);%s为原始信号x=awgn(s,10);%x为观测信号，即对原始信号加入白噪声，信噪比10dBfigure(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('观测彳

11、言号x');grid;%计算自相关函数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)A2)*r0;%一阶方差6/12%递推fo

12、rp=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);fork=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)A2)*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)<minmin=FPE(k);p=k;endend%

13、功率谱估计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);endPxx=var(p)./(abs(he).A2);%功率谱函数；F=W*fs/(pi*2);%将角频率坐标换算成HZ坐标，便于观察；重要！figure;plot(F,abs(Pxx)xlabel('频率/Hz'),ylabel('功率谱P'),title('AR模型的最佳阶数p='int2st

14、r(p);grid;2、仿真结果8/12ieAR膜型的最过阶款60SOJ0.150.20.250.30,350.40.450.5庭j率,Hn五、Burg法谱估计程序1、源程序fs=1;%设采样频率为1N=900；喊据长度改变数据长度会导致分辨率的变化；f1=0.2*fs;%第一个sin信号的频率,f1/fs=0.2f2=0.3*fs;%第二个sin信号的频率，f1/fs=0.2或者0.3M=512;%8波器阶数白最大取值,超过则认为代价太大而放弃n=1:N;s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%s为原始信号x=awgn(s,10);%x为观测信号，即对原

15、始信号加入白噪声，信噪比10dBfori=1:Nef(1,i)=x(i);eb(1,i)=x(i);endsum=0;fori=1:Nsum=sum+x(i)*x(i);9/12endr(1)=sum/N;%Burg递推forp=2:M%求解第p个反射系数sum1=0;forn=p:Nsum1=sum1+ef(p-1,n)*eb(p-1,n-1);endsum1=-2*sum1;sum2=0;forn=p:Nsum2=sum2+ef(p-1,n)*ef(p-1,n)+eb(p-1,n-1)*eb(p-1,n-1);endk(p-1)=sum1/sum2;%求解预测误差平均功率r(p)=(1-k

16、(p-1)*k(p-1)*r(p-1);%求解p阶白噪声方差q(p)=r(p);%系数aifp>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);10/12end%求解后向预测误差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;fori=1:p-1sum4=sum4+a(p-1,i)*sin(2*pi*i*j/N);endpxx=sqrt(sum3*sum3+sum4*sum4);pxx=q(M)/pxx;pxx