basicmusic.m

利用MATLAB实现基本的music算法

%基本Music算法，线阵
clear;
clc;
M=9;%阵元个数
%==================================================================
%产生信源，N=4
N=2;
L=1000;%每个信源的样本点数
Theta=[-40 60];%信源的入射方向
for ii=1:N
    SBasic(ii,:)=randint(1,L);
end
%对基带信号进行调制，FSK方式
for ii=1:N
    SModulate(ii,:)=fskmod(SBasic(ii,:),2,100,2,200);
%     SModulate(ii,:)=pskmod(SBasic(ii,:),2);
end
% YY=SModulate(1,:);
% YY_len=length(YY)
% YY_fft=fft(YY);
% Amp_YY=YY_fft.*conj(YY_fft);
% Amp_YY=fftshift(abs(Amp_YY))
% Amp_YY=10*log10(Amp_YY);
% fs_YY=400;
% freq_YY=-fs_YY/2:fs_YY/YY_len:fs_YY/2-fs_YY/YY_len;
% plot(freq_YY,Amp_YY);
L_SModulate=length(SModulate(1,:));
Fs=4;
Fd=1;
R=0.5;
Delay=5;
for ii=1:N
    S(ii,:)=rcosflt(SModulate(ii,:),Fd,Fs,'normal/Fs',R,Delay); %设置上升余弦滤波器
end
L_S=length(S(1,:));
% %计算S（1）的功率
% Power_S1=S(1,:)*S(1,:)'/L_S;

fs=1000000;%定义载频
c=3e8;
lamda=c/fs;
k=2*pi/lamda;
%阵列流型矢量
d=lamda/2;
for ii=1:N
    for jj=0:M-1
        A(jj+1,ii)=exp(-k*d*sin(Theta(ii)/180*pi)*jj);%阵列流型矢量
    end
end
%阵元上接收的数据为x
varn=0.001;%高斯噪声
N=varn*randn(size(A*S));
x=A*S+N;         %窄带信号数学模型
x
Rxx=0;          %预初始化x（t）的协方差矩阵Rxx里各元素为0
for ii=1:L_S
    Rxx=Rxx+x(:,ii)*x(:,ii)'; %计算Rxx
end
Rxx=Rxx/L_S;
Rxx
[V D]=eig(Rxx);%对RXX进行矩阵特征分解，v为特征向量，Dwei特征值
D
%确定信源数目
%以数量级确定信源的数目，D中的特征值从小到大排列
for ii=1:M
    Lamda_R(ii)=D(ii,ii)  
end 
Lamda_R
for ii=1:M
    for jj=-30:30
        if(Lamda_R(ii)/10^jj>=1 & Lamda_R(ii)/10^jj<10)  %？？？？？？？？  
        Num(ii)=jj;
        end
    end
end
Num
N_S=0;%N_S为信源数目
for ii=2:M
    if(Num(ii)~=Num(1))         %????
        N_S=N_S+1;
    end
end
%求解噪声子空间

U=V(:,1:M-N_S);
cs_U=U*U';
%谱估计
Theta0=-90:90;
Theta0=Theta0/180*pi;
L_Theta0=length(Theta0);
for ii=1:L_Theta0           %对整个线阵的方向角进行功率谱搜索，搜索出功率最大值对应的方向角即为信号源的DOA
    for jj=0:M-1
        a(jj+1)=exp(-j*k*d*sin(Theta0(ii))*jj);    
    end
    P(ii)=a*a'/(a*U*U'*a');
end
P=real(P);
P=10*log10(P);    %log10 为特定的matlab M函数
plot(Theta0/pi*180,P);