MUSIC算法的改进

采用共轭重排和传播算子对基本MUSIC算法做出改进，使得算法能分辨角度较近的相干信源。

tic
%-------将数据共轭重排，后进行正交运算。
clc;
clear;
close all;
f=4*10^9 ;       %1 .8750e+009;% 信源载波频率
c=3e+8;
w=c/f;% 信源波长
d=w*0.4;% 
md=10;% 阵元数
th=[10 20];% 来波方向角
nd=length （th）;
snr=0;%:2:30;% 信噪比
snr1=10^（snr/20）;
kp=200;
% b2c=0.1;
rad=pi/180;

%--------------------------------------------------------------------------
for m=1:md
    for n=1:nd
        a （mn）=exp （-j*2*pi*d*（m-1）*sin （th （n）*rad）/w）;
    end
end
%--------------------------------------------------------------------------

 for p=1:nd
    for q=1:kp
        p1=rand （11）;
        p2=rand （11）; % 生成两个随机数 （在0-1之间）
        sr （pq）=sqrt （-2*snr1*snr1*log （p1））*cos （2*pi*p2）;    % 实部
        si （pq）=sqrt （-2*snr1*snr1*log （p1））*sin （2*pi*p2）;    % 虚部
        s （pq）=sr （pq）+j*si （pq）;                         % 模拟生入射信号矩阵S （t）
    end
end
%-------------------------------------------------------------------------
n=wgn （mdkp0‘complex‘）;% 生成复信号高斯白噪声，md*kp阶，功率0dBW即功率为1，幅度为1
%--------------------------------------------------------------------------
x=a*s+n;% 生成模拟入射信号矩阵
%--------------------------------------------------------------------------
% X1=x（1:nd1:md）;
% X2=x（nd+1:end1:md）;
%%%%%%共轭PM
nm=eye（mdmd）;
jiaohuan=rot90（nm）;
ll=conj（x）;
bn=jiaohuan*ll;
z=[xbn];%为md*2kp维的
r=z*z‘/2*kp;% 构成协方差矩阵为md维的

% K=bn（1:nd1:md）;
% L=bn（nd+1:end1:md）;
% rr=size（X1*X1‘+K*K‘）
% tt=size（X2*X1‘+L*K‘）
% G=[X1*X1‘+K*K‘;X2*X1‘+L*K‘];
% H=[X1*X2‘+K*L‘;X2*X2‘+L*L‘];
% P=inv（G‘*G）*G‘*H;
% Q=[P‘-diag（ones（1md-nd））];
G=r（:1:nd）;
H=r（:nd+1:end）;
P=inv（G‘*G）*G‘*H;
Q=[P‘-diag（ones（1md-nd））];
% D=（Q*Q‘）^（-0.5）;
F=Q‘*（（Q*Q‘）^（-0.5））;
theta=-90:0.1:90;
for i=1:length （theta）
    AA=exp （-j*2*pi*d*（0:md-1）‘*sin （theta （i）*rad）/w）;
    WW=AA‘*F*F‘*AA;
    pu（i）=abs（1/WW）;

end
pu =10*log （pu ）;
theta=linspace （-90901801）;
figure（1）
plot（thetapu‘k‘‘linewidth‘2）;% thetapuu‘b‘
grid
xlabel （‘角度/\[Degree]‘） 
ylabel （‘空间谱/dB‘）;
title （‘PM算法均匀线阵谱‘）;

grid on

%-----------------------------------------
%-----------------------------------------
n=0;% 扫描出谱峰值并显示出估计角
for i=2:1799
    while （pu （i）>pu （i-1））&&（pu （i）>pu （i+1））&&（pu （i）>-16）
        n=n+1;
        dd （n）=pu （i）;
        angle （n）=（i-901）/10  % 将采样点数转化为度，一度为10个步长
        break
    end
end

%-----在谱峰输入搜索数值----------------------------------------------------
[maxdkd]=max （dd）;
for i=1:n
    doa=num2str （angle （i））;
    if i        text （angle （i）dd （i）doa）;% 在横坐标angle，纵坐标dd显示数值doa
    else
        text （angle （i）dd （i）doa）;
    end
end
toc

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2757  2012-06-08 11:03  gong-e-pm\gongepm.asv

     文件       2775  2012-06-08 08:28  gong-e-pm\gongepm.m

     文件       2893  2012-06-01 20:20  gong-e-pm\gongepmshuju.asv

     文件       2893  2012-06-01 20:21  gong-e-pm\gongepmshuju.m

     目录          0  2014-06-29 09:55  gong-e-pm

----------- ---------  ---------- -----  ----

                11318                    5