现代谱估计实验报告

杭州电子科技大学现代谱估计实验报告

现代谱估计实验报告

吴迪松 20040089 信息与信息处理

一、实验题目：

用MATLAB 编写一下算法的程序：

(1)修正协方差法

(2)多重信号分类（MUSIC ）算法

(3)ESPRIT 算法

(4)皮萨论科（Pisarenko ）谐波分解法

并对各算法进行分析。

二、实验原理

(1) 修正协方差法

1、得出输入信号x(n),此输入信号是与白噪声相加的。

xx R ?()()()x

n x n n t =+ 2、求自相关估计器(),xx c j k

()()11**01,()(2()N N P xx n p n c j k x n j x n k x n j x n k N p ???==)()??=??+++?????

∑∑ 3、求滤波器的系数 ?()a

k

1?()(,)(,1)p xx xx k a

k c l k c l ==?∑ 4、求白噪声方差的估计值2?σ

2

1

??(1,1)()(1,)p xx xx k c a k c σ==+∑k 5、求功率谱

()p f 2221?()?1()p j fk k p f a

k e πσ

?==+∑

(2)多重信号分类（MUSIC ）算法

1、估计样本自相关矩阵 ?xx r

2、对作特征值分解 ?xx r

3、确定的最小特征值分解的数目,求出这个最 ?xx r

E n E n 小特征值1,,p M λλ+"，令

()2121p p M E

n σλλλ++=+++" 并求出相对应的特征向量,构造噪声矩阵

1,,p v v +"M

21,,p M V v v +??=??"

4、构造函数：

()211

?MU M H i i p P f e v =+=∑

(3)ESPRIT 算法

1、利用已知信息求 {}(0),(1),()xx xx xx r r r M "

2、由{}(0),(1),()xx xx xx r r r M "构造自相关矩阵xx R 和互相关矩

阵xy R

3、计算xx R 的特征值分解。对于M>P ，最小特征值为

噪声方差2σ

4、利用2σ计算2xx xx C R I σ=?和2xy xy C R I σ=?

5、计算矩阵对{},xx xy C C 的广义特征值分解。从而确定谐

波频率。

(4)皮萨论科（Pisarenko ）谐波分解法

1、求x(n)的自相关函数

xx r 2、求的toeplitz 矩阵xx r xx R

3、求出xx R 的特征值，由此得出最小特征值

4、求出最小特征值对应的一列特征向量

5、求这个特征矢量形成的多项式的根。由这些根求出

它的角度，并且最后求出频率

三、各算法的MATLAB 程序

%AR model parameter estimate

clear

m=sqrt(-1);

delta=0.1;

a1=0.5;

sample=32; %number of sample spot

p=10; %number of sample spot in coef method;

f1=0.05; f2=0.20; f3=0.45; %三个频率分量

fstep=0.01;

fstart=-0.5;

fend=0.5;

f=fstart:fstep:fend; %频率的范围

nfft=(fend-fstart)/fstep+1;

%un=urn+juin

urn= normrnd(0,delta/2,1,sample);

uin= normrnd(0,delta/2,1,sample);

un=urn+m*uin;

%calculate zn 噪声

for n=1:sample-1

zn(1)=un(1);

zn(n+1)=a1*zn(n)+un(n+1);

end

%calculate xn 出入信号

for n=1:sample

xn(n)=2*cos(2*pi*f1*(n-1))+2*cos(2*pi*f2*(n-1))+2*cos(2*pi*f3*(n-1)) +sqrt(2)*real(zn(n));

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%calculate cxx(j,k) 公式3.102

for j=1:p

for k=1:p

s=0;

for n=p+1:sample

s=s+1/2/(sample-p)*(conj(xn(n-j))*xn(n-k)+xn(n-p+j)*conj(xn(n-p+k)));

end

cxx(j,k)=s;

end

end

%calculate cxx0(j,1) 公式 3.103

for j=1:p

s=0;

for n=p+1:sample

s=s+1/2/(sample-p)*(conj(xn(n-j))*xn(n)+xn(n-p+j)*conj(xn(n-p))); end

cxx0(j,1)=s;

end

%calculate a 公式3.104

a=-inv(cxx)*cxx0;

%a(k)*exp(-j2pifk)累加

for i=1:length(f)

sum=0;

for k=1:p

sum=sum+a(k)*exp(-m*2*pi*f(i)*k);

end

%a(k)cxx(1,k)累加

sun=0;

for k=1:p

sun=sun+a(k)*cxx(1,k);

end

delta1=cxx(1,1)+sun;

%AR谱计的形式

XIE(i)=delta1/(abs(1+sum))^2;

end

figure

semilogy(f,XIE);

title('修正协方差算法'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%xn= xn(:); x=xn';

%calculate the value of rxx

for k=0:1:sample-1

s=0;

for n=1:sample-k

s=s+conj(x(1,n))*x(1,n+k);

end

rxx(1,k+1)=(1/sample)*s;

end

%calculate the value of Rx

Rx=zeros(sample,sample);

Rx=toeplitz(rxx(1,1:32));

[U,S,V]=svd(Rx); %特征值分解

Pmusicf=zeros(1,1/fstep+1);

ei=zeros(1,sample);

for i=1:length(f)

for j=1:sample

ei(1,j)=exp(-2*pi*(j-1)*f(i)*m); %calculate the value of ei end;

sum=0;

for k=p+1:sample

sum=sum+abs(ei*V(:,k))^2; %累加e*vi的绝对值

end

Pmusicf(1,i)=(1/sum);

end

figure

semilogy(f,Pmusicf);

title('music算法'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pp=6; %6 ge fen liang

x=xn';

M=length(x);

rxx=xcorr(x,'biased'); % x de zi xiang guan

rxx=[rxx(M:end),0]; %qu hou 32 wei

R=toeplitz(rxx,rxx); %qiu Rxx,Rxy (公式5.158 5.159)

Rxx=R(1:M,1:M);

Rxy=R(1:M,2:end);

D=eig(Rxx); %Rxx de te zheng zhi

[delta2 i]=min(D); %qiu zhi xiao te zheng zhi

for i=pp+1:M

ey=eye(7,7);

Z=ey(1:pp,2:end); %qiu Z

Cxx=Rxx(1:pp,1:pp)-delta2*ey(1:pp,1:pp); %p182 (4)

Cxy=Rxy(1:pp,1:pp)-delta2*Z;

D1=eig(Cxx,Cxy); %p182 (5)

y=angle(D1)/2/pi;

end

figure

stem(y,ones(1,length(y)));

title('esprit算法'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xt=xn';

sin_num=3; %frequency fen liang

N=length(xt);

rxx=xcorr(xt,'biased');%x de zi xiang guang

rxx=rxx(N:(2*sin_num)+N);%rxx de qu zhi

%Frequencies estimation

Rxx=toeplitz(rxx); %qiu Rxx

ev=eig(Rxx); %qiu Rxx de te zheng zhi

[S i]=min(ev); %qiu zhui xiao te zheng zhi

[V D]=eig(Rxx); %te zheng zhi fen jie

a=V(:,i); %qiu zui xiao te zheng zhi dui ying de yi lie te zheng xiang liang

rts=roots(a); %qiu gen

w_est=[];

for i=1:2*sin_num

w_est(i)= angle(rts(i)); %qiu jiao du

end

F=(w_est/(2*pi))'; %qiu pin lv

figure

stem(F,ones(1,length(F)),'*');

title('Pisarenko算法');

四、各算法MATLAB仿真后的分析

四种谱估计算法中修正协方差法是对协方差法的一种改进，它可以得到高分辨率的统计稳定的谱估计值，但不能保证建立一个稳定的

全极点滤波器。而皮萨论科谐波分解法、MUSIC法都是噪声子空间频率估计，都是根据托布列兹自相关矩阵的噪声子空间特征矢量与信号矢量正交这一性质来求解的。ESPRIT法是借助旋转不变技术估计参数，是一种主要特征分解法。从上面四个图中可以明显的看出，它们对频率的估计都较精确。