精华内容
下载资源
问答
  • 原理探究——空间平滑MUSIC算法
    千次阅读 热门讨论
    2021-05-24 15:15:26

      空间平滑MUSIC算法的核心在于:信号相干性与协方差矩阵的联系;同时与矩阵的秩、行列向量间线性相关/无关的联系

      MUSIC算法中阵列接收信号的表达式为:X=AS+N,A为空间阵列流型矩阵,即天线阵列的数学表达;S为空间信号矢量,即信号的数学表达;N为噪声。当信号不相干时,求得X的协方差矩阵为:
    协方差矩阵

      此时发现Rs为满秩矩阵,且由于噪声N的存在,Rx也为满秩矩阵。这是因为协方差矩阵Rs中,每一个元素都在利用协方差公式计算不同信号间的相关程度!而当这些信号都是不相干的情况下,协方差矩阵Rs必然满秩。从而使得有噪声的情况下Rx也为满秩矩阵。而在证明了Rx矩阵的共轭转置等于它本身后, Rx同时也为厄米特矩阵

      由于厄米特矩阵有两个性质:1、所有特征值都为实数;2、不同特征值对应的特征向量是正交的。利用这些性质,可将Rx矩阵进行特征分解,通过特征值的属性来对特征向量进行筛选分类,划分出信号子空间与噪声子空间。最终利用两个子空间的正交性构造出谱函数,即可在天线接收的角度范围内进行谱峰搜索,找出各信号的来向角。

      讲到这里,再去看空间平滑算法处理相干信号的思路就很清晰透彻了。由于信号之间的相干性,如果直接计算阵列接收信号X(t)的协方差矩阵,会使得公式中的信号源矩阵Rs不再是满秩矩阵。而将Rx矩阵特征分解后划分的噪声子空间,也由于Rs矩阵出现的秩亏损,导致噪声子空间中混入了信号源对应

    更多相关内容
  • 空间平滑MUSIC算法

    2020-02-10 21:30:25
    空间平滑MUSIC算法 有声压阵和矢量阵的
  • 空间平滑music算法

    2014-05-22 10:50:06
    空间平滑music算法,用于估计相干信号入射角度
  • 面阵中二维角度估计 Unitary -ESPRIT算法MATLAB程序
  • 空间平滑MUSIC算法的MATLAB程序
  • 基于空间平滑music源程序,包括前向平滑,后向平滑和双向平滑
  • 基于经典music算法及其修正算法的DOA估计以及性能仿真
  • 这个算法我自己用过了好多遍,是可以实现的,很实用的,下载下来后就可以直接粘贴上去的 空间平滑MUSIC算法MATLAB程序
  • 空间平滑算法处理相干源, 两种实现方案
  • 空间平滑MUSIC算法,空间平滑music算法原理,matlab源码
  • 基于空间平滑MUSIC算法的相干信号DOA估计(1)

    千次阅读 多人点赞 2020-10-29 20:38:20
    空间平滑MUSIC算法(1) 1. 前言 在上一篇博客中有提到,当多个入射信号相干时,传统MUSIC算法的效果就会不理想。具体原因是多个入射信号相干时,有部分能量就会散发到噪声子空间,使得MUSIC算法不能对其进行有效...

    空间平滑MUSIC算法(1)

    1. 前言

    在上一篇博客中有提到,当多个入射信号相干时,传统MUSIC算法的效果就会不理想。具体原因是多个入射信号相干时,有部分能量就会散发到噪声子空间,使得MUSIC算法不能对其进行有效估计。
    针对这种情况,解相干方法被提出。本文主要讲解通过降维处理解相干,之所以叫降维处理是因为这种方法将原阵列拆分成很多个子阵列,通过子阵列的协方差矩阵重构接收数据协方差矩阵,阵列的自由度DOF会随之降低,即可分辨的相干信号数降低。
    先看看传统MUSIC算法对相干信号进行DOA估计的效果。

    MATLAB代码

    % fss.m
    % Code For Music Algorithm
    % The Signals Are All Coherent
    % Author:痒羊羊
    % Date:2020/10/29
    
    clc; clear all; close all;
    %% -------------------------initialization-------------------------
    f = 500;                                        % frequency
    c = 1500;                                       % speed sound
    lambda = c/f;                                   % wavelength
    d = lambda/2;                                   % array element spacing
    M = 20;                                         % number of array elements
    N = 100;                                        % number of snapshot
    K = 6;                                          % number of sources
    coef = [1; exp(1i*pi/6);... 
            exp(1i*pi/3); exp(1i*pi/2);... 
            exp(2i*pi/3); exp(1i*2*pi)];            % coherence coefficient, K*1
    doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals
    
    %% generate signal
    dd = (0:M-1)'*d;                                % distance between array elements and reference element
    A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
    S = sqrt(2)\(randn(1,N)+1i*randn(1,N));         % vector of random signal, 1*N
    X = A*(coef*S);                                 % received data without noise, M*N
    X = awgn(X,10,'measured');                      % received data with SNR 10dB
    
    %% calculate the covariance matrix of received data and do eigenvalue decomposition
    Rxx = X*X'/N;                                   % covariance matrix
    [U,V] = eig(Rxx);                               % eigenvalue decomposition
    V = diag(V);                                    % vectorize eigenvalue matrix
    [V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
    U = U(:,idx);                                   % reset the eigenvector
    P = sum(V);                                     % power of received data
    P_cum = cumsum(V);                              % cumsum of V
    
    %% define the noise space
    J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
    J = J(1);                                       % number of principal component
    Un = U(:,J+1:end);
    
    %% music for doa; seek the peek
    theta = -90:0.1:90;                             % steer theta
    doa_a = exp(-1i*2*pi*dd*sind(theta)/lambda);    % manifold array for seeking peak
    music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
    music = 10*log10(music/max(music));             % normalize the result and convert it to dB
    
    %% plot
    figure;
    plot(theta, music, 'linewidth', 2);
    title('Music Algorithm For Doa', 'fontsize', 16);
    xlabel('Theta(°)', 'fontsize', 16);
    ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
    grid on;
    

    在这里插入图片描述
    可以看到,对于相干信号,传统MUSIC算法DOA估计效果很差。

    2. 空间平滑算法

    降维处理解相干方法主要包括空间平滑类处理算法,而空间平滑类处理算法可以分成前向空间平滑算法(FSS)、后向平滑算法(BSS)、前后向平滑算法(FBSS),正如上面所说,这些算法的估计效果很好,但是损失了阵列孔径,导致可分辨的相干信号数降低。

    2.1 线性阵列信号模型

    线阵模型
    关于这个线性阵列中符号的具体说明可以参考上一篇博文。子空间方法——MUSIC算法

    2.2 前向空间平滑算法

    前向空间平滑算法将阵列分成多个互相重叠的子阵列,然后对子阵接收数据的协方差矩阵作平均。当子阵阵元数目大于等于相干信号数时,可以有效解相干。
    前向空间平滑
    如上图所示,我们把M元阵列均匀分成L个子阵列,每个子阵列有N=M-L+1个阵元。以最左边的子阵列为参考阵列,定义第J个子阵的接收数据为:
    X J f = [ x J , x J + 1 , ⋯   , x J + N − 1 ] X_{J}^{f}=\left[x_{J}, x_{J+1}, \cdots, x_{J+N-1}\right] XJf=[xJ,xJ+1,,xJ+N1]那么第J个子阵接收数据的协方差矩阵(也称作空间平滑矩阵)可以表示为
    R N j j = A 1 D j − 1 R s ( D j − 1 ) H A 1 H + σ 2 I N j = 1 , 2 , ⋯   , L R_{N}^{j j}=A_{1} D^{j-1} R_{s}\left(D^{j-1}\right)^{H} A_{1}^{H}+\sigma^{2} I_{N} \quad j=1,2, \cdots, L RNjj=A1Dj1Rs(Dj1)HA1H+σ2INj=1,2,,L其中, D = diag ⁡ [ e − j − 2 π d λ sin ⁡ θ 1 , e − j 2 π d λ sin ⁡ θ 2 , ⋯ e − j 2 π d λ sin ⁡ θ R ] D=\operatorname{diag}\left[\mathrm{e}^{-j\frac{-2 \pi d}{\lambda} \sin \theta_{1}}, \mathrm{e}^{-j\frac{2 \pi d}{\lambda} \sin \theta_{2}}, \cdots \mathrm{e}^{-j\frac{2 \pi d}{\lambda} \sin \theta_{R}}\right] D=diag[ejλ2πdsinθ1,ejλ2πdsinθ2,ejλ2πdsinθR],A1为第一个子阵即参考阵列的流型矩阵。
    因此,前向空间平滑后的协方差矩阵可以由各个子阵的协方差矩阵求平均得到。
    R ~ f = 1 L ∑ j = 1 L R N j j \tilde{R}^{f}=\frac{1}{L} \sum_{j=1}^{L} R_{N}^{j j} R~f=L1j=1LRNjj利用前向空间平滑协方差矩阵和MUSIC算法就可以分辨出多个相干信号的方位。可以证明,该方法最大可以检测M/2个相干的信号。

    MATLAB代码

    % Code For Music Algorithm
    % The Signals Are All Coherent
    % Author:痒羊羊
    % Date:2020/10/29
    
    clc; clear all; close all;
    %% -------------------------initialization-------------------------
    f = 500;                                        % frequency
    c = 1500;                                       % speed sound
    lambda = c/f;                                   % wavelength
    d = lambda/2;                                   % array element spacing
    M = 20;                                         % number of array elements
    N = 100;                                        % number of snapshot
    K = 6;                                          % number of sources
    L = 10;                                         % number of subarray
    L_N = M-L+1;                                    % number of array elements in each subarray
    coef = [1; exp(1i*pi/6);... 
            exp(1i*pi/3); exp(1i*pi/2);... 
            exp(2i*pi/3); exp(1i*2*pi)];            % coherence coefficient, K*1
    doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals
    
    %% generate signal
    dd = (0:M-1)'*d;                                % distance between array elements and reference element
    A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
    S = sqrt(2)\(randn(1,N)+1i*randn(1,N));         % vector of random signal, 1*N
    X = A*(coef*S);                                 % received data without noise, M*N
    X = awgn(X,10,'measured');                      % received data with SNR 10dB
    
    %% reconstruct convariance matrix
    %% calculate the covariance matrix of received data and do eigenvalue decomposition
    Rxx = X*X'/N;                                   % origin covariance matrix
    Rf = zeros(L_N, L_N);                           % reconstructed covariance matrix
    for i = 1:L
        Rf = Rf+Rxx(i:i+L_N-1,i:i+L_N-1);
    end
    Rf = Rf/L;
    [U,V] = eig(Rf);                                % eigenvalue decomposition
    V = diag(V);                                    % vectorize eigenvalue matrix
    [V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
    U = U(:,idx);                                   % reset the eigenvector
    P = sum(V);                                     % power of received data
    P_cum = cumsum(V);                              % cumsum of V
    
    %% define the noise space
    J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
    J = J(1);                                       % number of principal component
    Un = U(:,J+1:end);
    
    %% music for doa; seek the peek
    dd1 = (0:L_N-1)'*d;
    theta = -90:0.1:90;                             % steer theta
    doa_a = exp(-1i*2*pi*dd1*sind(theta)/lambda);   % manifold array for seeking peak
    music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
    music = 10*log10(music/max(music));             % normalize the result and convert it to dB
    
    %% plot
    figure;
    plot(theta, music, 'linewidth', 2);
    title('Music Algorithm For Doa', 'fontsize', 16);
    xlabel('Theta(°)', 'fontsize', 16);
    ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
    grid on;
    

    前向空间平滑MUSIC
    可以看到,在6个入射信号均相干的情况下,基于前向平滑的MUSIC算法能较好地对其进行DOA估计,但仍存在估计精度的问题,比如真实入射角度为75°的信号的方位被估计成了74.2°。
    后向空间平滑MUSIC算法和前/后向空间平滑MUSIC算法在下一篇博客中讲解。
    基于空间平滑MUSIC算法的相干信号DOA估计(2)

    参考论文:《基于空间平滑MUSIC算法的相干信号DOA估计,陈文锋,吴桂清》
    代码Code均为原创。
    欢迎转载,表明出处。

    展开全文
  • 空间平滑MUSIC算法(只适用于均匀线性阵列) 1、由于信号阵列会接收到不同方向上的相干信号,而相干信号会导致信源协方差矩阵的秩亏缺(不是满秩矩阵),从而使得信号特征向量发散到噪声子空间去(就是组成噪声子...

    空间平滑MUSIC算法(只适用于均匀线性阵列)

    1、由于信号阵列会接收到不同方向上的相干信号,而相干信号会导致信源协方差矩阵的秩亏缺(不是满秩矩阵),从而使得信号特征向量发散到噪声子空间去(就是组成噪声子空间的向量中混入了信号源对应的特征向量)。

    2、对于相干信号波达方向估计:就是想办法把信号协方差矩阵的秩恢复到等于信号源的个数(这样我们得到的信号子空间的秩才是和信号源的个数相等的)。

    3、其中一种思路:在进行空间阵列谱函数估计出波达方向之前,先进行预处理,将协方差的秩恢复到信号源数(即去相关)

    4、去相关中的一种思路:降维处理,比如空间平滑MUSIC算法,牺牲有效阵列孔径来实现信号源的去相干(因为毕竟子阵列的孔径比原阵列的孔径小)

                                            

    clear all;close all;clc
    %阵列的基础信息
    derad=pi/180;%角度转换成弧度
    redeg=180/pi;%弧度转换成角度
    twpi=2*pi;
    Melm=7;%阵元数目
    K=6;%阵元数目(后面会从这两种数目中选择一个)
    dd=0.5;%阵元间距
    d=0:dd:(Melm-1)*dd;%阵元总长度(7-1)*0.5=3,这里共7个元素!
    iwave=3;%信号源数目(有用信号,不包括噪声源)
    theta=[0 30 60];%这三个信号源的入射角度分别是0、30、60
    n=200;%采样点数目(即快拍数)
    A=exp(-j*twpi*d.'*sin(theta*derad));%构造方向矢量(阵列流型),注意d.’
     
    %构造相干信号源(了解相干信号的定义)
    S0=randn(iwave-1,n);%信号源符合高斯分布(正态分布),矩阵维度自行设定,主要考虑后面运算不要起冲突就行
    S=[S0(1,:);S0];%S0是2*n维,S0(1,:)是1*n维,而;是下一行,所以构成了相干信号源S是3*n维
    X0=A*S;%阵列接收信号(与信号源和方向向量有关)
    snr=10;%信噪比,也可以设定为好几个值,snr=[10 20 30],但要保证后面矩阵运算维数不起冲突
    X=awgn(X0,snr,'measured');%往接收信号中添加噪声,模拟真实环境(存在噪声)接收到的信号
    Rxxm=X*X'/n;%协方差矩阵(这只是待定的,最后要用的是Rxx)
    issp=1;%选用的平滑算法模式,我们设计了几种,从中选用一种,比如issp=1是第一种模式,issp=2是第二种
     
    %设计了几种空间平滑算法,以子程序中自定义函数的形式进行调用(ssp和issp)
    if issp == 1%注意,等于是==,而赋值是=
      Rxx = ssp(Rxxm,K);%自定义函数ssp的输出,经过在ssp函数中进行空间平滑算法,输出协方差矩阵
    elseif issp == 2%选择第二种空间平滑算法,注意是在另一个自定义函数issp中运算
      Rxx = issp(Rxxm,K);%将参数Rxxm和kelm代入自定义函数中,最后输出Rxx,比以前那些算法的要好
    else
    Rxx=Rxxm;%没有选用设计的空间平滑算法,所以沿用以前方法的协方差矩阵
    K=Melm;%因为沿用以前方法,所以阵元数目用7
    end
     
    %对协方差矩阵进行特征值分解
    [EV,D]=eig(Rxx);%分解后得到特征值组成的向量形式D和特征值对应特征向量组成的矩阵EV
    DA=diag(D)';%构造对角阵DA,特征值向量D的元素是对角线上的数值,并将矩阵转置(下下句有用)
    [DA,I]=sort(DA);%按照特征值大小对对角阵元素升序排列,原先的顺序为I
    DA=fliplr(DA);%翻转矩阵(中心对称旋转),所以对角阵元素此时是按照降序排列的(结合diag(D)’
    EV=fliplr(EV(:,I));%先使EV中列向量按照I的顺序排列,在进行翻转
     
    %阵列进行扫描
    for jiaodu =1:361%从1循环到361(而且循环数必须为正,>0)
    angle(jiaodu )=(jiaodu -181)/2;%使得阵列扫描角度是从-90到90
    phim=derad*angle(jiaodu );%提前精简,避免下面表达式过于臃肿
    a=exp(-j*twpi*d(1:K)*sin(phim)).';%方向向量,其中取d的前kelm个元素,不会取到Melm个元素。并且要转置.主要是后面矩阵运算,维度有要求
    L=iwave;%信源数目
    En=EV(:,L+1:K);%噪声空间,前L列为信号空间,而后面的则为噪声空间
    P(jiaodu )=(a'*a)/(a'*En*En'*a);%空间阵列谱函数估计
    end
    P=abs(P);%求绝对值或者复数实部与虚部的平方和的算术平方根
    Pmax=max(P);
    P=10*log10(P/Pmax);%归一化,结合信噪比的相关公式
    figure%创建一个用来显示图形输出的一个窗口对象
    h=plot(angle,P);
    set(h,'Linewidth',2)%将所有线宽设置为2
    xlabel('angle(degree)')
    ylabel('magnitude(dB)')
    axis([-90 90 -60 0])%设置坐标轴范围
    set(gca,'XTick',[-90:30:90],'YTick',[-60:10:0])%设置网格的显示格式,gca获取当前figure的句柄,选择x、y轴的要进行标注的刻度
    grid on 
    hold on
    legend('平滑MUSIC')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %子程序1
    %自定义函数ssp
    function Rxx=ssp(Rxxm,K)%自定义函数(参数1待定协方差矩阵,参数2阵元数目)
    [M,MM]=size(Rxxm);%返回Rxxm矩阵的维度大小
    N=M-K+1;%子阵列的数目
    Rxx=zeros(K,K);%初始化矩阵Rxx,类似于初始化为0
    for i=1:N%第1个子阵列到第N个子阵列
    Rxx=Rxx+Rxxm(i:i+K-1,i:i+K-1);%向前滑动,注意每个子阵列是k个阵元,不是k-1个
    end
    Rxx=Rxx/N;%子阵列协方差矩阵可以相加后平均取代原来意义上的协方差矩阵
    %子程序2
    %自定义函数issp
    function Rxx=issp(Rxxm,K)
    [M,MM]=size(Rxxm);
    N=M-K+1;
    J=fliplr(eye(M));%翻转对角阵,中心对称,维数为M*M
    Rxxb=(Rxxm+J*Rxxm.'*J)/2;%前后向平滑协方差矩阵(同时结合了Toeplitz的协方差矩阵估计值);共轭倒序得到的信号子空间与原来的一样
    Rxx=zeros(K,K);
    for i=1:N
    Rxx=Rxx+Rxxb(i:i+K-1,i:i+K-1);
    end
    Rxx=Rxx/N;

    将MUSIC算法的步骤归纳为:

    1、根据n个接收信号矢量,得到协方差矩阵的估计值Rxxm=X*X’/n

    2、若各子阵列的阵列流形相同(等距线阵满足这个假设),则子阵列协方差矩阵可以相加后平均取代原来意义上的协方差矩阵Rxx=Rxx/N

    3、对平均后得到的协方差矩阵进行特征值分解[EV,D]=eig(Rxx)

    4、按照特征值的大小顺序,把与信号个数kelm相等的最大特征值对应的特征向量看作信号子空间,剩下的(M-K)个特征值对应特征向量看作噪声子空间

    5、使θ变化,即扫描角度从-90到90,按照空间阵列谱函数公式来计算,通过寻求峰值来得到波达方向的估计值

     

    展开全文
  • 基于空间平滑MUSIC算法的相干信号DOA估计(2)

    千次阅读 多人点赞 2020-10-31 11:06:09
    空间平滑MUSIC算法(2) 继续上一篇博客,继续讲后向空间平滑和前/后向空间平滑MUSIC算法。 基于空间平滑MUSIC算法的相干信号DOA估计(1) 2.3 后向空间平滑算法 后向空间平滑更准确的说是共轭后向空间平滑,它是对后...

    空间平滑MUSIC算法(2)

    继续上一篇博客,继续讲后向空间平滑和前/后向空间平滑MUSIC算法。
    基于空间平滑MUSIC算法的相干信号DOA估计(1)

    2.3 后向空间平滑算法

    后向空间平滑
    后向空间平滑更准确的说是共轭后向空间平滑,它是对后向子阵列地共轭接收数据协方差矩阵进行平滑。定义第一个共轭后向子阵列由 { M , M − 1 , ⋯   , M − p + 1 } \{M, M-1, \cdots, M-p+1\} {M,M1,,Mp+1} 组成, 第二个子阵列由 { M − 1 , M − 2 , ⋯   , M − p } \{M-1, M-2, \cdots, M-p\} {M1,M2,,Mp} 组成,依次组成的子阵列个数为 L = M − p + 1 个 \mathrm{L}=\mathrm{M}-\mathrm{p}+1 个 L=Mp+1
    容易知道,共轭后向空间平滑协方差矩阵 R ~ b \tilde{R}^{b} R~b与前向空间平滑协方差矩阵 R ~ f \tilde{R}^{f} R~f的关系:
    R ~ b = J ∗ ( R ~ f ) ∗ ∗ J \tilde{R}^{b}=J*(\tilde{R}^{f})^**J R~b=J(R~f)J利用后向空间平滑协方差矩阵和MUSIC算法也可以分辨出多个相干信号的方位。可以证明,该方法最大也可以检测M/2个相干的信号。

    MATLAB代码

    % bss.m
    % Code For Music Algorithm Based On Backward Spatial Spectrum
    % The Signals Are All Coherent
    % Author:痒羊羊
    % Date:2020/10/29
    
    clc; clear all; close all;
    %% -------------------------initialization-------------------------
    f = 500;                                        % frequency
    c = 1500;                                       % speed sound
    lambda = c/f;                                   % wavelength
    d = lambda/2;                                   % array element spacing
    M = 20;                                         % number of array elements
    N = 100;                                        % number of snapshot
    K = 6;                                          % number of sources
    L = 10;                                         % number of subarray
    L_N = M-L+1;                                    % number of array elements in each subarray
    coef = [1; exp(1i*pi/6);... 
            exp(1i*pi/3); exp(1i*pi/2);... 
            exp(2i*pi/3); exp(1i*2*pi)];            % coherence coefficient, K*1
    doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals
    
    %% generate signal
    dd = (0:M-1)'*d;                                % distance between array elements and reference element
    A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
    S = sqrt(2)\(randn(1,N)+1i*randn(1,N));         % vector of random signal, 1*N
    X = A*(coef*S);                                 % received data without noise, M*N
    X = awgn(X,10,'measured');                      % received data with SNR 10dB
    
    %% reconstruct convariance matrix
    %% calculate the covariance matrix of received data and do eigenvalue decomposition
    Rxx = X*X'/N;                                   % origin covariance matrix
    H = fliplr(eye(M));                             % transpose matrix
    Rxxb = H*(conj(Rxx))*H;
    Rf = zeros(L_N, L_N);                           % reconstructed covariance matrix
    for i = 1:L
        Rf = Rf+Rxxb(i:i+L_N-1,i:i+L_N-1);
    end
    Rf = Rf/L;
    [U,V] = eig(Rf);                                % eigenvalue decomposition
    V = diag(V);                                    % vectorize eigenvalue matrix
    [V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
    U = U(:,idx);                                   % reset the eigenvector
    P = sum(V);                                     % power of received data
    P_cum = cumsum(V);                              % cumsum of V
    
    %% define the noise space
    J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
    J = J(1);                                       % number of principal component
    Un = U(:,J+1:end);
    
    %% music for doa; seek the peek
    dd1 = (0:L_N-1)'*d;
    theta = -90:0.1:90;                             % steer theta
    doa_a = exp(-1i*2*pi*dd1*sind(theta)/lambda);   % manifold array for seeking peak
    music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
    music = 10*log10(music/max(music));             % normalize the result and convert it to dB
    
    %% plot
    figure;
    plot(theta, music, 'linewidth', 2);
    title('Music Algorithm For Doa', 'fontsize', 16);
    xlabel('Theta(°)', 'fontsize', 16);
    ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
    grid on;
    

    后向空间平滑可以看到,在6个入射信号均相干的情况下,基于后向空间平滑的MUSIC算法能较好地对其进行DOA估计,且估计精度较高。

    2.3 前/后向空间平滑算法

    把前向和共轭后向空间平滑协方差矩阵定义为前向空间平滑协方差矩阵和共轭后向空间平滑协方差矩阵的平均值,即:
    R ‾ = 1 2 ( R ~ f + R ~ b ) \overline{\mathbf{R}}=\frac{1}{2}\left(\tilde\mathbf{R}^{f}+\tilde\mathbf{R}^{b}\right) R=21(R~f+R~b)那么只要空间平滑的次数大于等于相干信号源的个数,前向和共轭后向空间平滑协方差矩阵一般情况下都是满秩的。使用前/后向空间平滑的方法最多可以检测的相干信号源个数为 2 M / 3 2 \mathrm{M} / 3 2M/3 个。你可能很好奇这个最大相干信号源检测数是怎么得到的?
    假设: 阵列天线的阵元数为 M \mathrm M M个,前/后向空间平滑次数分别为 L \mathrm L L次,则每个子阵的阵元数为 N = M − L + 1 \mathrm N={\mathrm M}-\mathrm{L}+1 N=ML+1 个, 同时可以知道, 可以分辨的最大信号个数为 M − L \mathrm{M}-\mathrm{L} ML个, 也就是子阵的阵元个数减 1 ; 1 ; 1;前后向分别平滑 N \mathrm{N} N次可以解相干信号的个数为 2 L 2 \mathrm{L} 2L个, 最大情况下,两者相等所以 M − L = 2 L \mathrm M- \mathrm L=2 \mathrm L ML=2L,也就是 L = M / 3 ; \mathrm{L}=\mathrm{M} / 3 ; L=M/3; 所以 2 L = 2 M / 3 , 2 \mathrm{L}=2 \mathrm{M} / 3, 2L=2M/3, 所以前/后向空间平滑最大可以解相干的信号个数为 2 M / 3 2 \mathrm{M} / 3 2M/3 个。所以说采用前/后向空间平滑的改进技术可以很大地提高阵列孔径。

    MATLAB代码

    % fbss.m
    % Code For Music Algorithm Based On Forward And Backward Spatial Spectrum
    % The Signals Are All Coherent
    % Author:痒羊羊
    % Date:2020/10/29
    
    clc; clear all; close all;
    %% -------------------------initialization-------------------------
    f = 500;                                        % frequency
    c = 1500;                                       % speed sound
    lambda = c/f;                                   % wavelength
    d = lambda/2;                                   % array element spacing
    M = 20;                                         % number of array elements
    N = 100;                                        % number of snapshot
    K = 6;                                          % number of sources
    L = 10;                                         % number of subarray
    L_N = M-L+1;                                    % number of array elements in each subarray
    coef = [1; exp(1i*pi/6);... 
            exp(1i*pi/3); exp(1i*pi/2);... 
            exp(2i*pi/3); exp(1i*2*pi)];            % coherence coefficient, K*1
    doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals
    
    %% generate signal
    dd = (0:M-1)'*d;                                % distance between array elements and reference element
    A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
    S = sqrt(2)\(randn(1,N)+1i*randn(1,N));         % vector of random signal, 1*N
    X = A*(coef*S);                                 % received data without noise, M*N
    X = awgn(X,10,'measured');                      % received data with SNR 10dB
    
    %% reconstruct convariance matrix
    %% calculate the covariance matrix of received data and do eigenvalue decomposition
    Rxx = X*X'/N;                                   % origin covariance matrix
    H = fliplr(eye(M));                             % transpose matrix
    Rxxb = H*(conj(Rxx))*H;
    Rxxfb = (Rxx+Rxxb)/2;
    Rf = zeros(L_N, L_N);                           % reconstructed covariance matrix
    for i = 1:L
        Rf = Rf+Rxxfb(i:i+L_N-1,i:i+L_N-1);
    end
    Rf = Rf/L;
    [U,V] = eig(Rf);                                % eigenvalue decomposition
    V = diag(V);                                    % vectorize eigenvalue matrix
    [V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
    U = U(:,idx);                                   % reset the eigenvector
    P = sum(V);                                     % power of received data
    P_cum = cumsum(V);                              % cumsum of V
    
    %% define the noise space
    J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
    J = J(1);                                       % number of principal component
    Un = U(:,J+1:end);
    
    %% music for doa; seek the peek
    dd1 = (0:L_N-1)'*d;
    theta = -90:0.1:90;                             % steer theta
    doa_a = exp(-1i*2*pi*dd1*sind(theta)/lambda);   % manifold array for seeking peak
    music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
    music = 10*log10(music/max(music));             % normalize the result and convert it to dB
    
    %% plot
    figure;
    plot(theta, music, 'linewidth', 2);
    title('Music Algorithm For Doa', 'fontsize', 16);
    xlabel('Theta(°)', 'fontsize', 16);
    ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
    grid on;
    

    前/后向空间平滑因为前/后向空间平滑的改进技术很大地提高了阵列孔径,从以上的DOA结果图可以看出分辨率提高了。

    参考论文:《阵列信号处理相关技术研究,陈四根》
    代码Code均为原创。
    欢迎转载,表明出处。

    展开全文
  • 主要进行波束空间MUSIC算法的matlab仿真
  • 基于互素数组的时空空间平滑MUSIC算法进行DOA估计
  • 基本MUSIC和前后平滑MUSIC算法,matlab 代码,可处理相干与非相干信号,完整版本,可运行
  • 前向平滑MUSIC算法.m

    2020-06-22 08:41:32
    利用MATLAB实现了前向平滑MUSIC算法,通过仿真结果可以看出,该算法在低信噪比的情况下仍具有良好的分辨能力
  • 详细的MUSIC空间平滑解相干算法,程序能够运行,好使,能够根据程序改参数进行分析。
  • 比较详细的介绍了空间谱估计基础和DOA估计模型,研究了DOA估计中的MUSIC算法,给出了MUSIC算法的原理和步骤,并通过一些计算机仿真实验,得出了MUSIC算法的性能分析。
  • 阵列信号处理 MUSIC以及空间平滑处理,前、后、前后向平滑算法,经典的MUSIC算法不能对相干信号进行解算,经过空间平滑处理后,可以对相干信号进行解算,并进行波束形成,判断波达方向。
  • 基于特征空间 MUSIC算法的相干信号波达方向空间平滑估计
  • 空间平滑算法详解

    2022-03-10 10:06:46
      MUSIC算法能有效运行的前提是矩阵RSR_SRS​是非奇异的,即各条传播路径不相干。如果L路到达信号中存在Q路相干信号(Q≤L),则通过MUSIC算法能被检测到的信号数量为L-Q+1,能被解出的信号数量为L-Q。   ...
  • 双向平滑MUSIC算法.m

    2020-06-22 08:47:00
    利用MATLAB实现了双向平滑MUSIC算法,通过仿真结果可以看出,该算法在低信噪比的情况下仍具有良好的分辨能力
  • 利用空间平滑方法进行相干源DOA估计,采用MUSIC算法实现
  • 后向平滑MUSIC算法.m

    2020-06-22 08:44:29
    利用MATLAB实现了后向平滑MUSIC算法,通过仿真结果可以看出,该算法在低信噪比的情况下仍具有良好的分辨能力
  • 参考论文:《相干信源波达方向估计的加权空间平滑算法》王布宏 王永良 陈辉武汉空军雷达学院重点实验室 武汉430010摘要 提出了一种用于空间相干源 DOA 估计的加权空间平滑算法 WSS, weighted spatial smoothing常规...
  • 一种基于加权空间平滑的新MUSIC算法,赵娥,陈桂丽,本文针对基于空间平滑的新MUSIC算法对相干源的分辨率较低的问题,提出了一种基于加权空间平滑的新MUSIC算法。该算法充分利用了协方��

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,489
精华内容 595
关键字:

空间平滑music算法

友情链接: combobox.rar