mmusic_music.m

music高阶谱分析算法研究

%                        function     mmusic_music           
% antenna_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;                               %天线阵元个数
A0A=[-20   0   10];                                  %信号的方位角/度
A0A_radian=A0A*pi/180;                          %设置信号的角度为弧度。
signal_number=length(A0A);                      %信号源个数
signal_type='r';
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];                           %信号与噪声的信噪比

%#####设置信号的方向矩阵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;
         Jm=zeros(antenna_number,antenna_number);
                 for  row_3=1:antenna_number
                      Jm(row_3,antenna_number-(row_3-1))=1;
                 end; 
         R=1/samplepoint_number*X*X';                               %MUSIC法相关矩阵R
         Rb_m=Jm*conj(R)*Jm;
         Rmmusic=R+Rb_m;                                            %修正MUSIC法相关矩阵Rmmusic 
%调用求空间谱函数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_mmusic,search_azimuth,AOA_solved_mmusic,Pmu_solved_mmusic,signal_number_juded_mmusic,big_index_mmusic]=fromeigencomposition_lastsolved(Rmmusic,antenna_number,samplepoint_number,antenna_space,lambda_bochang);
        AOA_number_solved_music=length(AOA_solved_music);
        AOA_number_solved_mmusic=length(AOA_solved_mmusic);
%#####画出MUSIC、MMUSIC两种方法下的空间谱函数曲线         
         figure(1);
         plot(search_azimuth,Pmu_mmusic,'r:',search_azimuth,Pmu_music,'b');
         xlabel('入射角\theta / °'),ylabel('空间谱函数 /dB');grid;legend('MMUSIC','MUSIC');
         hold on;
         plot(AOA_solved_mmusic,Pmu_solved_mmusic,'m*',AOA_solved_music,Pmu_solved_music,'g*');
              figure(2)
              subplot(2,1,1);  plot(search_azimuth,Pmu_mmusic,AOA_solved_mmusic,Pmu_solved_mmusic,'r*');  
              grid;            legend('MMUSIC');             title('空间谱函数 /dB');
              x_label=num2str(AOA_solved_mmusic');
              y_label=num2str(Pmu_solved_mmusic');
               text(AOA_solved_mmusic,Pmu_solved_mmusic,[ones(AOA_number_solved_mmusic,1)*'  (',x_label,ones(AOA_number_solved_mmusic,1)*' , ',y_label,ones(AOA_number_solved_mmusic,1)*')']);      
              subplot(2,1,2);  plot(search_azimuth,Pmu_music,AOA_solved_music,Pmu_solved_music,'r*');  
              grid;            legend('MUSIC');
              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 / °');