中期报告
2020-03-02 12:07:13 来源：范文大全收藏下载本文
基于music算法的DOA仿真及探究
姓名：陈望桥 廖正松
院系：工学院
学号：1200011011
指导老师：黄迅
摘要：基于传统music算法进行DOA估计与仿真，并根据music算法提出二维波达方向角跟踪的初步模型；在最后对数字麦克风（sigma-delta）的在DOA中的使用进行模拟与初步的仿真探索，并提出假设。
关键词：DOA跟踪，music算法，sigma-delta工作原理，数字麦克风；
引言：波达方向估计，又称为谱估计、波达角估计。一个信源有很多可能的传播路径和到达角。如果几个发射机同时工作，每个信源在接收机处形成潜在的多径分量。因此，接收天线能估计出这些到达角就显得很重要，目的是破译出哪个发射机在工作以及发射机所处的可能位置。
通过测量辐射信号的波达方向来估测辐射源的位置。理论上这种估计只需要两个接收阵元就可以确定辐射源的位置，但在实际中，由于受到角度分辨率、多径和噪声限制，所需阵元通常要多于两个。
Music算法是20世纪中后期崛起的一种利用了信号子空间和噪声子空间的正交性，构造空间谱函数，通过谱峰搜索，检测信号的DOA算法．它具有此前传统列阵算法所不具备的优点，利用空间正交性可以大量减少运算并提高运算速度。并且也摆脱了传统阵列现行排布的局限。最具有代表性的就是基于music算法的圆形阵列DOA估计。在music算法问世后涌现了许多根据空间正交性进行DOA运算的算法。所以Music算法具有很强的代表性。
一 Music算法概述 假设在空间中某点存在一个声源，该声源向各个方向发出声音信号，并在空间中另一位置存在以一组以已知顺序排列的接受器，接受声源发出的声音信号。在声源据接收器组足够远时可以近似认为声源发出的声波在接收器位置为平面波。波达方向角是空间信号的到达方向（各个信号到达阵列参考阵元的方向角，简称波达方向）。估计波达方向角的算法称为DOA算法。
MUSIC算法是一种基于矩阵特征空间分解的方法。从几何角度讲，信号处理的观测空间可以分解为信号子空间和噪声子空间，显然这两个空间是正交的。信号子空间由阵列接收到的数据协方差矩阵中与信号对应的特征向量组成，噪声子空间则由协方差矩阵中所有最小特征值（噪声方差）对应的特征向量组成。MUSIC算法就是利用这两个互补空间之间的正交特性来估计空间信号的方位。噪声子空间的所有向量被用来构造谱，所有空间方位谱中的峰值位置对应信号的来波方位。MUSIC算法大大提高了测向分辨率，同时适应于任意形状的天线阵列，但是原型MUSIC算法要求来波信号是不相干的。
假设，空间中存在的信源点在O点（阵列参考点）产生的信号为s，第i个信源点产生信号为si，由于接收器阵列几何分布与波达方向角共同影响，导致各个接收器接收信号与参考信号之间具有相位差，这一相位差与波达方向角有关。设X为接受信号向量；A为阵列方向矩阵；S为信源向量；[1] 其中 A=[ a1,a2,a3,….aN]；ai阵元对第i个信号的响应；N为信源数； 有：X=A*S+N； 对于圆形阵列有aij=exp(i*w*r*cos(∅i)*cos(−θi)/c); 其中 θj=2∗n∗j θi 为经度角∅为纬度角
由于各阵元的噪声互不相关，且也与信号不相关，因此接收数据X的协方差矩阵为 R=E（x*x’）；因为R=A*P*A’+σ2*I；P为信号的协方差阵；P=E(s*s’)；σ2＞0，而APA’的特征值为正，R为满秩阵，因此R有n个正特征值u1,u2,….un，它们所对应的特征向量为l1,l2,l3,….ln，且各特征向量是相互正交的，这些特征向量构成n*n维空间的一组正交基。与信号有关的特征值有N个，且u1…uN为P的特征值与σ2之和，而矩阵的其余En=（n-N）个特征值为σ2，也就是说σ2R的最小特征值，它是（n-N）重的。因此只要将天线各阵元输出数据的协方差矩阵进行特征值分解，找出最小特征值的个数En，据此就可以求出信号源的个数N，即有
由于R的最小特征向量仅与噪声有关，因此由这En个特征向量所张成的子空间称之为噪声子空间，而与它正交的子空间，即由信号的方向向量张成的子空间则是信号子空间。将矩阵R所在的n*n维空间分解成两个完备的正交子空间，信号子空间和噪声子空间。且易知两空间正交。利用正交性，可以用En个小特征值对应的En个噪音特征向量与信源向量正交，进而得到：U=[lN+1,….ln]’；U*ai=0；[1] 由于协方差矩阵R是根据有限次观测数据估计得到的，对其进行特征分解时，最小特征值(噪声方差)和重数En的确定以及最小特征向量的估计都是有误差的，当U存在偏差时，式U*ai右边不是零向量。这时，可取使得U*ai(θ，∅)的2-范数为最小值的kˆ(θ，∅)作第k个信号源方向的估值。连续改变(θ，∅)值，进行谱峰搜索，由此得到N个最小值所对应的
π(θ，∅)就是N个信号源的位置角度。通常做法是利用噪声子空间与信号子空间的正交性，构造如下空间谱函数
=
1 ，∅ ∗∗′∗ ，∅ ′
[1] 对P进行峰值搜索即可得到N个信号的波达方向角。
二 music算法模拟
利用MATLAB软件，生成信号并进行模拟实验。 模拟生成信号：
s0=exp(j*2*pi*(f0*t+0.5*u0*t.^2)); s1=exp(j*2*pi*(f0*t+0.5*u1*t.^2)); s2=exp(j*2*pi*(f0*t+0.5*u2*t.^2)); s3=exp(j*2*pi*(10*t+0.5*u3*t.^2))； S=[s0,s1,s2,s3]; 对于圆形阵列，阵列方向矩阵A：
a_theta0=exp(j*2*pi*d*cos(thetaa(1)/180*pi)*cos((thetaa(4)-a*e)/180*pi)/lamda0); a_theta1=exp(j*2*pi*d*cos(thetaa(2)/180*pi)*cos((thetaa(5)-a*e)/180*pi)/lamda0); a_theta2=exp(j*2*pi*d*cos(thetaa(3)/180*pi)*cos((thetaa(6)-a*e)/180*pi)/lamda0); a_theta3=exp(j*2*pi*d*cos(thetaa(7)/180*pi)*cos((thetaa(8)-a*e)/180*pi)/lamda0); A=[a_theta0 a_theta1 a_theta2 a_theta3]; X=A*S+N; 根据music算法，得到普空间函数的图像为：
可清楚看到有四个峰值。并经过检验得这四个峰值对应的经度纬度角与模拟信号给出的经度纬度角相同。 对此我们针对music算法可以得到如下结论： 1 计算量大，对于全平面进行遍历，缺少表达式。 2不利于进行实时反馈不利于进行实际操作。 3 对于相干源缺少解决办法。 4 对于非相干源具有超分辨能力。 5 对于阵列排列要求很低
6具有比较良好的低信噪比下的分辨能力。
三 DOA跟踪算法（基于协方差矩阵跟踪算法）
3.1 选择协方差矩阵跟踪算法原因
由于实际中我们无法预知信源发声频率，进而难以确定均匀线阵的阵元间距。从而不防使用Music算法，采用均匀圆阵，采用协方差跟踪算法可以很好的根据协方差矩阵估计出信源运动方向，从而减少了大量的计算。
3.2算法原理
在每一个周期内拥有M个采样，每个周期时长为T，假设在每个周期内信源位置不变，则有对于KT与（K+1)T时刻信源位移为小量，可根据协方差差值计算得到该位移。
3．3二维协方差矩阵跟踪算法
算法原理（基于均匀圆阵）
利用每一个采样周期得到每一个周期内的协方差Rk，通过Rk差值的泰勒展开保留的一次项得到所需的线性方程组。进而得到估计的位置。 通过泰勒展开及矩阵运算有：
H=[H1,H2;H3,H4] H1(I,j)=γ1， +1
H2(I,j)=μH1(I,j)=γ1，+1 +1 ，1H2(I,j)=μμ，+1，1
=∗∗cos ∅ ∗ cos −∗ −cos −∗ ∗sin ∅ ∗(cos −∗
+cos −∗ )
γ=∗∗cos ∅ ∗ cos −∗ −cos −∗ ∗cos ∅ ∗(sin −∗
+sin −∗ )
α∗=∗/ ∆R=[R i+1 −R i ]
∆x= ∆∅1,∆∅2….∆∅N,∆θ1,∆θ2……∆θN ′ ,∆y=[∆R12,∆R13……∆R1n,∆R21,∆R22……∆R2n]′
∆x=(H′∗H)−1∗′∗∆
H通常为非奇异阵，可得：
∆x=H−1∗∆
由上可知：在得到初始波达角后，可以不再进行矩阵特征分解即可得到波达角变化量。并且波达角与协方差矩阵变化量存在线性关系。
根据上诉步骤我们可以得知，对于普遍的，经纬度变化率不大的情况下这种方法较为实用。对于经纬度变化率的变化较大，即泰勒展开时不能忽略二阶项时，上述过程应该进行更改，并不能写成线性表达。
四 基于sigma-delta ADC的阵列探索
越来越多的应用，例如过程控制、称重等，都需要高分辨率、高集成度和低价格的ADC、新型∑-△转换技术恰好可以满足这些要求。然而，很多设计者对于这种转换技术并不十分了解，因而更愿意选用传统的逐次比较ADC。∑-△转换器中的模拟部分非常简单(类似于一个1bit ADC)，而数字部分要复杂得多，按照功能可划分为数字滤波和抽取单元。由于更接近于一个数字器件，∑-△ADC的制造成本非常低廉。
4.1 sigma-delta ADC运用基础
针对DOA实现，我们通常采用模拟信号接收器，然后通过ADC转化成数字信号进行处理。但是如果我们抛弃模拟信号接收器，直接采用sigma-deltaADC就可以省去模数转换过程进而节省很多成本。
在数字麦克风的使用中，我们对于采集到的高速数字信号的处理的意义。我认为，数字麦克风与传统麦克风的主要区别和优势也在于它可以输出超高速信号的特性，但是数字麦克风输出的信号并没有位数的概念。数字麦克风它只是保留了频域信息，并且可以通过频域信息还原成初始信号的，那么，我们在使用中，没有任何的理由可以说将其认为成n位二进制信息，因为它本事并不是作为二进制信息生成的。
对于这一特点，进行了一些仿真和计算，发现，sigma-delta过程完美的将信噪比降低并将噪音信息推至高频，那么，可以通过它在频域中所保留的特性来进行信息的比对。对比了一下之前的传统Doa算法的过程，在这些算法中，无一例外都是通过采集到的信息的强度来进行分析而得到相位的差值。
但是根据FFT变换中延迟定理，同一信源点发出信息，在忽略幅值变化的时候，它的时延是可以在频域中体现的，我们可以根据频域特性来求出相位的变化，并且不再需要滤波过程和进制转换过程。
通过MAtlAB进行了一些仿真，实现一个模拟圆形阵列数字麦克风的在频域中分析相位差的程序，并且计算结果和设定的相位差基本相同。
4.2 MATLAB模拟sigma-delta 阵列DOA 1 基准初始信号：s(m)=exp(1i*pi*m*2*w/n)；保留实部即为cos函数 2 sigma-delta ADC转换： 转换后信号图像如下（w=10）
3 模拟接收信号，并进行相位差运算
（麦克风数 10 经度角 40 纬度角 20 半径 100 频率 1000） 接收器还原出的相位差：
106.2394
37.9350 114.4446 128.3194
99.0509 106.2394 114.4446 128.3194
99.0509；
实际模拟生成相位差：
105.8438
37.5729 114.1146 127.9471
98.6451 105.8438 114.1146 127.9471
98.6451；
构造范数并输出图像：
37.9350 37.5729
可知存在一个声源点，且经纬度与设立相同。
4.3 利用sigma-deltaADC DOA与传统DOA差别
对于静止声源点，sigma-delta ADC DOA算法中缺失了经度角（-360,-180）即，目前只能完成半个空间的DOA估计。并且对于相位差较小（r较小）的情况误差大于传统算法。但是对于较大的相位差（r较大）精确度相比传统算法大幅增高。并由于sigma-deltaADC自身特性，不再存在迭代更新的过程。对于动声源点，sigma-deltaADC DOA可以快速准确的得到动态经纬度。这一特点是传统DOA算法不具备的。并且sigma-deltaADC DOA 采用FFT算法进行运算，摆脱了大规模矩阵运算的运算量，使得运算时间减少，并且可以根据FFT算法直接得到各个声源点的声源频率。对于多个相干声源的DOA计算仍然缺乏较好的解决办法。
五 预期与假设
由于sigma-delta ADC具有的高速数据传递特性，频域完整的特性以及噪音频率特性我么可以利用频域进行DOA估计。并且这些特性保证了sigma-delta ADC 可以实时传递信息，采样频率极高，保证了动态DOA估计的快速性与准确性。
由于sigma-delta ADC 输出信号是一维数字信号，没有位数与分辨率的概念，导致了使用sigma-delta ADC 如果采用传统算法的数据意义的模糊。所以我认为采用频域分析是一个正确的选择。并且有某种方式使得缺失的经度角被补充回来。并且有某种方式可以避免低相位差时高误差现象的出现。
处于成本以及时间等等方面考虑，数字麦克风或sigma-deltaADC 在DOA中的使用时不可避免的，我们应该专注于sigma-deltaADC 的工作过程以及它的种种特性来设计并完善现有DOA算法。最终达到快速，准确，稳定的目标。 参考文献：
[1]Schmidt R.O.Multiple emitter location and signal parameter estimation [J].IEEE Trans AP,1986
附录一 music算法仿真程序（MATLAB） %阵列接收数据： clear; clc; close all; tic;
M=10;L=1024;N=4; thetaa=[15 20 40 0 25 21 30 46]; snr=20;
f0=10e6;C=3e8;fs=200;Ts=1/fs; lamda0=C/f0;d=5*lamda0;
%%%%%%%%%%%% 产生信号%%%%%%%%%%%%%%%% u0=-15;u1=0;u2=15;u3=30; t=(0:L-1)*Ts;
s0=exp(j*2*pi*(f0*t+0.5*u0*t.^2)); s1=exp(j*2*pi*(f0*t+0.5*u1*t.^2)); s2=exp(j*2*pi*(f0*t+0.5*u2*t.^2)); s3=exp(j*2*pi*(10*t+0.5*u3*t.^2)) S=[s0;s1;s2;s3;];
%%%%%%%%%%% 生成导向矢量矩阵%%%%%%%% a=[0:M-1]\'; e=360/M; a_theta0=exp(j*2*pi*d*cos(thetaa(1)/180*pi)*cos((thetaa(4)-a*e)/180*pi)/lamda0); a_theta1=exp(j*2*pi*d*cos(thetaa(2)/180*pi)*cos((thetaa(5)-a*e)/180*pi)/lamda0); a_theta2=exp(j*2*pi*d*cos(thetaa(3)/180*pi)*cos((thetaa(6)-a*e)/180*pi)/lamda0); a_theta3=exp(j*2*pi*d*cos(thetaa(7)/180*pi)*cos((thetaa(8)-a*e)/180*pi)/lamda0); A=[a_theta0 a_theta1 a_theta2 a_theta3]; X0=A*S;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% real_noise=randn(size(X0)); imag_noise=randn(size(X0));
noise=(real_noise+j*imag_noise)/2^0.5; noise=1^(-snr/20)*noise;
%%%%% 阵元接受数据 %%%%%%%%%%%% X0=X0+noise;
%%%%%%%%无噪声经典MUSIC算%%%%%%%%%% Rx0=X0*X0\'; [V0,D0]=eig(Rx0);%%%%%%%% 求矩阵的特征值%%%%%%% E0=diag(D0);%%%%%%%% 提取D矩阵的对角线构成一个一列矩阵%%%%%%%%
Un0=V0(:,1:M-N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta=-5:0.05:60;%%%%% 从-90--90，每格取0.01 %%%% for ff=1:length(theta)
thetab=-5:0.05:60;
for gg=1:length(thetab) a_theta=exp(j*2*pi*d/lamda0*cos(theta(ff)/180*pi)*cos((thetab(gg)-a*e)/180*pi)); Pmusic0(ff,gg)=100*log10(1/abs((a_theta\'*Un0*Un0\'*a_theta)));
end end %生成三维图； xa=-5:0.05:60; ya=xa; [theta,thetab] = meshgrid(xa,ya); surf(theta,thetab,Pmusic0); colormap(\'jet\'); shading interp; toc;
附录二 sigma-delta ADC DOA仿真程序 function [y]= ADC( f,w ) %UNTITLED 进行sigma-delta转换 %
此处显示详细说明 y(1)=1;u(1)=0;n=100000;T=1; for m=2:n;
s(m)=exp(1i*pi*m*2*w/n)*exp(1i*f/180*pi);
u(m)=u(m-1)+(s(m)-y(m-1))/n;
if u(m)
y(m)=-1;
else
y(m)=1;
end end plot(y) end function [ f ] = phase( N,v,u,r,w) %UNTITLED3 此函数描述给定圆矩阵，各个接收器的相位差 %
此处显示详细说明
v为纬度角，u为经度角
for i=1:1:N; p(i)=2*pi*i/N; f(i)=w*r*cos(v)*cos(p(i)-u)/340*180/pi; while(f(i)>360)
f(i)=f(i)-360; end while(f(i)
f(i)=f(i)+360; end if f(i)>180;
f(i)=360-f(i); end end end function [o] = micGet( f,w ) %UNTITLED5 接收器接收信号计算相位差值 %
此处显示详细说明 y(1)=1;u(1)=0;y1(1)=1;u1(1)=0;
%%%%%%%%进行sigma-delta转换 y=ADC(0,w); y1=ADC(f,w); %%%%%%进行频域分析
y1fft=(real(fft(y1,100000))); yfft=(real(fft(y,100000))); if max(y1fft)>abs(real(min(y1fft)))
if max(yfft)>abs(real(min(yfft))) [y1max,tp]=max(y1fft); [ymax,tp2]=max(yfft);
else if max(yfft)
[y1max,tp]=max(y1fft);
[ymax,tp2]=min(yfft);
end
end
else if max(y1fft)
if max(yfft)>abs(real(min(yfft))) [y1max,tp]=min(y1fft); [ymax,tp2]=max(yfft);
else if max(yfft)
[y1max,tp]=min(y1fft);
[ymax,tp2]=min(yfft);
end
end
end end t1=real(y1max/ymax);%%%%%两个极值点的比值，为exp（j*2*pi*f)
o=acos(t1)/pi*180%%%%%得到f*360.单位为度 end
function Get( N,u,v,r,w ) %UNTITLED7 遍历查找声源位置 %
此处显示详细说明 f=phase(N,v,u,r,w);
for i=1:1:N; o(i)=micGet(f(i),w); end [x,y]=meshgrid(1:1:180); for i=1:1:180;
for j=1:1:180;
t=phase(N,j,i,r,w);
z(i,j)=(abs(t)-o)*(abs(t)-o)\';
z(i,j)=1/z(i,j);
end end mesh(x,y,z); o i=phase(N,v,u,r,w) end
人人范文网 m.inrrp.com.cn 手机版