MUSIC算法原理分析与实现

素颜马尾好姑娘i 2022-07-14 10:15 204阅读 0赞

很多人私信博主本人关于阵列信号处理的相关问题，顾在此开一个专栏，写一写个人关于阵列信号处理相关算法的一些文章，因一些内容为本人已发表和待发表论文（注：本文内容已收入知网），为了方便，其中有发杂公式的部分会采用截图的形式说明。

**正文**

学模型，如图2.1所示。  
![这里写图片描述][SouthEast]

![这里写图片描述][SouthEast 1]  
![这里写图片描述][SouthEast 2]  
![这里写图片描述][SouthEast 3]  
![这里写图片描述][SouthEast 4]  
![这里写图片描述][SouthEast 5]  
![这里写图片描述][SouthEast 6]  
![这里写图片描述][SouthEast 7]  
![这里写图片描述][SouthEast 8]  
![这里写图片描述][SouthEast 9]  
![这里写图片描述][SouthEast 10]  
![这里写图片描述][SouthEast 11]  
![这里写图片描述][SouthEast 12]  
![这里写图片描述][SouthEast 13]  
![这里写图片描述][SouthEast 14]  
![这里写图片描述][SouthEast 15]

仿真实现代码：

% DOA estimation by MUSIC 
    % Developed by Megustas (DLUT）
    % Blog ：http://blog.csdn.net/megustas_jjc
    clear all
    close all
    derad = pi/180;        % deg -> rad
    radeg = 180/pi;
    twpi = 2*pi;
    kelm = 8;               % number
    dd = 0.5;               % space 
    d=0:dd:(kelm-1)*dd;     % 
    iwave = 3;              % number of DOA
    theta = [10 30 60];     % angle
    snr = 10;               % input SNR (dB)
    n = 500;                 % 
    A=exp(-j*twpi*d'*sin(theta*derad));%%%% direction matrix
    S=randn(iwave,n);
    X=A*S;
    X1=awgn(X,snr,'measured');
    Rxx=X1*X1'/n;
    InvS=inv(Rxx); %%%%
    [EV,D]=eig(Rxx);%%%% 
    EVA=diag(D)';
    [EVA,I]=sort(EVA);
    EVA=fliplr(EVA);
    EV=fliplr(EV(:,I));
    
    % MUSIC
    for iang = 1:361
            angle(iang)=(iang-181)/2;
            phim=derad*angle(iang);
            a=exp(-j*twpi*d*sin(phim)).';
            L=iwave;    
            En=EV(:,L+1:kelm);
            SP(iang)=(a'*a)/(a'*En*En'*a);
    end
    
    % 
    SP=abs(SP);
    SPmax=max(SP);
    SP=10*log10(SP/SPmax);
    h=plot(angle,SP);
    set(h,'Linewidth',2)
    xlabel('angle (degree)')
    ylabel('magnitude (dB)')
    axis([-90 90 -60 0])
    set(gca, 'XTick',[-90:30:90])
    grid on

关于ESPRIT算法及均匀圆阵和多数据集、宽带信号处理等问题会在后续博客讲到。

[SouthEast]: /images/20220714/2d7b4e53f5da4d31ad5cac14cdaef1de.png
[SouthEast 1]: /images/20220714/1151350bfa5b425daa18f0b5a32be14a.png
[SouthEast 2]: /images/20220714/c9d0e0d231204b429b941b151df95778.png
[SouthEast 3]: /images/20220714/c2de497cda844236bf81d72a01760cc6.png
[SouthEast 4]: /images/20220714/cab93f61ed374e86bc0fb365c22f5517.png
[SouthEast 5]: /images/20220714/e2f78e77a2bc40b684089436d89d95a9.png
[SouthEast 6]: /images/20220714/3e655b7dcb64488282e751331fbd8980.png
[SouthEast 7]: /images/20220714/ffbfb4e221b94ef4b740b29fae1a4fd6.png
[SouthEast 8]: /images/20220714/fab2c0a9dbf94d2e9710cf451a2ce257.png
[SouthEast 9]: /images/20220714/79e80f2f44454aaa8cf3443e539bea04.png
[SouthEast 10]: /images/20220714/2e39a094a6ee4e3987efa837577c09f0.png
[SouthEast 11]: /images/20220714/f1f2761d06014bb0ad5059bdf7d98ee2.png
[SouthEast 12]: /images/20220714/198b6cf2a0594d948a883c6988418821.png
[SouthEast 13]: /images/20220714/a059802b549342d3ae9e0db4022a7c3c.png
[SouthEast 14]: /images/20220714/50488767d3474ce88acba68c79dc1d85.png
[SouthEast 15]: /images/20220714/fb2bcd817e3f405eb30043cc2e8128ac.png