%改進(jìn)雙向空間平滑music算法  7陣元等距陣波達(dá)方向估計（d=λ/2）,2信源， 4個子陣（4陣元） 最小范數(shù)
clear all
clc
d=1;
lma=2;
N=1024;                                               %設(shè)定采樣點(diǎn)數(shù)
SNR=[2 4 6];                         
fs=8;                             
f=[0.9 0.9 0.5];                                     %設(shè)兩個信號元的頻率，單位為MHz
w=[20 16 40]*pi/180;                                 %入射方向角向量
phi=[60 20 30]*pi/180;                               %兩入射信號初始相位
M=9;                                                 %設(shè)定陣元數(shù)
AMP=10.^(SNR./20);              
m=0:M-1;
n=0:N-1;
S=repmat(AMP',1,N).*exp(-j*(2*pi*f.'/fs*n+phi'*ones(1,N)));          %產(chǎn)生兩個遠(yuǎn)場正弦信號
A=exp(-j*2*pi*d*m'*(sin(w)/lma));                       %方向矩陣
X=A*S;
noise=randn(M,N)+j*randn(M,N);
X=X+noise;                                            %接受數(shù)據(jù)矩陣
r=X*(X)'/N;
J=rot90(eye(5));
for k=1:5
    b=sprintf('r%d',k);
    c=sprintf('rb%d',k);
    rk=sprintf('rk%d',k);
    F=[zeros(5,k-1) eye(5) zeros(5,5-k)];
    f=F.';
    eval([b,'=','F','*r*','f;']);
    eval([c,'=J*conj(',b,')*J;']);
    eval([rk,'=(',b,'+',c,')/2;']);
end
rx=(rk1+rk2+rk3+rk4+rk5)/5;
[V,D]=eig(rx);
[bn,index]=sort(diag(D));
Vn=V(:,index);
Un=Vn(:,[1:2]);
II=Un*(Un)';
w=[-pi/2:pi/180:pi/2];
for l=1:length(w)
    ao=[exp(-j*pi*sin(w(l))*[0:4])].';
    pm(l)=(ao)'*ao*(inv((ao)'*II*ao));
     p(l)=10*log10(abs(pm(l)));
end
figure(1)
polar(w,p);title('P(f/b)')
figure(2)
t=w*180/pi;
plot(t,p);grid;
title('修正雙向空間平滑算法');
ylabel('P(\theta) (dB)');
xlabel('\theta (deg)');
hold on