% function     mssmethod           
% antenna_number                                天線陣元個數
% subarray_length                               子陣陣元數
% subarray_number                               子陣個數
% A0A                                           信號的方位角/度
% A0A_radian                                    設置信號的角度為弧度。
% signal_number                                 信號源個數
% signal_type                                   選擇信號源種類
% signal_frequency                              信號的中心頻率
% sample_frequency                              抽樣頻率
% lambda_bochang                                信號源波長
% antenna_space                                 天線陣元的間距
% samplepoint_number                            對天線陣輸出X的采樣個數
% SNR                                           信號與噪聲的對數信噪比
clc;                                          % MATLAB語言實現MUSIC算法在天線測向上的應用
clear;
close all;
antenna_number=8;                               %天線陣元個數
subarray_length=5;                              %子陣陣元數
subarray_number=antenna_number-subarray_length+1;    %子陣個數
A0A=[-18    -10   0   10];                      %信號的方位角/度
A0A_radian=A0A*pi/180;                          %設置信號的角度為弧度。
signal_number=length(A0A);                      %信號源個數
signal_type='g';
signal_frequency=6*1e9;
sample_frequency=2.5*signal_frequency;
lambda_bochang=3*1e11/signal_frequency;    
antenna_space=0.5*lambda_bochang;               %天線陣元的間距
samplepoint_number=1000;                        %采樣個數
SNR=[7   7   5    5];                           %信號與噪聲的信噪比
%#####設置信號的方向矩陣azimuth_matrice
for row_1=1:antenna_number;                 
    for column_1=1:signal_number;
        t=(row_1-1)*2*pi*antenna_space/lambda_bochang*sin(A0A_radian(column_1));
        azimuth_matrice(row_1,column_1)=exp(-j*t);
    end;
end;
if   signal_type=='r'
        %調用產生實余弦信號源和高斯噪聲的子程序signal_noise_real
        [signal,noise]=signal_noise_real(antenna_number,signal_number,samplepoint_number,SNR,signal_frequency,sample_frequency);
elseif   signal_type=='c'
        %調用產生復正弦信號源和高斯噪聲的子程序signal_noise_complex
        [signal,noise]=signal_noise_complex(antenna_number,signal_number,samplepoint_number,SNR,signal_frequency,sample_frequency);
else
        %調用產生復高斯信號源和高斯噪聲的子程序signal_noise_gauss
        [signal,noise]=signal_noise_gauss(antenna_number,signal_number,samplepoint_number,SNR,signal_frequency,sample_frequency);
end;     
        %調用求信號源（噪聲）相關系數矩陣的子程序corr_coefficient    
        [signal_corr,noise_corr]=corr_coefficient(signal,noise,signal_number,antenna_number);      
         co_12=(mean(conj(signal(1,:)).*signal(2,:))-mean(signal(1,:))*mean(signal(2,:)))/sqrt(var(signal(1,:))*var(signal(2,:)));

%#####求天線陣輸出信號X、及其協方差矩陣R
         X=azimuth_matrice*signal+noise;
         Rf=zeros(subarray_length,subarray_length);
         J=zeros(subarray_length,subarray_length);
                 for  row_3=1:subarray_length
                      J(row_3,subarray_length-(row_3-1))=1;
                 end; 
         for   subarray_out_index=1:subarray_number
               Xf_1{subarray_out_index}=X(subarray_out_index:subarray_out_index+subarray_length-1,:);
               Rf_1{subarray_out_index}=1/samplepoint_number*Xf_1{subarray_out_index}*Xf_1{subarray_out_index}';
               Rf=Rf+Rf_1{subarray_out_index};
          end;
         Rf=Rf/subarray_number;                               %前向平滑相關矩陣/SSf
         Rb=J*conj(Rf)*J;                                     %后向平滑相關矩陣/SSb
         Rmss=Rf+Rb;                                          %修正空間平滑法相關矩陣 
         R=1/samplepoint_number*X*X';                         %MUSIC相關矩陣
%調用求空間譜函數Pmu、信號方位角signal_azimuth_solved及對應譜值Pmu_solved的子程序fromeigencomposition_lastsolved        
        [Pmu_music,search_azimuth,AOA_solved_music, Pmu_solved_music, signal_number_juded_music,big_index_music]=fromeigencomposition_lastsolved(R,antenna_number,samplepoint_number,antenna_space,lambda_bochang);
        [Pmu_ss,   search_azimuth,AOA_solved_ss,    Pmu_solved_ss,    signal_number_juded_ss,   big_index_ss]=fromeigencomposition_lastsolved(Rf,subarray_length,samplepoint_number,antenna_space,lambda_bochang);
        [Pmu_mss,  search_azimuth,AOA_solved_mss,   Pmu_solved_mss,   signal_number_juded_mss,  big_index_mss]=fromeigencomposition_lastsolved(Rmss,subarray_length,samplepoint_number,antenna_space,lambda_bochang);
         AOA_number_solved_music=length(AOA_solved_music);
         AOA_number_solved_ss=length(AOA_solved_ss);
         AOA_number_solved_mss=length(AOA_solved_ss);
%#####畫出MUSIC、MMUSIC兩種方法下的空間譜函數曲線         
         figure(1);
         plot(search_azimuth,Pmu_mss,'r-.',search_azimuth,Pmu_ss,'g:',search_azimuth,Pmu_music,'b');
         xlabel('入射角\theta / °'),ylabel('空間譜函數 /dB');grid;legend('MUSIC','SS','MSS');
         hold on;
         plot(AOA_solved_mss,Pmu_solved_mss,'m*',AOA_solved_ss,Pmu_solved_ss,'b*',AOA_solved_music,Pmu_solved_music,'g*');
              figure(2)
              subplot(3,1,1);  plot(search_azimuth,Pmu_mss,AOA_solved_mss,Pmu_solved_mss,'r*');  
              grid;            legend('MUSIC');             title('空間譜函數 /dB');
              x_label=num2str(AOA_solved_mss');
              y_label=num2str(Pmu_solved_mss');
               text(AOA_solved_mss,Pmu_solved_mss,['('*ones(AOA_number_solved_mss,1),x_label,ones(AOA_number_solved_mss,1)*' , ',y_label,ones(AOA_number_solved_mss,1)*')']);      
              subplot(3,1,2);  plot(search_azimuth,Pmu_ss,AOA_solved_ss,Pmu_solved_ss,'r*');  
              grid;            legend('SS');             
              x_label=num2str(AOA_solved_ss');
              y_label=num2str(Pmu_solved_ss');
               text(AOA_solved_ss,Pmu_solved_ss,[ones(AOA_number_solved_ss,1)*'  (',x_label,ones(AOA_number_solved_ss,1)*' , ',y_label,ones(AOA_number_solved_ss,1)*')']);      
               subplot(3,1,3);  plot(search_azimuth,Pmu_music,AOA_solved_music,Pmu_solved_music,'r*');  
              grid;            legend('MSS');
              x_label=num2str(AOA_solved_music');
              y_label=num2str(Pmu_solved_music');
               text(AOA_solved_music,Pmu_solved_music,[ones(AOA_number_solved_music,1)*'  (',x_label,ones(AOA_number_solved_music,1)*' , ',y_label,ones(AOA_number_solved_music,1)*')']);      
               xlabel('入射角\theta / °');