基于正交偶极子的四元数MUSIC算法及其Matlab代码

引言
信源数估计
MUSIC算法
基于正交偶极子的MUSIC算法
- 正交偶极子模型
- 正交偶极子的阵列接受模型
- 基于正交偶极子的MUSIC算法
- - 模值约束法求极化信息
基于正交偶极子的四元数MUSIC算法
- 四元数的阵列接受模型
- 四元数MUSIC算法

引言

本文介绍了空间谱估计中的信源数估计、MUSIC算法、正交偶极子的MUSIC算法、极化参数估计的谱峰搜索法、模值约束法和四元数MUSIC算法，并给出相应的代码。

信源数估计

信息论方法的一般的数学形式如下在这里插入图片描述
式中，L(k)是对数概似函数，p(k)是障碍函数。通过对两者的不同的选择可以得到不一样的准则。下面介绍EDC信息论准则，其为 (2-2)
式中，n为被估计的雷达阵列接收信号的信源数(也称为自由度)，L为采样数，其A(n)中为概似函数，且
在这里插入图片描述
式中的C(L)需满足如下条件：
(1)
(2) (2)

当C(L)等于以上情况时，准则EDC便成为了一致性估计。
在EDC信息论准则中选择C(L)分别为1，(lnL)/2及(ln lnL)/2时，就可以得到最小信息、最小描述长度及HQ准则，即
在这里插入图片描述

根据文献，可以得出以下的结论：
（1）最小信息准则不是一致性估计(显而易见，C(L) =1时以上的第2个条件并不成立)，即在大采样数时，它仍然有较大的信号源估计的误差概率；而最小描述长度准则相对较好；HQ准则的估计雷达阵列接收信号的信源数性能在这两者之间，这是由该准则中的障碍函数项导致的。
（2）最小描述长度准则是一致性估计，也可以说在高信噪比时该准则有较好的特性，此时信号源估计的错误概率比最小信息准则小。但是，在小信噪比时该准则相比于最小信息准则有较高的错误概率。
（3）EDC准则中C(L)取下式时：
在这里插入图片描述
EDC准则也就是最小描述长度准则，从某种角度来说，可以说最小描述长度准则是它的一种特殊情况。
（4）EDC准则中C(L)取下式时：

EDC准则等价于HQ准则，从某种角度来说，HQ准则是它的一种特殊情况。在低信噪比的情况下来看，这三种算法中HQ准则性能最好，其次便是最小信息算法。
综合考虑，该文及其对应代码使用MDL准则。信源数估计的代码（完整代码见文末）：

 for s=0:(M-1)
     a=0;
     b=1;
 for m=(s+1):M
    a=a+D2(m);
    b=b*D2(m); 
 end
 aa1=(1/(M-s))*a;
 aa2=b^(1/(M-s));
 ld(s+1)=((1/(M-s))*a)/(b^(1/(M-s)));
 MDL(s+1)=kp*(M-s)*log(ld(s+1))+0.5*s*(2*M-s)*log(kp);
 end
[q1,hq]=min(MDL);
qq=hq-1;   %信号个数
z=P(:,k);
Vn= z(:,qq+1:M);%噪声子空间 ；

MUSIC算法

部分代码（完整代码见文末）：

X=A*ss+w;   %接受信号模型
XZ=X';
Rx=X(:,1)*XZ(1,:); 
Rx0=zeros(size(Rx,1),size(Rx,2));
for m=1:kp
    Rx=X(:,m)*XZ(m,:);    %求协方差矩阵
    Rx0=Rx0+Rx;
end
Rx1=Rx0/kp;  %求协方差矩阵除以快拍数的均值
y=size(Rx1,1);%协方差矩阵的行数
P=eye(y);   %产生y阶单位阵
%计算A中上半对角线零的个数，确定是否跳出循环
for x=1:1:100     %100为循环足够大的次数
q=0;
    for n=1:1:y-1      %按照论文方法进行循环
        for m=n+1:1:y  
            if round(Rx1(n,m))==0%判断A（n,m）是否为零
                q=q+1;
            end
Y = eye(y);        %生成一个y阶对角矩阵，为U（p，k）做准备
B=Rx1([n m],[n m]);   %取主子阵
b1=B(1,2);           
b2=B(1,1);
b3=B(2,2);
t1=sqrt((real(b1)^2)+(imag(b1)^2));    %分子
t2=real(b1)-imag(b1)*1i;           %分母
T=[1,0;0,t2/t1];               %所求二阶实矩阵
phi=asin(imag(b1)/t2);          %求角phi，未用到
tan=(2*t1/(b2-b3));
thet=atan(tan);            %求角2*theta
theta=thet/2;
Y(n,n)=cos(theta);
Y(n,m)=-sin(theta);
Y(m,n)=sin(theta)*(t2/t1);
Y(m,m)=cos(theta)*(t2/t1); %Y即U（p,k）
dds=Y'*Rx1;
dds1=dds*Y;
Rx1=Y'*Rx1*Y;     %An
P=P*Y;
       end
    end
 if (2*q-y*y+y)==0
     break;       %通过q值来判断是否跳出循环
 end
end
D1=real(diag(Rx1));
[D2,k]=sort(D1,'descend'); %对特征值进行排序并返回结果

结果展示：
在这里插入图片描述

基于正交偶极子的MUSIC算法

正交偶极子模型

与四元数模型区别，这里采用长矢量模型。假设信号源为无限远处单一频率TEM电磁波，其波达方向为 -r，如下图所示。
在这里插入图片描述

承载信号的复基带信号为s(t)，载波频率为,其空间到达角为(θ,φ),且假设该信号为完全极化电磁波，极化信号的两个正交方向的信号幅度比和相位差为(γ,η)。单位矢量(θ,φ,r)构成空间球极坐标系，则该TEM信号可以完全描述为：
空间信号入射正交偶极子示意图
式中：Eφ和Eθ分别为φ和θ方向的极化分量；r为空间接收信号的传播方向的矢量，代表了信号源；k为传播矢量。
空间球极坐标系与空间直角坐标系单位向量之间的转化从关系如下:
在这里插入图片描述
令：u、v1和v2分别为r、θ和φ在空间直角坐标系中的坐标矢量。根据不同坐标参照点下不同坐标之间的转化关系可以得出电场向量在空间笛卡尔坐标中的坐标向量为：

因为磁场与电场正交，因此可以仅用电场表示全部信息，该论文的研究只基于接收电磁波的电场分量。正交偶极子仅可以接收x和y方向的信息，因此最终接收信号坐标矢量矩阵为：
在这里插入图片描述

正交偶极子的阵列接受模型

正交偶极子阵列由相互垂直放置的M个正交偶极子对组成，阵元的位置结构如下图所示，两个相互垂直的正交偶极子分别沿轴和轴方向放置，可以分别接收来自方向和方向的电场分量(可以接收电磁波，但是因为只需要电场信息，因此当作只接收了电场分量，以下的正交偶极子阵列均是如此)。各个阵元构成位置平均分布的圆形阵列，正交偶极子阵元中心距离圆阵的圆心距离为r。
在这里插入图片描述
该阵列的阵列流型矢量为：

正交偶极子的极化导向矢量为：

包含极化信息的表达式为：

空间极化导向矢量为：

该式表示克罗内特积。
假设有M个不相关远场单一频率的完全极化信号入射到该正交偶极子阵列上，接收信号的模型为
在这里插入图片描述
其中A为正交偶极子阵列的极化导向矩阵；S(t)为信号矢量；n(t)为代表噪声的复数矩阵，假设噪声为高斯白噪声。
阵列接收信号的协方差矩阵为

部分代码（完整代码见文末）：

delay=ones(M,Num);%定义延迟矩阵
As=ones(M,Num);%定义阵列流型矢量矩阵
for m=1:M      %阵列流型矢量矩阵赋值
    for n=1:Num  
        delay(m,n)=(cos(rfw(n)-(m-1)*pi/4)*r*cos(rfy(n)))/c;%列向量延迟赋值
        As(m,n)=exp((-2*pi*f*1i).*delay(m,n));%阵列流型矢量赋值
    end      
end
A=ones(2*M,Num);%定义空间极化导向矢量
for n=1:Num   %空间极化导向矢量赋值
wfw=rfw(n);
wfy=rfy(n);
V=[-sin(wfw) cos(wfw)*cos(wfy);cos(wfw) cos(wfy)*sin(wfw)];%导向矢量
wfd=rfd(n);
wxw=rxw(n);
E=[cos(wfd);sin(wfd)*exp(1i*wxw)];%定义极化矢量
Ap=V*E;%极化导向矢量
A(:,n)=kron(As(:,n),Ap);%柯洛内特积
end
AH=zeros(2*M,1);
for m=1:2*M
for n=1:Num
AH(m,1)=AH(m,1)+A(m,n);
end
end

基于正交偶极子的MUSIC算法

已知矩阵为厄尔米特矩阵，对其进行特征值分解，求取其特征向量。取出厄尔米特矩阵的N个特征值，用这些特征值所匹配的特征向量，构造噪声等价矩阵。此处的特征值为较小的N个特征值。
空间谱估计函数为
在这里插入图片描述
令

因为

且

则可得

定义

根据信号子空间的原理，我们可以得到，信号子空间与噪声子空间之间的关系是正交的，因此代表了信号的空间位置信息和电磁波极化信息与噪声子空间之间的关系也是相互正交的,即可得
在这里插入图片描述
当γ∈[0,pi/2]，ap满秩。为了满足上式成立，又由于G不是满秩的，所以可得

因此，信号源的空间位置参数公式为

极化信息可以通过将接收信号的空间位置参数带入阵列接收信号的协方差矩阵进行极化信息的二维谱峰搜索得到。
结果展示：
在这里插入图片描述
极化信息搜索：

模值约束法求极化信息

此外，极化信息可以通过模值约束方法得到。通过正交关系去估计极化信息可以等价为求解一个优化问题。其约束条件可以表示为：
在这里插入图片描述
代价函数为：

要使上式对E求导数使结果恒等于0，通过计算得：

显然，E是G相对应于特征值u的特征向量。由于：

因此，为了让目的函数的值最小，就相当于使特征值u达到最小的预期的E是G的最小的特征值所匹配的特征向量。即：
在这里插入图片描述
按照极化信息的表达式，我们可以根据接下来的等式得到该接收信号的极化信息：

其中，

同时，接收信号的的每一对空间位置参数和极化参数都是相互照应的，不需要进行额外匹配。
部分代码（完整代码见文末）：

fd_j=1:num_i;
xw_j=1:num_i;
fd_jr=1:num_i;
xw_jr=1:num_i;
delay_j=zeros(M,1);%定义极化信息估计时的延迟
As_j=zeros(M,1);%定义极化信息估计时的导向矢量
for k=1:num_i
for m=1:M    %导向矢量赋值
    delay_j(m,1)= (cos(rfw_j(k)-(m-1)*pi/4)*r*cos(rfy_j(k)))/c ;  %定义列向量延迟
    As_j(m,1)=exp((-2*pi*f*1i).*delay_j(m,1));     %定义A列向量
end
wfw_j=rfw_j(k);
wfy_j=rfy_j(k);
V_j=[-sin(wfw_j) cos(wfw_j)*cos(wfy_j);cos(wfw_j) cos(wfy_j)*sin(wfw_j)]; 
A_j=kron(As_j(:,1),V_j);%空间极化导向矢量
G_j=A_j'*(Vn*Vn')*A_j;%空间导向矢量与噪声子空间的乘
y_j=size(G_j,1);%协方差矩阵的行数
P_j=eye(y_j);   %产生y阶单位阵
%计算A中上半对角线零的个数，确定是否跳出循环
for x=1:1:1000     %100为循环足够大的次数
q=0;
    for n=1:1:y_j-1      %按照论文方法进行循环
        for m=n+1:1:y_j  
            if round(G_j(n,m))==0;%判断A（n,m）是否为零
                q=q+1;
            end
Y = eye(y_j);        %生成一个y阶对角矩阵，为U（p，k）做准备
B_j=G_j([n m],[n m]);   %取主子阵
b1=B_j(1,2);           
b2=B_j(1,1);
b3=B_j(2,2);
t1=sqrt((real(b1)^2)+(imag(b1)^2));    %分子
t2=real(b1)-imag(b1)*1i;           %分母
phi=asin(imag(b1)/t2);          %求角phi，未用到
tan=(2*t1/(b2-b3));
thet=atan(tan);            %求角2*theta
theta=thet/2;
Y(n,n)=cos(theta);
Y(n,m)=-sin(theta);
Y(m,n)=sin(theta)*(t2/t1);
Y(m,m)=cos(theta)*(t2/t1); %Y即U（p,k）
dds=Y'*G_j;
dds1=dds*Y;
G_j=Y'*G_j*Y;     %An
P_j=P_j*Y;
       end
    end
 if (2*q-y_j*y_j+y_j)==0
     break;       %通过q值来判断是否跳出循环
 end
end
T_j=real(diag(G_j));
[T_sort_j,ss_j]=sort(T_j,'descend'); %对特征值进行排序并返回结果
E_j=P_j(:,ss_j(size(ss_j,1)));%
fd_jr(k)=atan(abs(E_j(2,1)/E_j(1,1)));%计算极化幅度信息
fd_j(k)=fd_jr(k)/radian;
xw_jr(k)=angle(E_j(2,1)/E_j(1,1));%计算极化相位信息
xw_j(k)=xw_jr(k)/radian;
end

结果展示：
在这里插入图片描述

基于正交偶极子的四元数MUSIC算法

四元数的阵列接受模型

假设空间有一个无限远处的单一频率TEM信号入射到阵列上，承载信号的复基带信号为s(t)，载波频率为f0,其空间到达角为(θ,φ ),且假设该信号为完全极化波，极化信号的两个正交方向的信号幅度比和相位差为(γ,η)，对于如图2.2摆放的电磁偶极子阵列，其两个相互垂直的阵元分量的极化导向矢量可以表示为
在这里插入图片描述
两个相互垂直的阵元分量接收信号的两个量可以表示为

为了保持两路极化信号的分量之间稳定的正交关系(因为是两个正交方向的信号幅度比和相位差稳定的信号)，用一个四元数来表示两个相互垂直的正交偶极子阵元的接收信号，如下：
在这里插入图片描述
式中，P为极化导向矢量，式中的数全是以四元数形式存在。通过以上推导，运用以四元数形式存在的接收信号模型可以得出正交偶极子均匀圆阵的接收信号矩阵如下：

式中，as为正交偶极子均匀圆阵的空间导向矢量，φ表示接收信号入射至两个相邻的正交偶极子阵元的到达时间，也称作空间相移因子；xt表示正交偶极子均匀放置的圆阵的接收信号的四元数矢量矩阵。当空间有Ｋ个信号源且仅考虑高斯白噪声的情况下：
在这里插入图片描述
式中，dk为以四元数形式存在，由极化导向矢量和空间导向矢量组成，被称为极化域-空域联合导向矢量；D为极化域-空域联合导向矩阵；st为入射信号复基带矩阵；et为以四元数形式存在的高斯白噪声分量，e1t为方向上的复数形式的噪声矢量，e2t为方向上的复数形式的噪声矢量。
雷达阵列的入射电磁波的相关矩阵可表示为
在这里插入图片描述
式中，三角符号的形式代表了共轭转置矩阵；Re表示雷达阵列接收的总噪声的相关系数矩阵；Rs表示输入信号的自相关系数矩阵。可见自相关系数矩阵R的秩与Rs的秩有很大的关系，被Rs的秩所约束。

部分代码（完整代码见文末）：

delay=ones(M,Num);%定义延迟矩阵
As=ones(M,Num);%定义阵列流型矢量矩阵
for m=1:M      %阵列流型矢量矩阵赋值
    for n=1:Num  
        delay(m,n)=(cos(rfw(n)-(m-1)*pi/4)*r*cos(rfy(n)))/c;%列向量延迟赋值
        As(m,n)=exp((-2*pi*f*1i).*delay(m,n));%阵列流型矢量赋值
    end      
end
A=ones(2*M,Num);%定义空间极化导向矢量
for n=1:Num   %空间极化导向矢量赋值
wfw=rfw(n);
wfy=rfy(n);
V=[-sin(wfw) cos(wfw)*cos(wfy);cos(wfw) cos(wfy)*sin(wfw)];%导向矢量
wfd=rfd(n);
wxw=rxw(n);
E=[cos(wfd);sin(wfd)*exp(1i*wxw)];%定义极化矢量
Ap=V*E;%极化导向矢量
A(:,n)=kron(As(:,n),Ap);%柯洛内特积
end
AH=zeros(2*M,1);
for m=1:2*M
for n=1:Num
AH(m,1)=AH(m,1)+A(m,n);
end
end

四元数MUSIC算法

对四元数模型下阵列接收信号的协方差矩阵进行四元数特征值分解：
在这里插入图片描述
由空间谱估计理论可知，噪声的全部信息都保存在特征向量UN之中，同时不含有入射信号的信息，接收信号的全部信息包含在特征向量US之中。因此，可以通过UN估计出入射信号的全部信息。可得，空间谱估计函数为
在这里插入图片描述
可以通过四维谱峰搜索得到接收信号的空间位置信息，但是四维谱峰搜索算法所需运算量过大，工程可实现性有限。因此，需要对其进行降维。
定义

设

,可以明显地发现T中包含入射电磁波的空间位置和极化信息，而C中只有入射电磁波的空间位置信息。此外，γ∈[0,pi/2]时，P满秩，则此时C=0。
因此，信号的空间位置信息的空间谱估计函数为
在这里插入图片描述
此时，通过二维谱峰搜索即可求得接收信号的空间位置信息。
但是，在四元数模型下极化导向矢量由复数域内的复数矢量变为四元数域的四元数矢量，极化导向矢量被信号矢量吸收，之后才进行了四元数矩阵分解，因此不能通过此方法求得极化信息。但是可以通过构建长矢量MUSIC算法的空间谱估计函数来求得极化信息，将接收信号的空间位置参数带入空间谱估计函数，进行二维谱峰搜索便可以得到极化信息。该算法利用了不同电场方向的正交特性，相比于长矢量MUSIC算法提升了估计性能。

部分代码（完整代码见文末）：

Xa=A*s+w;   %接受信号模型
%XX=sr+wr
Xx=ones(M,kp);
Xy=ones(M,kp);
for n=1:M
Xx(n,:)=Xa(2*n-1,:);
Xy(n,:)=Xa(2*n,:);
end
X=zeros(2*M,2);
Rx=X*X';%协方差矩阵
Rx_p=zeros(size(Rx,1),size(Rx,2));%定义协方差矩阵的和
for n=1:kp
    X=[Xx(:,n),Xy(:,n);-conj(Xy(:,n)),conj(Xx(:,n))];
    Rx=X*X';  %求一次的协方差矩阵  
    Rx_p=Rx_p+Rx;      %协方差矩阵求和
end
Rx_j=Rx_p/kp;  %求协方差矩阵除以快拍数的均值
y=size(Rx_j,1);%协方差矩阵的行数
P=eye(y);   %产生y阶单位阵
%计算A中上半对角线零的个数，确定是否跳出循环
for x=1:1:1000     %1000为循环足够大的次数
q=0;
    for n=1:1:y-1      %按照论文方法进行循环
        for m=n+1:1:y  
            if round(Rx_j(n,m))==0;%判断A（n,m）是否为零
                q=q+1;
            end
Y = eye(y);        %生成一个y阶对角矩阵，为U（p，k）做准备
B=Rx_j([n m],[n m]);   %取主子阵
b1=B(1,2);           
b2=B(1,1);
b3=B(2,2);
t1=sqrt((real(b1)^2)+(imag(b1)^2));    %分子
t2=real(b1)-imag(b1)*1i;           %分母
phi=asin(imag(b1)/t2);          %求角phi，未用到
tan=(2*t1/(b2-b3));
thet=atan(tan);            %求角2*theta
theta=thet/2;
Y(n,n)=cos(theta);
Y(n,m)=-sin(theta);
Y(m,n)=sin(theta)*(t2/t1);
Y(m,m)=cos(theta)*(t2/t1); %Y即U（p,k）
dds=Y'*Rx_j;
dds1=dds*Y;
Rx_j=Y'*Rx_j*Y;     %An
P=P*Y;          %特征向量
       end
    end
 if (2*q-y*y+y)==0
     break;       %通过q值来判断是否跳出循环
 end
end
T=real(diag(Rx_j));%取协方差矩阵的特征值
[T_sort,ss_T]=sort(T,'descend'); %对特征值进行排序并返回结果

结果展示：
在这里插入图片描述