rootmusic.m

Root-MUSIC 算法进行了深入的研 究

%root-music
clear
w1=pi/4;
w2=-pi/3;
w3=pi/6;
w4=-pi/5;
kp=2048;
a1=[exp(-j*pi*sin(w1)*[0:9])].';
a2=[exp(-j*pi*sin(w2)*[0:9])].';
a3=[exp(-j*pi*sin(w3)*[0:9])].';
a4=[exp(-j*pi*sin(w4)*[0:9])].';
A=[a1,a2,a3,a4];
k=1:kp;
s=[1.3*cos(k*0.015);sin(k*0.05);cos(k*0.02);sin(k*0.035)];
n=randn(10,kp)+j*randn(10,kp);
x=A*s+n;
r=x*(x)'/kp;
[V,D]=eig(r);
[b,index]=sort(diag(D));
Vn=V(:,index);
Un=Vn(:,[1:6]);
II=Un*(Un)';
z =sym('z','unreal');
for i=1:10
pz(i)=z^(i-1);
pzt(i)=z^(10-i);
end
p=sym2poly(pzt*II*(pz).');
zz=roots(p);
zz1=abs(abs(zz)-ones(18,1));
[zb,ix]=sort(zz1);
zzn=zz(ix);
zm=zzn(1:8);
ww=angle(zm);
theta=asin(-ww/pi);
t=theta*180/pi;
stem(t,ones(1,8));grid;
axis([-90 90 0 2]);
title('root-music算法');
ylabel('P(\theta) (dB)');
xlabel('\theta (deg)');
[Wn Pn]=ginput;
hold on
plot(Wn,Pn,'or');
for k=1:length(Wn)
    text(Wn(k),Pn(k),['  \theta=',num2str(Wn(k)),'P=',num2str(Pn(k))]);
end
hold off