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

1、功率谱密度估计方法的功率谱密度估计方法的 MATLAB 实现实现 功率谱密度估计方法的 MATLAB实现 在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density,PSD）或者谱功率分布（spectral power distribution,SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即

2、每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中，例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要，借助于 Matlab 工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大熵法、Levinson 递推法和 Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。以下程序运行平台：Matlab R2015

3、a（8.5.0.197613）一、周期图法谱估计程序 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(幅度);grid on;window=boxcar(length(xn);%矩形窗 nfft=N/4;%采样点数 Pxx f=periodogram(Y,window,nfft,Fs

4、);%直接法 subplot(2,1,2);plot(f,10*log10(Pxx);grid on;title(周期图法谱估计,int2str(N),点);xlabel(频率（Hz）);ylabel(功率谱密度);2、仿真结果 二、修正周期图法（加窗）谱估计程序 1、源程序 Fs=100000;%采样频率 100kHz N=512;%数据长度 M=32;%汉明窗宽度 n=0:N-1;t=n/Fs;xn=sin(2000*2*pi*t);%正弦波,f=2000Hz Y=awgn(xn,10);%加入信噪比为 10db的高斯白噪声 subplot(2,1,1);subplot(2,1,1);pl

5、ot(n,Y)title(信号)xlabel(时间);ylabel(幅度);grid on;window=hamming(M);%汉明窗 Pxx f=pwelch(Y,window,10,256,Fs);subplot(2,1,2);plot(f,10*log10(Pxx);grid on;title(修正周期图法谱估计 N=,int2str(N),M=,int2str(M);xlabel(频率（Hz）);ylabel(功率谱密度);2、仿真结果 三、最大熵法谱估计程序 1、源程序 fs=1;%设采样频率 N=128;%数据长度 改变数据长度会导致分辨率的变化；f1=0.2*fs;%第一个 s

6、in 信号的频率,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

7、(幅度),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;%设采样频率为 1 N=1000;%数据长度 改变数据长度会导致分辨率的变化；f1=0.2*fs;%第一个 sin 信号的频率,f1/fs=0.2 f2=0.3*fs;%第二个 sin 信号的频率，f

8、1/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(0 100-3 3),xlabel(时间),ylabel(幅度),title(原始信号 s);grid;subplot(2,1,2);plot(x,

9、r),axis(0 100-3 3),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

10、;%一阶模型参数 var(1)=(1-(abs(a(1,1)2)*r0;%一阶方差 FPE(1)=var(1)*(M+2)/(M);%递推 for p=2:M sum=0;for k=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);for k=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%确定 A

11、R 模型的最佳阶数 min=FPE(1);%求出 FPE 最小时对应的阶数 p=1;for k=2:M if FPE(k)2 for i=1:p-2 a(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);%求解前向预测误差 for n=p+1:N ef(p,n)=ef(p-1,n)+k(p-1)*eb(p-1,n-1);end%求解后向预测误差 for n=p:N-1 eb(p,n)=eb(p-1,n-1)+k(p-1)*ef(p-1,n);end end%计算功率谱 for j=1:N sum3=0;sum4=0;for i=1:p-1 sum3=sum3+a(p-1,i)*cos(2*pi*i*j/N);end sum3=1+sum3;for i=1:p-1 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;end%画出功率谱 ff=1:N;ff=ff/N;figure;plot(ff,pp),axis(0 0.5-20 10),xlabel(频率),ylabel(幅度（dB）),title(功率谱 P);grid;2、仿真结果