modifiedmusic.m

利用MATLAB实现修正MUSIC算法的DOA估计的仿真

%修正的Music算法,前向空间平滑技术
%多源，未知源个数,信号源相干
%产生信号源
clear;
clc;
N=7;%信号源个数为N
%我们这样设计信号源，频率分别为f0,2f0,…，(N-1)f0,相位为随机初始相位
f0=1000;
fs=10000;%采样频率
t=0:1/fs:0.1-1/fs;
L=length(t);
for ii=1:5
    s(ii,:)=sin(2*pi*ii*f0*t+rand*2*pi);
end
s(6,:)=s(3,:);
s(7,:)=s(4,:);
%仍用均匀阵列接收
M=15;
fc=1000000;
c=3e8;
lamda=c/fc;
d=lamda/2;
%阵列流型矢量
% theta=[-70 -50 -30 -10 20 40 55]/180*pi;
theta=[-60 -50 -40 -20 0 20 40]/180*pi;
for ii=1:N
    for jj=1:M
        A(jj,ii)=exp(-j*2*pi/lamda*(jj-1)*d*sin(theta(ii)));
    end
end
varn=0.1;
x=A*s+varn*randn(size(A*s));
%========================================================================
%利用修正的Music算法
%求自相关函数Rxx
Rxx=0;
for ii=1:L
    Rxx=Rxx+x(:,ii)*x(:,ii)';
end
Rxx=Rxx/L;
%定义一个交换矩阵，M*M维
J=0;
for ii=1:M
    J(ii,M+1-ii)=1;
end
J
Ryy=J*conj(Rxx)*J;
R=(Rxx+Ryy)/2;
%特征值分解
[V D]=eig(R);
D
% %噪声子空间,这里我们假定我们已知信源个数N
% U=V(:,1:M-N);
%如果未已知信源个数，则需要通过特征值来判断
%建立一种规则，D中的特征值是从小到大排列，我们把数量级最小的特征值认为是噪声空间对应的特征值
%其余的特征值认为是信号空间对应的特征值
for ii=1:M
    Lam(ii)=real(D(ii,ii));
end
for ii=1:M
    for jj=-30:30
        if(Lam(ii)/10^(jj)>=1 & Lam(ii)/10^(jj)<10)
            Num(ii)=jj;
        end
%         if(Num(ii)==-30 | Num(ii)==30)
%             disp('error!特征值不在范围内');
%         end
    end
end
Num
N_N=1;
for ii=2:M
    if(Num(ii)<Num(1)+2)
        N_N=N_N+1;
    end
end
U=V(:,1:N_N);
the=[-90:90]/180*pi;
L_the=length(the);
for ii=1:L_the
    for jj=1:M
        a(jj,1)=exp(-j*2*pi/lamda*(jj-1)*d*sin(the(ii)));
    end
    P(ii)=(a'*a)/(a'*U*U'*a);
end
P=10*log10(real(P));
plot(the/pi*180,P);
grid on;
hold on;
xlabel('波达方向');
ylabel('空间谱');
grid on;