%改進雙向空間平滑music算法，10陣元 等距陣波達方向估計（d=λ/2）,2信源  5個子陣（6陣元）
v=2;
N=2000;                                                %設定采樣點數
SNR=[-5 8];                         
fs=8;                             
f=[0.9 0.9];                                          %設兩個信號元的頻率，單位為MHz
w=[20 30]*pi/180;                                     %入射方向角向量
phi=[60 20]*pi/180;                                   %兩入射信號初始相位
M=10;                                                 %設定陣元數
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)));          %產生兩個遠場正弦信號
A=exp(-j*2*pi*d*m'*(sin(w)/v));                       %方向矩陣
X=A*S;
noise=randn(M,N)+j*randn(M,N);
X=X+noise;                                            %接受數據矩陣
r=X*(X)'/N;
J=rot90(eye(5));
for k=1:6
    b=sprintf('r%d',k);
    c=sprintf('rb%d',k);
    rk=sprintf('rk%d',k);
    F=[zeros(5,k-1) eye(5) zeros(5,6-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:3]);
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