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

1、功率谱密度预计方法的MATLAB实此刻应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用观点，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号往常是波的形式，比如电磁波、随机振动或许声波。当波的频谱密度乘以一个适合的系数后将获得每单位频次波携带的功率，这被称为信号的功率谱密度（powerspectraldensity,PSD）或许谱功率散布（spectralpowerdistribution,SPD）。功率谱密度的单位往常用每赫兹的瓦特数（W/Hz）表示，或许使用波长而不是频次，即每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的安稳过程

2、的时候才存在。假如信号不是安稳过程，那么自有关函数必定是两个变量的函数，这样就不存在功率谱密度，可是能够使用近似的技术预计时变谱密度。信号功率谱的观点和应用是电子工程的基础，特别是在电子通讯系统中，比如无线电和微波通讯、雷达以及有关系统。所以学习怎样进行功率谱密度预计十分重要，借助于Matlab工具能够实现各样谱预计方法的模拟仿真并输出结果。下边对周期图法、修正周期图法、最大熵法、Levinson递推法和Burg法的功率谱密度预计方法进行程序设计及仿真并给出仿真结果。以下程序运转平台：MatlabR2015a（）一、周期图法谱预计程序1、源程序Fs=100000;%采样频次100kHzN=10

3、24;%数据长度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(幅度);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(周期图法谱预计,in

4、t2str(N),点);xlabel(频次（Hz）);ylabel(功率谱密度);2、仿真结果二、修正周期图法（加窗）谱预计程序1、源程序Fs=100000;%采样频次100kHzN=512;%数据长度M=32;%汉明窗宽度n=0:N-1;2/12t=n/Fs;xn=sin(2000*2*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=p

5、welch(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.2f2=0.3*fs;%第二个sin信号的频次，f2/fs=0.2或许0.3P=10;%滤波器阶数n=1:N;s=sin(2*p

6、i*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(幅度),title(观察信号x);Pxx1,f=pmem(x,P,N,fs);%最大熵谱预计figure(2);plot(f,10*log10(Pxx1);xlabel(频次(Hz);yl

7、abel(功率谱(dB);title(最大熵法谱预计模型阶数P=,int2str(P),数据长度N=,int2str(N);4/122、仿真结果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*

8、pi*f2*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(观察信号x);grid;%计算自有关函数rxx=xcorr(x,x,M,biased);%计算有偏预计自有关函数，长度为-M到M，%共2M+1r0

9、=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:Msum=0;fork=1:p-1%求a(p,p)sum=sum+a(p-1,k)*R(p-k);enda(p,p)

10、=-(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)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);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