%8元均勻線陣，6個入射信號分別以W1，W2，W3，W4，W5，W6的角度入射，利用MUSIC（多重
%信號分類法）進行波達方向估計
%function y = MUSIC(u)
N = 2000;%采樣數
%6個入射波方向
w1 = 10*pi/180;
w2 = 50*pi/180;
w3 = 30*pi/180;
w4 = 60*pi/180;
w5 = 45*pi/180;
%w6 = 80*pi/180;
%陣列導向矢量
A = [1                   1                    1                    1                    1                    ;%1 
    exp(-j*pi*sin(w1))   exp(-j*pi*sin(w2))   exp(-j*pi*sin(w3))   exp(-j*pi*sin(w4))   exp(-j*pi*sin(w5))   ;%exp(-j*pi*sin(w6))
    exp(-j*2*pi*sin(w1)) exp(-j*2*pi*sin(w2)) exp(-j*2*pi*sin(w3)) exp(-j*2*pi*sin(w4)) exp(-j*2*pi*sin(w5)) ;%exp(-j*2*pi*sin(w6)) 
    exp(-j*3*pi*sin(w1)) exp(-j*3*pi*sin(w2)) exp(-j*3*pi*sin(w3)) exp(-j*3*pi*sin(w4)) exp(-j*3*pi*sin(w5)) ;%exp(-j*3*pi*sin(w6)) 
    exp(-j*4*pi*sin(w1)) exp(-j*4*pi*sin(w2)) exp(-j*4*pi*sin(w3)) exp(-j*4*pi*sin(w4)) exp(-j*4*pi*sin(w5)) ;%exp(-j*4*pi*sin(w6)) 
    exp(-j*5*pi*sin(w1)) exp(-j*5*pi*sin(w2)) exp(-j*5*pi*sin(w3)) exp(-j*5*pi*sin(w4)) exp(-j*5*pi*sin(w5)) ;%exp(-j*5*pi*sin(w6))
    exp(-j*6*pi*sin(w1)) exp(-j*6*pi*sin(w2)) exp(-j*6*pi*sin(w3)) exp(-j*6*pi*sin(w4)) exp(-j*6*pi*sin(w5)) ;%exp(-j*6*pi*sin(w6))
    exp(-j*7*pi*sin(w1)) exp(-j*7*pi*sin(w2)) exp(-j*7*pi*sin(w3)) exp(-j*7*pi*sin(w4)) exp(-j*7*pi*sin(w5)) ;%exp(-j*7*pi*sin(w6)) 
    ];
f = 9e+8;%載頻
t = linspace(0,1,N);%采樣點
s = [sin(10*pi*t).*cos(2*pi*f*t);
    2*cos(100*pi*t).*cos(2*pi*f*t);
    8*exp(t).*cos(2*pi*f*t);
    100*cos(2*pi*t).*cos(2*pi*f*t);
    25*sin(60*pi*t).*cos(2*pi*f*t);
    %10*cos(45*pi*t).*cos(2*pi*f*t)
    ];%信號
n = [randn(1,N);randn(1,N);randn(1,N);randn(1,N);randn(1,N);randn(1,N);randn(1,N);randn(1,N)].*0.0001;
sum = zeros(8,8);
for i = 1:N
x(:,:,i) = A*s(:,i);%+n(:,i);
sum = sum+x(:,:,i)*x(:,:,i)';
end
R = sum./N;
[q,r] = eig(R);
q = [q(:,1) q(:,2) q(:,3)];
w = linspace(0,90,2000);
u = w*pi/180;
for m = 1:2000
a = [1;exp(-i*pi*sin(u(m)));
    exp(-i*2*pi*sin(u(m)));exp(-i*3*pi*sin(u(m)));
    exp(-j*4*pi*sin(u(m)));exp(-j*5*pi*sin(u(m)));
    exp(-j*6*pi*sin(u(m)));exp(-j*7*pi*sin(u(m)))];
f(m) = 1/(a'*q*q'*a);
cap(m) =1/(a'*inv(R)*a);
end
plot(w,abs(f));
hold on;
plot(w,abs(cap),'r-');