clear all;
%-----------------------------------------------------------------
% 參數設置
SNR=20;
silta1=0;
silta2=0.15;
K=1000;
%-----------------------------------------------------------------
N=10;
bi=0.5;
AI=10^(SNR/20);
v = zeros(N,1);
u = zeros(N,1);
n = zeros(N,1);
vy= zeros(N,1);
for i = 1:N
   v(i) = exp(sqrt(-1)*(-(N+1)/2+i)*2*pi*bi*sin(silta1));
   u(i) = exp(sqrt(-1)*(-(N+1)/2+i)*2*pi*bi*sin(silta2));
end 
%----------------------------------------------------------------------
Mon=4;
capon=zeros(K,1);
music=zeros(K,1);
rmusic=zeros(K,1);
for mon=1:Mon
    mon
    %-------------------------------------------------------------------
    Rx=zeros(N);
    for k=1:20
        %_______________________________________
        s=zeros(N);
        for mm=1:10  
            if rand>0.5
                tmp1=1; 
            else
                tmp1=-1;   
            end 
            if rand>0.5
                tmp2=1;
            else
                tmp2=-1;
            end
            s1=AI*(tmp1+j*tmp2)/(sqrt(2));
            if rand>0.5
                tmp1=1;
            else
                tmp1=-1;
            end
            if rand>0.5
                tmp2=1;
            else
                tmp2=-1;
            end
            s2=AI*(tmp1+j*tmp2)/(sqrt(2));
            s(1:N,1)=s1*v+s2*u;
            n=(randn(N,1)+sqrt(-1)*randn(N,1))/(sqrt(2)); 
            tmp=s1*v+s2*u+n;
            Rx=(Rx*(k-1)+tmp*tmp')
            
        end
        n=(randn(N,N)+sqrt(-1)*randn(N,N))/(sqrt(2)); 
        x=s+n;
        Rx=((k-1)*Rx+x*x')/k+1*eye(10);
        [V,D]=eigs(Rx,2);
        count=0;
        for y=-0.3:0.1:0.3
            count=count+1;
            for i = 1:N
               vy(i)=exp(sqrt(-1)*(-(N+1)/2+i)*2*pi*bi*sin(y));
            end
            p(count)=abs((vy'*Rx*vy)/(vy'*vy));
            mpdr(count)=abs(vy'*inv(Rx)*vy);
            music(count)=abs(vy'*(eye(10)-V*V')*vy);
            % root music
            z=exp(j*y);
            for wyyy=1:N
               Vz(wyyy)=z^(i-1);
               Vzz(wyyy)=(1/z)^(i-1);
           end
           rmusic=abs(Vzz*(eye(10)-V*V')*Vz');
       end
%        [tmp1,tmp2]=max(p);
%        p(tmp2)=-1000;
%        jiao1(k)=-pi/2+(tmp2-1)*pi/200;
%        [tmp1,tmp2]=max(p);
%        jiao2(k)=-pi/2+(tmp2-1)*pi/200;
             
       [tmp1,tmp2]=min(mpdr);
       mpdr(tmp2)=10000;
       jiao3(k)=-pi/2+(tmp2-1)*pi/200;
%        [tmp1,tmp2]=min(mpdr);
%        jiao4(k)=-pi/2+(tmp2-1)*pi/200;
       
       [tmp1,tmp2]=min(music);
       music(tmp2)=10000;
       jiao5(k)=-pi/2+(tmp2-1)*pi/200;
%        [tmp1,tmp2]=min(music);
%        jiao6(k)=-pi/2+(tmp2-1)*pi/200;
       
       [tmp7,tmp8]=min(rmusic);
       rmusic(tmp2)=10000;
       jiao7(k)=-pi/2+(tmp2-1)*pi/200;
%        [tmp1,tmp2]=min(rmusic);
%        jiao8(k)=-pi/2+(tmp2-1)*pi/200;
       
       %------------------------
       tmp_capon=min(abs(jiao3(k)),abs(jiao3(k)-0.15));
       tmp_musicr=min(abs(jiao5(k)),abs(jiao5(k)-0.15));
       tmp_rmusicr=min(abs(jiao7(k)),abs(jiao7(k)-0.15));
       
       capon(k)=capon(k)+tmp_capon;
       musicr(k)=musicr(k)+tmp_musicr;
       rmusicr(k)=rmusicr(k)+tmp_rmusicr;
         
   end
end
capon=capon/K;
musicr=musicr/K;
rmusicr=rmusicr/K;

capon=log10(capon);
musicr=log10(musicr);
rmusicr=log10(rmusicr);


% figure(1)
% plot(1:length(jiao1),jiao1,'r');
% grid on;
% hold on;
% plot(1:length(jiao2),jiao2,'g');
% jiao1(length(jiao1))
% jiao2(length(jiao2))
% hold off;
% 
% figure(2);
% plot(1:length(jiao3),jiao3,'r');
% grid on;
% hold on;
% plot(1:length(jiao4),jiao4,'g');
% jiao3(length(jiao3))
% jiao4(length(jiao4))
% hold off;
% 
% figure(3);
% plot(1:length(jiao5),jiao5,'r');
% grid on;
% hold on;
% plot(1:length(jiao6),jiao6,'g');
% jiao5(length(jiao5))
% jiao6(length(jiao6))
% hold off;
% 
% figure(4);
% plot(1:length(jiao7),jiao7,'r');
% grid on;
% hold on;
% plot(1:length(jiao8),jiao8,'g');
% jiao7(length(jiao7))
% jiao8(length(jiao8))
% hold off;
figure(1);
plot(1:length(capon),capon,'r');
hold on;
plot(1:length(musicr),musicr,'g');
plot(1:length(rmusicr),rmusicr,'b');
grid on;
hold off;

figure(2);
step=1:K;
lgs=log10(step);
semilogx(10.^lgs,capon,'r')
hold on;
semilogx(10.^lgs,musicr,'g');
grid on;
semilogx(10.^lgs,rmusicr,'b');