随机信号 处理各种功率谱估计方法及matlab仿真实现

一、随机信号原理分析随机信号的古典法谱估计广义平稳随机过程的功率谱是自相关函数的福利叶变换，它取决于无数个自相关函数值。

但对于许多实际问题，可以利用的数据往往是有限的，所以要准确的计算功率谱是不可能的。

比较合理的目标是设法得出功率谱的一个好的估计，这就是功率谱估计。

功率谱估计有两类大的方法：古典谱估计和现代谱估计。

古典谱估计又有相关法估计，周期图法估计，WOSO法，Bartlett法；现代谱估计有Levinson-Durbin算法和Burg算法。

下面对他们分别作介绍。

1．相关法谱估计这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。

这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N的有限长序列X N（n）。

第二步：由N长序列求（2M-1）点的自相关函数序列。

即这里，m=-(M-1),。

,-1,0,1,。

,M-1，M N，R^在此处键入公式。

x(m)是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。

,M-1的，另一半也就知道了。

第三步：由相关函数的傅式变换求功率谱。

即以上过程中经历了两次截断，一次是将x(n)截成N长，称为加数据窗，一次是将截成（2M-1）长，称为加延迟窗。

因此所得的功率谱仅是近似值，也叫谱估计，式中的，代表估值。

一般取M〈〈N，因为只有当M较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。

因此，在FFT问世之前，相关法是最常用的谱估计方法。

当FFT问世后，情况有所变化。

因为截断后的可视作能量信号，由相关卷积定理可得这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。

我们可对上式两边取（2N-1）点DFT，则有于是将时域卷积变为频域乘积。

用快速相关求的完整方案如下：1.对N长的充（N-1）个零，成为（2N-1）长的。

2.求（2N-1）点的FFT，得。

3.求。

由DFT性质，是纯实的，满足共轭偶对称，而一定是实偶的，且以（2N-1）为周期。

4.求（2N-1）点的IFFT：这里是实偶的，m=-(N-1),。

,0,。

,N-1。

本来IFFT求和范围是0--2N-2，由于的实偶性与周期性，求和范围改为-（N-1）--（N-1）不影响计算结果。

同理可将m的范围改为-（N-1）--（N-1）。

上述的快速相关中，充零的目的是为了能用圆周卷积代替线性卷积，以便进一步采用快速卷积算法。

整个相关法的计算框图见图A。

快速相关输出是-（N-1）--（N-1）的2N-1点，加窗后截取的是-（M-1）--（M-1）的，最后作（2M-1）点FFT，得。

2.周期图法周期图法又称直接法。

1.从随机信号x(n)中截取N长的一段，把它视为能量有限信号，直接取的傅式变换，得频谱，2.取其幅值的平方，并除以N，作为对x(n)真实功率谱的估计的估计，即古典谱估计的改进相关法和周期图法都是用获得的N个数据对随机过程进行功率谱估计，它隐含着对无限长数据加了一个矩形窗截断。

时域中与矩形窗函数相乘对应于频域中与矩形窗频谱相卷积，由于矩形窗频谱不同于冲击函数，所以所得结果的谱估计不同于真实的功率谱。

为了使估计逼近真实谱，应设法使窗函数谱幅度函数逼近冲击函数。

矩形窗频谱与冲击函数相比存在两方面差异，一是主瓣不是无限窄，二是有旁瓣。

由于主瓣不是无限窄，使得功率向附近扩散，使信号模糊，降低了分辨率，主瓣越宽分辨率越低。

由于旁瓣的存在又会产生两个结果，一是功率谱主瓣内的能量泄露到旁瓣使谱估计方差增大，二是与旁瓣卷积后得的功率谱完全属于干扰。

严重情况下会使弱信号淹没在强信号的干扰中，无法检测出来，这是古典谱估计的主要缺点。

为了减少这些缺点造成的影响，对古典谱估计做了改进，现在比较常用的改进方法是Bartlett法和Welch法。

下面分别作介绍。

1.Bartlett法平滑是用一个适当的窗函数与算出的功率谱进行卷积，使谱线平滑。

这种方法得出的谱估计是无偏的，方差也小，但分辨率下滑。

平均就是将截取的数据段再分成L个小段，分别计算功率谱后取功率谱的平均，这种方法叫Bartlett法。

因为L个平均的发差比随机变量的单独方差小L倍，所以当时，L哥平方的方差趋于零，可以达到一致谱估计的目的。

2．Welch法Welch法又叫加权交叠平均法，简记为WOSA法，这种方法以加窗（加权）求取平滑，以分段重叠求得平均.其主要步骤是：1.将N长的数据段分成L个小段，每小段M点，相邻小段间交叠M/2点（即2:1分段）。

因为L（M/2）+M/2=N,所以段数2.对各小段加同样的平滑窗后起傅氏变换3.用下式求各小段功率谱的平均这里，代表窗函数平均功率，MU是M长窗函数的能量。

Welch法得出的是无偏谱估计，当段数 L增大时方差减小，趋于一致估计，分段平均减小了由数据样本本身的随机性带来的方差，段数越多方差越小，但分辨率下降。

分段多了每段的点数必然少了，所以采用2:1分段，拉长了每小段长度，有利于平滑过渡。

随机信号的现代法谱估计1.自相关法自相关法是通过Yule-Walker方程获得AR模型参数。

Levinson-Durbin算法是一种解Yule-Walker方程的高效算法。

该算法用信号p+1个自相关函数值解出AR模型的P+1个参数2．Burg算法Burg算法与自相关法不同, 它是使序列x ( n)的前后向预测误差功率之和:最小。

利用Burg法求解AR模型参数的步骤:第一步:由初始条件e f0(n)=x(n)和e b0(n)=x(n),根据式,,求出反射系数k1 ;第二步:根据序列x ( n)的自相关函数,求出阶次m = 1时的AR模型参数与前后向预测误差功率之和;第三步:由式求出前向预测误差( n)与后向预测误差( n) , 然后由式估计出反射系数k2 ;第四步:由式(10) Levinsion递推关系,求出阶次m =2时的AR模型参数和 ,以及;第五步:重复上述过程, 直到阶次m = p, 这样就求出了所有阶次的AR模型参数。

二、实验程序及实验结果%时域信号clear all;n=0:511; %设置点数xn=sin(2*pi*0.1*n)+2*sin(2*pi*0.4*n)+randn(1,512); %时域信号axis([0,512,-5,5]) %坐标轴范围subplot(5,1,1);plot(n,xn) %作图xlabel('n');ylabel('xn'); %坐标轴意义title('时域信号')grid on %显示方格% 相关法谱估计X1k=fft(xn,1023); %2N-1点福利叶变换X2k=abs(X1k.^2)/512;Rm=ifft(X2k,1023); %求出相关函数m=-255:255Sxk=fft(Rm,511); %信号功率谱估计S=10*log(Sxk); %功率谱幅值（用分贝表示）subplot(5,1,2); %作图k=0:510; %频域自变量范围k=k/511; %对自变量做归一化处理plot(k,S); %显示图形xlabel('Frequency')ylabel('10lg(PSD)');title('相关法谱估计结果')grid on%周期图表法谱估计Fx1=10*log(abs(fft(xn,512).^2)/512); %周期法求功率谱密度f1=(0:length(Fx1)-1)/length(Fx1);subplot(5,1,3)plot(f1,Fx1)xlabel('Frequency/Hz');ylabel('PSD');title('周期法功率谱估计结果')grid on;%Bartlett法 (l=2）N=512for i=1:N/2 %把信号分成2段yn1(1,i)=cos(2*pi*0.1*i)+cos(2*pi*0.4*i);yn2(1,i)=cos(2*pi*0.1*(i+N/2))+cos(2*pi*0.4*(i+N/2));if i==N/2break;endendyn1=yn1+randn(1,N/2); %第1段加噪声yn2=yn2+randn(1,N/2); %第2段加噪声xk3=fft(yn1,N/2); %做FFTxk4=fft(yn2,N/2); %做FFTxk5=(abs(xk3).^2+abs(xk4).^2)./(N); %模平方取均值k=1:N/2;f=(k-1)./(N/2);subplot(5,1,4)plot(f,10*log(abs(xk5)));title('Bartlett法 (l=2）');%Welch(汉宁窗，L=7,64点交叠）M=0;for i=1:7for j=1:128 %做128次循环wh(i,j)=0.5*(1-cos(2*pi*j/(127))); %汉宁窗M=M+wh(i,j).^2/7; %窗函数平均功率yn7(i,j)=(cos(2*pi*0.1*(j+64*(i-1)))+cos(2*pi*0.4*(j+64*(i-1)))).*wh( i,j); %信号分为重叠64点的7段endendfor i=1:7yn7(i,:) = yn7(i,:)+randn(1,128); %对每段加噪声xk_7(i,:)=fft( yn7(i,:),N/4); %对每段做fftendxk8=zeros(1,128);for i=1:7xk8=(xk8+ abs(xk_7(i,:)).^2./M); %求功率谱估计endxk8=xk8/7;k=1:N/4;ff=(k-1)./(N/4);subplot(515)plot(ff,10*log(abs(xk8)));title('Welch(汉宁窗，L=7,64点交叠');clear allN=512;n=0:N-1;y=sin(2*pi*0.1*n)+sin(2*pi*0.4*n)+randn(size(n));%-----------自相关（yule-walker,levinson）算法-----------% p=20; %定义阶数a=zeros(p,p);s=zeros(1,p+1);for m=1:Nr=0;for n=1:(N-m)r=r+y(n)*y(n+m-1);endrx(m)=r/N;enda(1,1)=-rx(2)/rx(1);q(1)=rx(1)*(1-a(1,1)^2);for m=2:p;s=0;for j=1:m-1;s=s+a(m-1,j)*rx(m-j+1);endk(m)=-(rx(m+1)+s)/q(m-1);a(m,m)=k(m);for i=1:m-1;a(m,i)=a(m-1,i)+k(m)*a(m-1,m-i); q(m)=q(m-1)*(1-k(m)^2);endends=[1,a(p,:)];P=abs(fft(s,1000)).^2;Pw=10*log(q(p)./P);w=(0:N-1)/N;f=0:0.001:0.999;M=length(f);plot(f(1:M/2),Pw(1,1:M/2)); ylabel('10log(PSD)');title('自相关算法功率谱');%-----------------基于burg方法的功率谱估计-------------%N =512;p1=3;f1 = 0.1;f2 = 0.3;n= 0:(N-1);x = sin(2*pi*f1*n) + sin(2*pi*f2*n)+randn(1,N);% xcorr(x,9,'biased')%ef=zeros(p1,N);eb=zeros(p1,N);ef(1,:)=x; %-----ef(0,:)-----% eb(1,:)=x; %-----eb(0,:)-----% cov(1)=x*x'/N;k(1)=0; %-------k(0)------% a=zeros(p1+1,p1+1);for p=2:p1+1num=0;den=0;for n=p:Nnum= num+ (-2)*ef(p-1,n)*eb(p-1,n-1);den=den+(ef(p-1,n))^2+(eb(p-1,n-1))^2;endk(p)=num./den;cov(p)=(1-abs(k(p))^2).*cov(p-1);a(p,p)=k(p);for n=p:Nef(p,n)=ef(p-1,n)+k(p).*eb(p-1,n-1); eb(p,n)=eb(p-1,n-1)+k(p).*ef(p-1,n); endendk=k(1,2:p1+1);a=a(2:p1+1,2:p1+1);for p=2:p1for i=1:p-1a(p,i)=a(p-1,i)+k(p)*a(p-1,p-i);endends=a(p1,:);s=[1,s];n=(0:511)/512;y=fft(s,512);figure(3);plot(n,cov(p1+1)./abs(y).^2);title('基于burg方法的功率谱估计');。