多种功率谱估计郑州大学随机信号处理大作业

随机信号处理大作业

多种功率谱估计的算法实现及性能比较

一、引言

频谱分析是信号处理的基石，为我们提供了时域以外的另一种信号研究手段——频域，使得很多在时域看起来很复杂的问题，用频域来分析就变得十分简单。对于随机信号而言，由于不存在傅里叶变换，我们通过对其功率谱的分析来研究其频域特性。功率谱估计问题就是根据一组有限观测值来估计该过程谱的内容，对于平稳随机过程而言，所有的功率谱估计方法都是根据有限的观测值来逼近真实值，估计结果的好坏与估计方法密切相关。功率谱估计的方法可分为古典法和现代法，古典法基于傅里叶变换，包括直接法和间接法，现代谱估计包括直接解Yule-Walker方程法、Levinson-Durbin快速递推法、Burg算法、MUSIC算法、本文将对上述功率谱估计的方法进行分析。

二、原理及过程

1、古典法

这里采用古典法中的直接法（周期图法）进行功率谱估计，其具体步骤如下。

第一步：由获得的N点数据构成的有限长序列直接求傅里叶变换，得频谱

（1.1）

第二步：取频谱幅度的平方，并处以N,以此作为对真实功率谱的估计，

即

（1.2）

2、Yule-Walker方程法

①假定所研究的随机过程是由一白噪声序列激励一因果稳定的可逆线性系统的输出

②由观测获得的数据记录估计的参数

③由的参数估计的功率谱

由上可知，可以将平稳随机信号的功率谱表示为

（2.1）

其中，是白噪声的功率谱（为常数），是系统的频谱。这样

谱估计问题就转化为模型参数的估计问题，在AR、MA和ARMA三种模型中，求AR模型的参数是解线性方程，易于求解，并且MA模型和ARMA模型都可以用高阶的AR模型近似，所以这里我们采用AR模型来进行功率谱估计。

阶AR模型的系统函数为

（2.2）

阶AR模型有+1个待定参数：，，…,和系统增益G。

自相关函数，（2.3）“*”代表取共轭。

Yule-Walker方程：

（2.4）

可表示成下面的矩阵形式：

（2.5）

上式用到了自相关函数的偶对称性质，由这个方程，可以求出个参数。有了参数()，就可以根据尤勒-沃克方程有自相关函数和参数求系统增益G。然

后可以根据（2.6）求功率谱，或者由用freqz()函数来求功率谱

下面进行功率谱估计的关键就成了解Yule-Walker方程，下面可分为两种方法解方程。

⑴、直接解尤勒-沃克方程

①由式(2.3)计算自相关函数

②根据式(2.5)列写矩阵方程

③通过矩阵求逆解矩阵方程，得()

④求系统增益G

⑤代入公式求功率谱

⑵、Levinson-Durbin快速递推法

Yule-Walker方程是线性方程组，它的一般解法是矩阵求逆或高斯消去法，通常避免矩阵求逆，因其运算量太大，在数量级，Levinson-Durbin快速递推法是解Yule-Walker方程的快速有效的算法，这种算法利用方程组系数矩阵所具有的一系列好的性质，使运算量大大减少，在数量级，这是一种按阶次进行递推的算法，Levinson-Durbin快速递推法是从一阶开始，由阶模型的参数递推求解阶模型的参数。

M阶AR模型参数的Levinson-Durbin快速递推算法表示如下：

（2.7）

其中，，为反射系数，为预测误差功率

自相关函数估计：

（2.8）

Levinson-Durbin快速递推法解Yule-Walker方程的步骤：

①由式(2.8)计算自相关函数R(0)

②解一阶时的Yule-Walker方程得一阶参数

③由Levinson-Durbin快速递推式，即式(2.7)逐次推出各参数

④将各参数带入式（2.6），或由freqz()函数和即可估计的功率谱

3、Burg算法

Burg算法是不需要自相关函数直接由观测数据求解反射系数的方法，其特点是不对自相关函数进行估计，而是利用前、后向线性预测系数之间的递推关系，直接求出反射系数，再利用Levinson关系式求得AR模型参数。

P阶前向预测表示为：

P阶前向预测误差表示为：

（3.1）

P阶后向预测表示为：

P阶后向预测误差表示为：

（3.2）

借助Levinson关系式，

分别代入式（3.1）、（3.2）导出预测误差的递推公式：

计算前向预测误差的递推公式：

计算后向预测误差的递推公式：

递推的初始值为，即零阶预测时预测误差等于信号值，因为零阶预测相当于信号直接通过。

一般将递推过程中不断变化的阶次用表示，于是可将计算阶预测误差的递推公式表示如下：

（3.3）

Burg递推算法是直接从已知的数据序列计算参数反射系数的方法，估计反射系数依据的准则是，使前向预测误差功率和后向预测误差功率的平均值达到最小。Burg 算法是将前向预测误差功率和后向预测误差功率的算术平均作为平均预测误差。

式中中各项由式（3.2）推出。

所得反射系数估计公式：

（3.4）

再代入Levinson关系式，即可得AR模型参数。

Brug算法的步骤：

①由初始条件和式（3.4）求出

②

③由和式（3.3）求出，再由式（3.4）估计

④由式（2.7）的Levinson递推式求出m=2时的

⑤重复③、④过程，求出所有所有阶次的AR参数

4、MUSIC算法（多信号分类法）

因为信号子空间与噪声子空间相互正交，所以信号向量

都正交，与它们的线性组合也正交，有：

式中为个正弦信号的频率。令

当时，应有

即

因此，可以定义一种类似与功率谱的函数

若令则

上式取峰值的个就是个正弦信号频率的估计

基本MUSIC算法的步骤如下：

1，求相关矩阵R

2，对相关矩阵进行特征值分解，并计算信号特征值的个数

3，将参数代入

5、信噪比（SNR）与均方误差（MSE）的关系

人为的设置信噪比范围为-15db~15db，通过分别调用上面的求功率谱的函数，求特定信噪比下1000次功率谱估计的均方误差，分别绘制采用上述功率谱估计时，均方误差与信噪比的关系。

求信噪比与均方误差的关系：

①求古典法在-15db下的1000次功率谱估计的均方误差

②同上求-14db~15db下的1000次功率谱估计的均方误差

③绘制均方误差与信噪比的关系

④同理，求出采用上文中现代谱估计时，均方误差与信噪比的关系

结果分析

原始信号为，噪声为，

为均值为0，方差为1的加性高斯白噪声，待观测信号为，

分别为0.2，0.213，要求求出待测信号的128点功率谱密度，

信噪比和均方误差的关系，并分析各种功率谱估计方法的性能优劣。

1，古典法

N不小于130时可以识别出两个波峰

n=128 f1=0.2000 n=130 f1=0.2000 f2=2.2126

2，直接解yole-walker方程

P=10时，无论N取多少，均不能识别两个波峰

P=10 N=128 f1=0.2020 P=10 N=256 f1=0.2015 P=20时，N=480以上，可以识别出两个波峰

P=20 N=128 f1=0.1992 P=20 N=256 f1=0.1995 P=40时，N取65以上均可识别出两个波峰

P=40 N=128 f1=0.1998 f2=0.2122 p=40 n=256 f1=0.1995 f2=0.2135 3、levinson-durbin递推解尤勒-沃克方程

P=10时，无论N取多少，均不能识别两个波峰

P=10 N=128 f1=0.2003 P=10 N=256 f1=0.2031 P=20时，N=450以上时可以识别两个波峰

P=20 N=128 f1=0.1992 P=20 N=256 f1=0.1995 P=40时，N取65以上均可识别出两个波峰

P=40 N=128 f1=0.1998 f2=0.2140 P=40 N=256 f1=0.1995 f2=0.2135 4Burg算法

P=10时，N取任何值都不能识别出两个波峰

P=10 N=128 f1=0.2031 P=10 N=256 f1=0.2027 P=20时，N不小于55时可识别出两个波峰

P=20 N=128 f1=0.2003 f2=0.2138 P=20 N=256 f1=0.2005 f2=0.2125 P=40时，N大于50时可识别出两个波峰

P=40 N=128 f1=0.2005 f2=0.2138 P=40 N=256 f1=0.2008 f2=0.2130 5MUSIC算法

F1=0.2031 f2=0.2148 f1=0.2031 f2=0.2148

结果分析:

在相同的条件下，f1=0.2,f2=0.213,N=128时,观察图形发现，用古典法只能估计一个频率点，用AR模型估计时，当阶次P达到一定时，仍能分辨出两频率点，所以只要选择合适的阶次，就能用AR模型估计两接近的频率点，而N越大，图像越精细，越易于分辨两频率点，对于古典法，当N合适时，仍能分辨相近的频率点。由图易知，分辨率的关系为：古典法﹤直接解Yole-walker方程法﹤快速递推法

基于AR模型的现代功率谱估计质量明显优于经典谱估计，基于AR模型的现代谱估计是逐步改进，性能递增，精度提高的，直接解Yole-Walker方程运算量大，且不一定有解，采用Levinson递推可以大大减小运算量，但由于自相关函数的估计默认取样点以外的值为零，引入了一定的误差，而Burg算法引入前后向预测误差，来估计功率谱，避免了自相关函数带来的误差，因而精度较高，因此，Burg算法求得的AR模型最稳定，而且Burg算法不需要自相关函数，所以性能优于自相关法，由图可知，在P=20,N=55时即可分辨出两个峰值，所以在短数据时，Burg算法优势明显，具有较高的谱分辨率。

但是Burg算法进行谱估计会出现谱线分裂、谱峰偏移等问题，而且对于低信噪比情况，用AR模型很难准确估计出淹没在噪声中的正弦波频率，而改进的MUSIC算法采用特征分解技术，在估计正弦信号和噪声叠加的信号中明显优于AR模型功率谱估计。6信噪比与均方误差的关系

心得体会：

起初在刚接触这门课时，甚至有些反感老师的严格要求，大概是习惯了其他老师的“温柔乡”吧，现在回过头，在开始写这篇报告时，对老师有一种说不出的感觉，我想那就是感谢的羞涩表达吧。虽然老师的表达有时不靠谱，但老师比那些看似靠谱的老师要强很多，他牢记灌输知识的使命，又深谙人心之道，不时来点心灵鸡汤，努力的鞭挞我们前进，将知识灌输给我们，这让我不由得想起了：“小时候，不爱吃饭的我被妈妈训斥着吃饭的事”，想到这，对老师有一种莫名的亲切，简直深得我心。天啊，我也不知道要写这些，也许，这tm是发自肺腑，不由而衷吧。

首先，让我对Matlab和Matlab的语法有了更深的了解，提高了我的编程能力，虽然以前开设过Matlab的实验课，却从没有真正入门Matlab,可以说，陈老师是我Matlab的启蒙老师啊！

然后，让我真正理解了一部分信号处理信号分析的知识，加深了对概念的理解，以前学信号分析信号处理时，只是记忆公式和概念，很少去分析为什么，甚至也不去管为什么，陈老师是一位真正有学识的人，他不照本宣科，总能说出自己的理解，而这些很多是其他老师未曾提到的（或许是我没听到其他老师提），这对于融会贯通整个学科至关重要。

最后，让我对一些问题有了更深层次的看法。

程序附录

1经典法(周期图法)

N=130;

n=0:(N-1);

w=n/N;

wn=randn(1,N);

xn=10*sin(2*pi*0.2*n+(pi/3))+5*sin(2*pi*0.213*n+(pi/4))+wn;

Xk=fft(xn);

psd=abs(Xk).*abs(Xk);

subplot(2,1,1);

plot(n,xn);

grid;

xlabel('n');

title('两个正弦信号与白噪声叠加的时域波形');

subplot(2,1,2);

plot(w,psd);

title('经典法功率谱估计');

grid;

xlabel('频率');

ylabel('幅度');

[pks,locs]=findpeaks(Hf,'minpeakheight',10);

f=(locs-1)/(2*N)

2直接解Yole-Walker方程

N=256;

n=1:N;

p1=40;

f1=0.2;

f2=0.213;

x=10*sin(2*pi*f1*n+pi/3)+5*sin(2*pi*f2*n+pi/4)+randn(1,N);

subplot(2,1,1);

plot(n,x);

title('两个正弦信号与白噪声叠加的时域波形');%产生两个正弦信号与白噪声的叠加图像R0=0;sum=0;

for n=1:N

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

end

R0=sum/N;

sum=0;

for m=1:N

R(m)=0;

for n=1:N-m

sum=x(n)*x(n+m)+sum;

end

R(m)=sum/N;

sum=0;

end

%求自相关函数

temp=0;

for m=1:p1

temp=[temp R(m)];

A=temp;

end

A(:,1)=[];%删除第一行,求[R1；R2；R3；......Rp]

a=zeros(p1,p1);

for i=2:p1

for k=1:i-1

a(i,k)=R(i-k);

end

end%a为p行p列的矩阵，对角线及以上部分为0

B=a+a';%topelite阵的性质

for i=1:p1

B(i,i)=R0;

end %B为自相关系数矩阵

X=B\(-A');%解矩阵方程，X为[a1;12;13;......ap]

b=X'*A';

G=(R0+b)^(1/2);

C=[G zeros(1,p1)];

D=[1 X'];

[H,w]=freqz(C,D,2000);

Hf=abs(H);

subplot(212);

plot(w/(2*pi),Hf);

title('直接解Yole-Walker方程的功率谱估计');

[pks,locs]=findpeaks(Hf,'minpeakheight',10);

f=(locs-1)/(2*2000)

3 Levinson-Durbin快速递推法

N=256;n=1:N;p1=40;

f1=0.2;f2=0.213;

x=10*sin(2*pi*f1*n+pi/3)+5*sin(2*pi*f2*n+pi/4)+randn(1,N);

subplot(2,1,1);

plot(n,x);

title('两个正弦信号与白噪声叠加的时域波形');%产生两个正弦信号与白噪声的叠加图像R0=0;sum=0;

for n=1:N

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

end

R0=sum/N;

sum=0;

for m=1:N

R(m)=0;

for n=1:N-m

sum=x(n)*x(n+m)+sum;

R(m)=sum/N;

sum=0;

end

%求自相关函数

aa(1,1)=-R(1)/R0;

p(1)=R0*(1-abs(aa(1,1))^2); %尤勒-沃克方程一阶系数

m=2;sum=0;

for m=2:p1

for i=1:m-1

sum=sum+aa(m-1,i)*R(m-i);

end

aa(m,m)=-(R(m)+sum)/p(m-1); %定义levinson-durbin递推式中的km

p(m)=p(m-1)*(1-aa(m,m)^2); %定义levinson-durbin递推式中的pm

for i=1:m-1

aa(m,i)=aa(m-1,i)+aa(m,m)*aa(m-1,m-i); %定义levinson-durbin递推式中的am(i)

end

sum=0;

end

%由H(z)用freqz()求解功率谱

G=p(p1)^(1/2);

B=[G,zeros(1,p1)];

A=[1,aa(p1,:)];

[H,w]=freqz(B,A,2000);%将pi分成2000份

Hf=abs(H);

subplot(2,1,2);

plot(w/(2*pi) ,Hf);

title('基于levinson-durbin算法的功率谱估计');

[pks,locs]=findpeaks(Hf,'minpeakheight',10);

f=(locs-1)/(2*2000)

4 Brug算法

N=256; %定义采样点数

f1=0.2;f2=0.213;

P=20;%滤波器阶数的最大取值

n=1:N;

xn=10*sin(2*pi*n*f1+pi/3)+5*sin(2*pi*n*f2+pi/4);

wn=randn(1,N);%产生高斯白噪声

xn=xn+wn;%信号加噪声

subplot(2,1,1);

plot(n,xn);

title('信号加噪声');

%Burg算法

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

eb(1,i)=xn(i);%ef(0,n)=eb(0,n)=x(n)因matlab中数组索引从1开始

end

sum=0;

for i=1:N

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

end %计算R(0)

r(1)=sum/N;%因matlab中数组索引从1开始

%Burg递推

suma=0;sumb=0;

for m=2:(P+1);%循环的阶次，因matlab中数组索引从1开始

for n=m:N

suma=suma+ef(m-1,n)*eb(m-1,n-1);%初步构造KM中的分子

sumb=sumb+ef(m-1,n)^2+eb(m-1,n-1)^2;%构造KM中的分母end

k(m-1)=-2*suma/sumb;%完全构造KM,根据前后向预测均方误差之和最小来求反射系数Km

suma=0;sumb=0;%每次循环前赋零值

for n=2:N

ef(m,n)=ef(m-1,n)+k(m-1)*eb(m-1,n-1);%前向误差

eb(m,n)=eb(m-1,n-1)+k(m-1)*ef(m-1,n);%后向误差

end

end

aa(1,1)=k(1);p(1)=r(1)*(1-k(1)^2);%赋初值

for m=2:P

for i=1:m-1

aa(m,i)=aa(m-1,i)+k(m)*aa(m-1,m-i);%利用levison算法求出的参数

end

aa(m,m)=k(m);

p(m)=p(m-1)*(1-k(m)^2);

end

%由H(z)用freqz()求解功率谱

G=p(P)^(1/2);

B=[G,zeros(1,P)];

A=[1,aa(P,:)];

[H,w]=freqz(B,A,2000);

Hf=abs(H);

subplot(2,1,2);

plot(w/(2*pi),Hf);

title('Burg算法功率谱');

[pks,locs]=findpeaks(Hf,'minpeakheight',10);

f=(locs-1)/(2*2000)

5 MUSIC算法

N=128;n=1:N;

x=10*sin(2*pi*0.2*n+pi/3)+5*sin(2*pi*0.213*n+pi/4)+randn(1,N);

%生成Toeplitz自相关矩阵矩阵

R0=0;sum=0;

for n=1:N

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

end

R0=sum/N;%求R(0)

for m=1:N

temp=0;

for n=1:N-m

temp=temp+x(n)*x(n+m);

end

R(m)=temp/N;

end %求自相关函数

Rx=zeros(N,N);

for m=1:N

for n=1:m-1

Rx(m,n)=R(m-n);%将值赋值给矩阵

end

end

Rx=Rx+Rx';%topelite阵的性质

for m=1:N

Rx(m,m)=R0;%将主对角线的值变为R(0)

end

[V,D]=eig(Rx); %特征分解V是长度为N的列向量,返回特征值的对角矩阵D 和矩阵V，其列是对应的右特征向量，使得A*V = V*D

p=0;%特侦知的个数

for i=1:N,

if D(i,i)/D(N,N)>0.05

p=p+1;

end

end %求p

p

w=[0:2000-1]/1000*pi;

e=exp(-1i*w'*[0:N-1]); %e矩阵

ev=e*V(:,1:N-p); %e*V

Pw=1./real(diag(ev*ev')'); %功率谱

plot(w/(2*pi),Pw);

title('MUSIC算法求功率谱')

xlabel('频率')

axis([0 0.5 0 10])

[pks,locs]=findpeaks(Pw,'minpeakheight',0);

f=(locs-1)/2000 6 信噪比与均方误差的关系

将求功率谱的各种算法改写成函数以便直接调用，由与采用不同功率谱估计方法时，仅需替换代码中的一部分，关系图谱已附，以下仅提供采用Burg时的代码和子程序代码：b=256;n=1:b;

MSE1=0;MSE2=0;

V=26;s=[1:26];

c=-15;10;

M1=zeros(1,v); %定义1行V列的全零矩阵

M2=zeros(1,v); %

for c=1:v %求30次信噪比和均方误差

for i=1:500 %循环500次求均方误差

x=awgn(x1,c,'measured'); %把信噪比为c的高斯白噪声加到信号x1中

pw1=Q(x);%我将burg算法改写成了Q函数，此处直接调用。

[pks,locs] = findpeaks(pw1,'minpeakheight',300); %峰值函数locs为数组中峰值对应的位置

f1=(locs(1,1)-1)/500; %将峰值对应的横坐标转化为归一化频率

f2=(locs(1,2)-1)/500;

MSE1=MSE1+(0.2-f1)^2/500; %求均方误差

MSE2=MSE2+(0.213-f2)^2/500;

end

M1(1,s)=MSE1;

M2(1,s)=MSE2;

MSE1=0;

MSE2=0;

S(1,s)=c;

c=c+1;

end

subplot(1,2,1);

plot(S,M1);

xlabel('SNR/dB');

ylabel('MSE');

title('music算法f1的SDR&MS的关系');

subplot(1,2,2);

plot(S,M2);

xlabel('SNR/dB');

ylabel('MSE');

title(' music算法f1的SDR&MS的关系');

Q函数（burg子程序）的代码：

function y=Q(x)

N=512; %定义采样点数

f1=0.2;f2=0.213;

P=40;%滤波器阶数的最大取值

n=1:N;

x1=10*sin(2*pi*n*f1+pi/3)+5*sin(2*pi*n*f2+pi/4);

x=awgn(x1,10,'measured'); %把信噪比为c的高斯白噪声加到信号x1中

%Burg算法

for i=1:N

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

eb(1,i)=x(i);%ef(0,n)=eb(0,n)=x(n)因matlab中数组索引从1开始

end

sum=0;

for i=1:N

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

end %计算R(0)

r(1)=sum/N;%因matlab中数组索引从1开始

%Burg递推

suma=0;sumb=0;

for m=2:(P+1);%循环的阶次，因matlab中数组索引从1开始

for n=m:N

suma=suma+ef(m-1,n)*eb(m-1,n-1);%初步构造KM中的分子

sumb=sumb+ef(m-1,n)^2+eb(m-1,n-1)^2;%构造KM中的分母end

k(m-1)=-2*suma/sumb;%完全构造KM,根据前后向预测均方误差之和最小来求反射系数Km

suma=0;sumb=0;%每次循环前赋零值

for n=2:N

ef(m,n)=ef(m-1,n)+k(m-1)*eb(m-1,n-1);%前向误差

eb(m,n)=eb(m-1,n-1)+k(m-1)*ef(m-1,n);%后向误差

end

end

aa(1,1)=k(1);p(1)=r(1)*(1-k(1)^2);%赋初值

for m=2:P

for i=1:m-1

aa(m,i)=aa(m-1,i)+k(m)*aa(m-1,m-i);%利用levison算法求出的参数

end

aa(m,m)=k(m);

p(m)=p(m-1)*(1-k(m)^2);

end

%由H(z)用freqz()求解功率谱

G=p(P)^(1/2);

B=[G,zeros(1,P)];

A=[1,aa(P,:)];

[H,w]=freqz(B,A,256);

y=abs(H);

%plot(w/(2*pi),y);

%[pks,locs]=findpeaks(y,'minpeakheight',10); %f=(locs-1)/(2*256)

end