msc.m

此程序是在MATLAB上实现用于短波定位的MUSIC算法仿真

%该程序完成主要功能:利用music算法对均匀圆阵天线阵进行仿真
%修改时间2008年6月27日



N = 100000;
M = 9;              %天线根数
D = 4;              %信号个数
w0 = 3 * 10^6;      %载频
beamlength=3*10^8/w0;%波长
rate=1.5;            %圆半径与波长之比 
r = rate*beamlength;        %圆阵半径
s0 = 100;           %基频
c=3*10^8;           %光速
xl = 360;
yl = 90;

Rx = zeros(M,M);   %生成一个M*M零矩阵 
caiyangbeishu=4;  %采样频率与载频之比
t = 1:1/(caiyangbeishu*w0):1+(N-1)/(caiyangbeishu*w0);

for i = 1 : N
    X = matrix_X( M, D, r, t(i), w0, s0 );
    
    Rx = Rx + X*X';
end

Rx = Rx/N;        %N组数据的平均协方差矩阵 

[K , sa] = eig(Rx); %利用函数eig()分析方阵的特征值和特征向量,分别放于sa和k中,特征向量是先噪声所对应的特征向量后信号对应的特征向量

K_Nis = K( : , [1:M-D] );%将噪声特征向量取出

k = 1 : M;

for i = 0 : xl-1
    st = 2*pi/xl*i;%将角度转换成弧度
    for j = 0 : yl-1
        fi = pi/2/yl*j;%将角度转换成弧度
        aw = r/c * cos( 2*pi*(k-1)/M - st) * cos( fi );
        aw = exp( aw*w0*1j );%1*M的矩阵
        
        P_msc(i+1,j+1) = 1 / ( aw*(K_Nis*K_Nis')*aw' );
    end
end
P_msc = 10*log10( abs(P_msc) );

x = 1:xl;
y = (1:yl)';

mesh(y,x,P_msc);%画三维图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%以下主要完成找出峰值点工作
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


line = max( max(P_msc) )/3;

num = 0;

for i=1:xl
    for j =1:yl
        if i-2 < 1
            i1 = 1;
        else
            i1 = i-2;
        end
        if i+2 > xl
            i2 = xl;
        else
            i2 = i+2;
        end
        if j-2 < 1
            j1 = 1;
        else
            j1 = j-2;
        end
        if j+2 > yl
            j2 = yl;
        else
            j2 = j+2;
        end
        if P_msc(i,j)>= P_msc(i1,j) & P_msc(i,j)>= P_msc(i1,j1) & P_msc(i,j)>= P_msc(i1,j2) & P_msc(i,j)>= P_msc(i2,j1) & P_msc(i,j)>= P_msc(i2,j2) & P_msc(i,j)>= P_msc(i2,j) & P_msc(i,j)>= P_msc(i,j1) & P_msc(i,j)>= P_msc(i,j2) &  P_msc(i,j) > line
            num = num +1;
            out(num,1) = i/xl*360;
            out(num,2) = j/yl*90;           
        end
    end
end

for i = 1 : num
    for j = i+1 : num
        if ( rem(abs( out(i,1)-out(j,1) ),360) + abs( out(i,2)-out(j,2) ) <= 10 & out(i,1) ~= -1 & out(i,2) ~= -1 ) 
            out(j,1) = -1 ;
            out(j,2) = -1 ;
        end
    end
end

out0 = [];

j = 1;

for i = 1 : num
    if out(i,1) ~= -1 ;
        out0(j,1) = out(i,1);
        out0(j,2) = out(i,2);
        j = j + 1;
    end
end
j-1
out0