%cyclic MUSIC算法
clear all
hold on
sensor_number=10;
signal_direction=[10 15  25];signal_number=length(signal_direction);
Fc = [10 10  10  17  17  17    17 17 100*rand(1,signal_number)];%載頻
signal_snr=10;
snr=sqrt(2*10^(signal_snr/10));
signal_direction=signal_direction*pi/180;
F0=3e+2;
light_vector=3e+8; %light_velocity

% 產生信號
Fd=5;%x的采樣頻率,有幾個周期
Fs=100;%y的采樣
Ts=1/Fs;
M=4;% 2是BPSK,4是QPSK
num_rand=50;%共產生幾組
snapshot_number=Fs*num_rand;ND=snapshot_number;
for ij=1:signal_number
    serier=floor(rand(1,num_rand)*M);   % 向上取整
    x=ones(Fd,1)*serier;
    x=x(:);
    y=dmod(x,Fc(ij),Fd,Fs,'psk',M);y=y.';%產生 Mpsk信號
    signal_s(ij,:)=snr*y;
end
% signal_s(2,:)=signal_s(1,:);
% signal_s(3,:)=signal_s(1,:);
% signal_s(5,:)=signal_s(1,:);
% signal_s(6,:)=signal_s(1,:);

% 產生數據
ddd=0.5*light_vector/F0;r1=ddd;r2=0;
for i1=1:sensor_number
  ant_x1(i1)=(i1-1)*r1;
  ant_y1(i1)=(i1-1)*r2;
end
ant_x=ant_x1';ant_y=ant_y1';sensor_position=[ant_x,ant_y];
delay=sensor_position*[sin(signal_direction);cos(signal_direction)];
steering_matrix=exp(-1*i*2*pi*delay*F0/light_vector);
MA=sensor_number;ND=snapshot_number;
noise=gennoise([1 -0.2 0.6],[1 0.2 0.7],sensor_number,snapshot_number);
%noise=randn(size(steering_matrix*signal_s)); 
X=steering_matrix*signal_s+noise;

% 產生共軛循環矩陣
Fcj=10;%選擇循環頻率
flagg=1;%循環平穩算法
% flagg=2;%共軛循環平算法
for ij1=1:MA % (xcorr(x1,x2)).*phas(-ND:ND)=xcorr((x1.*phas(0:Nd),x2))
    for ij2=1:MA
        x1=X(ij1,:);x2=X(ij2,:);
        if flagg==1
            xx=xcorr(x1,x2);              %循環平穩算法
        elseif flagg==2
            xx=xcorr(x1,conj(x2));        %共軛循環平算法
        end
        phas=exp(i*2*pi*Fcj*(-ND+1:ND-1)*Ts);
        R(ij1,ij2)=mean(xx.*phas);
    end
end
% R=X*X';
R=R*R';ss=MA;
% RR=R;
R=inverdiag(MA)*conj(R)*inverdiag(MA)+R;kk=3;ss=MA-kk+1;RR=zeros(ss);for ij=1:kk   RR=RR+R(ij:ij+ss-1,ij:ij+ss-1);end

% MUSIC算法估計
NX=signal_number;
[U value]=eig(RR);
value=diag(value);
[value,index]=sort(value); 
UN=U(:,index(1:ss-NX));
ddoa=-90:0.5:90;
delay=sensor_position*[sin(ddoa*pi/180);cos(ddoa*pi/180)];
steering=exp(-1*i*pi*delay*(2*F0+Fcj)/light_vector);
steering=steering(1:ss,:);
spectrum1=steering'*UN*UN'*steering;
spectrum2=diag(spectrum1)+eps;
spe1=steering'*steering;
spe2=diag(spe1);
clear spectrum1 spe1
spectrum=spe2./spectrum2;
spectrum=abs(spectrum);
plot(ddoa,10*log10(spectrum))
if flagg==1
    title('循環平穩算法')
elseif flagg==2
    title('共軛循環平穩算法');
end