mssmethod.m

三种DOA方法：经典的MUSIC,前后向空间平滑和改进的空间平滑算法

                      % 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 / °');