银川网站建设那家好,自建小程序商城,南阳专业做网站公司哪家好,h5页面制作工具是什么MUSIC算法原理及仿真 DOA波达方向估计MUSIC算法概述MUSIC算法原理MUSIC算法MATLB仿真 DOA波达方向估计
DOA#xff08;Direction Of Arrival#xff09;波达方向是指通过阵列信号处理来估计来波的方向#xff0c;这里的信源可能是多个#xff0c;角度也有多个。DOA技术主要… MUSIC算法原理及仿真 DOA波达方向估计MUSIC算法概述MUSIC算法原理MUSIC算法MATLB仿真 DOA波达方向估计
DOADirection Of Arrival波达方向是指通过阵列信号处理来估计来波的方向这里的信源可能是多个角度也有多个。DOA技术主要有ARMA谱分析、最大似然法、熵谱分析法和特征分解法特征分解法主要有MUSIC算法、ESPRIT算法WSF算法等。其中MUSIC是一种经典的DOA估计算法。
MUSIC算法概述
MUSIC(multiple signal classification algorithm)算法是一种基于矩阵特征空间分解的方法。从几何角度讲信号处理的观测空间可以分解为信号子空间和噪声子空间显然这两个空间是正交的。信号子空间由阵列接收到的数据协方差矩阵中与信号对应的特征向量组成噪声子空间则由协方差矩阵中所有最小特征值噪声方差对应的特征向量组成。
MUSIC算法原理
MUSIC算法是空间谱估计测向理论的重要基石。算法原理 如下 1 不管测向天线阵列形状如何也不管入射来波入射角的维数如何如图1所示假定阵列由 N N N个阵元组成并且符合窄带、远场模型。 图1. 均匀线性阵模型 则阵列输出模型的矩阵形式都可以表示为 Y ( t ) A S N ( t ) Y(t)ASN(t) Y(t)ASN(t) 其中 Y Y Y是观测到的阵列输出数据复向量 X X X是未知的空间信号复向量 N ( t ) N(t) N(t)是阵列输出向量中的加性噪声 A A A是阵列的方向矩阵又称之为阵列流形矩阵当我们以最左侧的阵元为参考位此处 A A A矩阵表达式表示为 A [ α ( θ 1 ) , α ( θ 2 ) , ⋯ , α ( θ M ) ] [ 1 1 ⋯ 1 e − j ϕ 1 e − j ϕ 2 ⋯ e − j ϕ D e j − 2 ϕ 1 e − j 2 ϕ 2 ⋯ e − j 2 ϕ D ⋮ ⋮ ⋱ ⋮ e − j ( N − 1 ) ϕ 1 e − j ( N − 1 ) ϕ 2 ⋯ e − j ( N − 1 ) ϕ D ] A[\alpha(\theta_1),\alpha(\theta_2),\cdots,\alpha(\theta_M)]\begin{equation} \begin{bmatrix} 1 1 \cdots 1 \\ e^{-j\phi_1} e^{-j\phi_2} \cdots e^{-j\phi_D} \\ e^{j-2\phi_1} e^{-j2\phi_2} \cdots e^{-j2\phi_D} \\ \vdots \vdots \ddots \vdots \\ e^{-j(N-1)\phi_1} e^{-j(N-1)\phi_2} \cdots e^{-j(N-1)\phi_D} \end{bmatrix} \end{equation} A[α(θ1),α(θ2),⋯,α(θM)] 1e−jϕ1ej−2ϕ1⋮e−j(N−1)ϕ11e−jϕ2e−j2ϕ2⋮e−j(N−1)ϕ2⋯⋯⋯⋱⋯1e−jϕDe−j2ϕD⋮e−j(N−1)ϕD 可以看到 A A A是一个范德蒙矩阵其中 ϕ k 2 π s i n θ k / λ \phi_k2{\pi}sin\theta_k/\lambda ϕk2πsinθk/λ且 θ k , 1 ≤ k ≤ D \theta_k,1\le k \le D θk,1≤k≤D为第 k k k个信源的角度共有 D D D个角度。 ϕ k \phi_k ϕk在这里是指波程差引起的相移如下图所示相邻阵元之间的波程差为 d s i n θ dsin\theta dsinθ将此转换为相位单位为弧度则为 2 π s i n θ k / λ 2{\pi}sin\theta_k/\lambda 2πsinθk/λ λ \lambda λ为信号波长。由下图的等相位面可知信号到达最左侧的阵元需要走更远的距离因此以最左侧的阵元为参考即左侧阵元接收到的信号为 x 1 ( t ) x_1(t) x1(t)则其他阵元接收到的信号都滞后于它此时 x 1 ( t ) s ( t ) , x 2 ( t ) s ( t ) e − j ϕ x_1(t)s(t),x_2(t)s(t)e^{-j\phi } x1(t)s(t),x2(t)s(t)e−jϕ。如果我们以最右侧信号为参考即最右侧阵元接收到的信号为 x 1 ( t ) x_1(t) x1(t)则 x 1 ( t ) s ( t ) , x 2 ( t ) s ( t ) e j ϕ x_1(t)s(t),x_2(t)s(t)e^{j\phi } x1(t)s(t),x2(t)s(t)ejϕ。注意角度的正负号。 图2. 波程差示意图 接收信号写为向量形式 Y ( t ) [ y 1 ( t ) y 2 ( t ) y 3 ( t ) ⋮ y N ( t ) ] Y(t)\begin{equation} \begin{bmatrix} y_1(t) \\ y_2(t) \\ y_3(t) \\ \vdots \\ y_N(t) \end{bmatrix} \end{equation} Y(t) y1(t)y2(t)y3(t)⋮yN(t) y i ( t ) y_i(t) yi(t)为第 i i i个阵元接收到的信号共有 N N N个阵元。 S [ s 1 ( t ) , s 2 ( t ) , ⋯ , s D ( t ) ] S[s_1(t),s_2(t),\cdots,s_D(t)] S[s1(t),s2(t),⋯,sD(t)] s i ( t ) , 1 ≤ i ≤ D s_i(t),1\le i \le D si(t),1≤i≤D表示第 i i i个信源发出的信号。 N ( t ) [ n 1 ( t ) , n 2 ( t ) , ⋯ , n N ( t ) ] T N(t)[n_1(t),n_2(t),\cdots,n_N(t)]^T N(t)[n1(t),n2(t),⋯,nN(t)]T表示噪声符合独立同分布的假设。 MUSIC算法的处理任务就是设法估计出入射到阵列的空间信号的个数 D D D以及空间信号源的强度及其来波方向 θ k \theta_k θk。 2 在实际处理中 Y Y Y得到的数据是有限时间段内的有限次数的样本也称快拍或快摄样本的长度取决于采样率和采样时间在这段时间内假定来波方向不发生变化且噪声为与信号不相关的白噪声则定义阵列输出信号的二阶矩 R y R_y Ry。 R y E [ Y Y H ] R_yE[YY^H] RyE[YYH] 其中 [ ∙ ] H [\bullet]^H [∙]H表示共轭转置值得注意的是在MATLAB中我们常用的转置符号‘其实是共轭转置而普通转置为(.)当然他们对于实数矩阵而言没有任何区别。 E [ ∙ ] E[\bullet] E[∙]表示统计学上的期望。 由于噪声信号符合独立同分布的假设信号与噪声互不相关噪声为零均值白噪声则上式又可以写为 R y E [ Y Y H ] E [ ( A S N ) ( A S N ) H ] A E [ S S H ] A H E [ N N H ] A R S A H R N R_yE[YY^H]E[(ASN)(ASN)^H]AE[SS^H]A^HE[NN^H]AR_SA^HR_N RyE[YYH]E[(ASN)(ASN)H]AE[SSH]AHE[NNH]ARSAHRN 其中 R S E [ S S H ] R_SE[SS^H] RSE[SSH]为信号相关矩阵 R N E [ N N H ] σ 2 I R_NE[NN^H]\sigma^2I RNE[NNH]σ2I为噪声相关矩阵 I I I为 N ∗ N N^*N N∗N的单位矩阵。 3 MUSIC算法的核心就是对 R y R_y Ry进行特征值分解利用特征向量构建两个正交的子空间即信号子空间和噪声子空间。对 R y R_y Ry进行特征分解。 4 U U U是非负定的厄米特矩阵所以特征分解得到的特征值均为非负实数有 D D D个大的特征值和 M − D M-D M−D个小的特征值大特征值对应的特征向量组成的空间Us为信号子空间小特征值对应的特征向量组成的空间 U n U_n Un为噪声子空间。 5 将噪声特征向量作为列向量组成噪声特征矩阵 并张成 M − D M-D M−D维的噪声子空间 U n U_n Un噪声子空间与信号子空间正交。而 U s U_s Us的列空间向量恰与信号子空间重合所以 U s U_s Us的列向量与噪声子空间也是正交的由此可以构造空间谱函数。 6 在空间谱域求取谱函数最大值其谱峰对应的角度即是来波方向角的估计值。 P m u s i c 1 α H ( θ ) U n U n H α ( θ ) P_{music}\frac{1}{\alpha^H(\theta)U_nU^H_n\alpha(\theta)} PmusicαH(θ)UnUnHα(θ)1
MUSIC算法MATLB仿真
% MUSIC空间谱估计算法仿真
close all; clear all; clc;
Jsqrt(-1);
source_number4;
source_doa[30 40 50 70];
sensor_number7;
snapshot_number2000;
snr10; Aexp(-J*(0:sensor_number-1)*pi*sin(source_doa*pi/180));
s(randn(source_number,snapshot_number)J*randn(source_number,snapshot_number))/sqrt(2);
xA*s;
yawgn(x,snr,measured);
Ry*y/snapshot_number;
[V,D]eig(R);
UnV(:,1:sensor_number-source_number);
GnUn*Un; searching_doa0:0.1:90;
for i1:length(searching_doa) a_thetaexp(-J*(0:sensor_number-1)*pi*sin(pi*searching_doa(i)/180))P_con(i)abs(a_theta*R*a_theta);P_BF(i)abs((a_theta*R*a_theta)./(a_theta*a_theta)); P_capon(i)1./abs((a_theta*inv(R)*a_theta)); P_music(i)1./abs((a_theta*Gn*a_theta));
end
plot(searching_doa,P_con/max(P_con),k);hold on;
plot(searching_doa,P_BF/max(P_BF),r); hold on;
plot(searching_doa,P_capon/max(P_capon),g); hold on;
plot(searching_doa,P_music/max(P_music),b); hold off;grid on;
xlabel(ang);
ylabel(功率谱估计);
legend(conditional spectrum,Bartlett spectrum,Capon spectrum,Music spectrum);仿真结果如下 图3. 仿真结果 还有一种参考程序
clear; close all;
%%%%%%%% MUSIC for Uniform Linear Array%%%%%%%%
derad pi/180; %角度-弧度
N 8; % 阵元个数
M 3; % 信源数目
theta [-30 0 60]; % 待估计角度
snr 10; % 信噪比
K 512; % 快拍数dd 0.5; % 阵元间距
d0:dd:(N-1)*dd;
Aexp(-1i*2*pi*d.*sin(theta*derad)); %方向矢量%%%%构建信号模型%%%%%
Srandn(M,K); %信源信号入射信号
XA*S; %构造接收信号
X1awgn(X,snr,measured); %将白色高斯噪声添加到信号中
% 计算协方差矩阵
RxxX1*X1/K;
% 特征值分解
[EV,D]eig(Rxx); %特征值分解
EVAdiag(D); %将特征值矩阵对角线提取并转为一行
[EVA,I]sort(EVA); %将特征值排序 从小到大
EVfliplr(EV(:,I)); % 对应特征矢量排序% 遍历每个角度计算空间谱
for iang 1:361angle(iang)(iang-181)/2;phimderad*angle(iang);aexp(-1i*2*pi*d*sin(phim)).; EnEV(:,M1:N); % 取矩阵的第M1到N列组成噪声子空间Pmusic(iang)1/(a*En*En*a);
end
Pmusicabs(Pmusic);
Pmmaxmax(Pmusic)
Pmusic10*log10(Pmusic/Pmmax); % 归一化处理
hplot(angle,Pmusic);
set(h,Linewidth,2);
xlabel(入射角/(degree));
ylabel(空间谱/(dB));
set(gca, XTick,[-90:30:90]);
grid on;仿真结果如下 图4. 仿真结果
