esprit与music算法性能

③ WCDMA数据卡的关键技术

空时处理技术通过在空间和时间上联合进行信号处理可以非常有效地改善系统特性。随着第三代移动通信系统对空中接口标准的支持以及软件无线电的发展，空时处理技术必将融入自适应调制解调器中，从而达到优化系统设计的目的。采用空时处理的方法，系统的发送端或接收端使用多个天线，同时在空间和时间上处理信号，它所达到的效果是仅靠单个天线的单时间处理方法所不能实现的：可以在一个给定BER质量门限下，增加用户数；在小区给定的用户数下，改善BER特性；可以更有效地利用信号的发射功率等等。
1、空时处理方法
在单用户的情况下，空时处理方法的分类如图1所示。
图1
由于移动台一般不适于用多天线接收，在基站采用多个天线进行发射分集，可以使移动台的接收效果和移动台用多个接收天线时的效果相比拟，所以本文主要围绕基站的空时处理技术展开讨论。
2、波束成形技术
波束成形技术（Beamforming, BF）可分为自适应渣慧波束成形、固定波束和切换波束成形技术。固定波束即天线的方向图是固定的，把IS－95中的3个120°扇区分割即为固定波束。切换波束是对固定波束的扩展，将每个120°的扇区再分为多个更小的分区，每个分区有一固定波束，当用户在一扇区内移动时，切换波束机制可自动将波束切换到包含最强信号的分区，但切换波束机制的致命弱点是不能区分理想信号和干扰信号。自适应波束成形器可依据用户信号在空间传播的不同路径，最佳地形成方向图，在不同到达方向上给予不同的天线增益，实时地形成窄波束对准用户信号，而在其他方向尽量压低旁瓣，采用指向性接收，从而提高系统的容量。由于移动台的移动性以及散射环境，基站接收到的信号的到达方向是时变的，使用自适应波束成形器可以将频率相近但空间可分离的信号分离开，并跟踪这些信号，调整天线阵的加权值，使天线阵的波束指向理想信号的方向。
自适应波束成形的关键技术是如何较精确地获得信道参数呢？对于上行链路，根据形成波束所用的信息可以将波束成形技术分成以下3类。
(1)基于空间结构的BF
基于空间结构的BF如基于输入信号到达方向的BF（DOB），包括3类：基于最大信干噪比（SINR）的BF；基于最大似然（ML）准则的BF；基于最小均方误差（MMSE）准则的BF。多址干扰的抑制依赖于信号的到达方向如陪答（DOA），所以DOB中的一个重要部分是信号的DOA估计。DOA估计方法有离散付里叶变换、MVDR（Minimum Variance Distortionless Response）估计器、线性预测、最大包络法（MEM）、ML滤波器以及可变特征结构的方法，其中包乱举括MUSIC（Multiple Signal Classification）和ESPRIT法（Estimation of Signal Parameters via Rotational Invariance Technique）。
(2)基于训练序列的BF
基于训练序列的BF即时间参考BF（TRB），适用于多径丰富且信道特性连续变化的环境，根据算法可以分为块自适应算法（BAA）和采样自适应算法（SAA）两类。BAA算法包括特征滤波器（EF）法、Stanford法、最大比合并（MRC）法和第一维纳滤波器解（FWFS）、第二维纳滤波器解（SWFS）。SAA算法包括最小均方（LMS）算法、归一最小均方（NLMS）算法、递归最小平方（RLS）算法和共轭梯度法（CGM）。TRB技术要求同步精确，当时延扩展小时可以得到较好的性能。
(3)基于信号结构的BF（SSBF）
基于信号结构的BF（SSBF）即利用接收信号的时间或空间结构和特性来构造BF，可利用SSBF需要存储例如恒包络调制信号的恒模（CM）特性、信号的周期平稳性或数字调制信号的FA（Finite Alphabet）特性等知识，这种BF方法可以应用于不同的传播条件，但需要考虑收敛性和捕获问题。
对于下行链路而言，不同的复用方式可采用不同的解决方法：TDD方式，由于上下行链路采用相同的频率，在保证信道参数在相邻的上下行数据帧中几乎没有变化的情况下可以直接利用上行估计得到的信道参数，但这只适用于慢速移动的系统；FDD方式，由于上下行链路的频率间隔一般都大于相关带宽，因此上下行的瞬时信道几乎是不相关的，此时采用反馈信道是最好的方法。
需要强调的一点是发送机的波束成形技术和接收机的波束成形技术是截然不同的，接收波束成形可在每个接收机独立实现而不会影响其他链路，而发送波束成形会改变对其他所有接收机的干扰，所以要在整个网络内部联合使用发送波束成形技术。
3、接收分集
由于CDMA系统通常有较多的多址干扰分量，而天线阵可以去除M－1个（M为天线数）干扰的特性并不能明显地改善接收机 的SINR，所以在一般情况下，更好的方法是利用接收分集的方法，估计接收信号的形式，并确定匹配滤波器的加权系数。接收分集技术中的分集天线其实是空间域内的分集合并器，而不是BF。对于宽带CDMA信号，信号带宽一般大于信道相干带宽，所以在时间域采用RAKE接收机，将信号在空间/时间上利用各种合并准则进行合并，这就是所谓的2D－RAKE接收机。一般的合并方式有：选择合并（SC）即选择具有最大信号功率的多径；最大比合并（MRC）即每一路有一加权，根据各支路信噪比（SNR）来分配加权的权重，SNR大的支路权重大，SNR小的支路权重小。当每个分离多径上的干扰不相关时，MRC方法可使合并信号的SINR最大；等增益合并（EGC）即选择每一路的加权值都相等；Wiener滤波（OPT）即无论多径之间的干扰是否相关，均可抑制干扰并使合并器输出端的SINR最大，因此Wiener滤波的方法要好于最大比合并法，又称为优化合并。
在空间和时间上利用不同的合并准则可以对系统起到不同的改善效果，理论证明，在理想功率控制和理想信道估计的条件下，空时联合域优化合并方式对系统性能的改善最好。
4、发送分集技术
当发送方不能获得信道参数时，空时发送分集可改善前向链路性能，这种机制是将发送天线的空间分集转化为接收机可以利用的其他形式的分集，如延迟发送分集和空时编码技术。空时编码技术是同时从空间和时间域考虑设计码字，它的基本原理是在多个天线上同时发送信息比特流所产生的向量，利用发送天线所发送序列的正交性，用两个发送天线、一个接收天线所获得的分集增益与一个发送天线、两个接收天线的MRC接收机的一样。
根据是否需要从接收机到发射机的反馈电路，发送分集技术可以分为开环和闭环两种类型，前者发射机不需要任何信道方面的知识。开环发送分集方式有空时发送分集（STTD）、正交发送分集（OTD）、时间切换发送分集（TSTD）、延迟发送分集（DTD）以及分层的空时处理和空时栅格编码；闭环发送分集方式有选择发送分集（STD）。发送分集各方式具体如下。
(1)正交发送分集（OTD）
经过编码和交织后的数据分成两个不同的子流在两个不同的天线上同时发送。为保证正交性，这两个子流所用的Walsh码是不同的。
(2)时间切换发送分集（TSTD）
在某一时刻每个用户只使用一个天线，使用伪随机码机制在两个天线之间切换。
(3)选择发送分集（STD）
由于在TSTD方式中，瞬时使用的发送天线并不一定能在接收端得到最大的信噪比，所以使用一个反馈电路来选择能提供使接收端得到最大信噪比的天线。
(4)空时发送分集（STTD）
空时发送分集是按如图2所示的方法将数据编码之后在两个天线上发送出去。
(5)延迟发送分集（DTD）
用多个天线在不同时刻发送同一原始数据信号的多个复本，人为地产生多径。
(6)分层空时结构(Bell LAyered Space－Time architecture，BLAST）
首先将原始信息比特分解成n个并行的数据流(称为层)，送入不同的编码器，再将编码器的输出调制以后使用相同的Walsh码通过不同的天线发送出去。接收机侧使用一个BF（迫零或MMSE准则）来分离不同的编码数据流，然后将数据送入不同的解码器，解码器的输出再重新组合建立原始的信息比特流。由于在波束成形处理中，MMSE和迫零方法都没有充分利用接收机天线阵的分集潜力，所以提出了改进方案将接收处理也进行分级。即首先使用ViterbiMLSE算法译出最强的信号，然后将该强信号从接收的天线信号中去除后再检测第二强的信号，如此反复直到检测出最弱的信号。
该机制中，层到天线的映射并不是固定的，而是每np个码符号之后周期性地改变，如图3所示。这种映射关系保证了这些数据流最大可能地在不同的天线上被发送出去。
(7)空时栅格编码
根据秩准则和行列式准则设计码字，使设计出的码字得到最大分集增益和编码增益。以四进制相移键控（QPSK）四状态空时栅格编码为例，假定使用两根天线发射，则星座图和格形如图4所示。
最右边的元素编号S1S2的涵义是：从第一根天线发射出去的字符为S1，从第二根天线发射的字符为S2。
5、结束语
空时处理技术已显示出非常诱人的发展前景，第三代移动通信标准中也支持空时处理技术，标准的出台为我们继续研究物理可实现的空时处理技术提供了可能性，但将此技术实用化还存在许多亟待解决的方法和技术问题，有待于我们进一步研究。

④ 什么叫窄带，宽带，单音信号

信号带宽远小于中心频率的是窄带信号，带宽能和中心频率相比拟或着是远大于中心频率的信号是宽带信号。
窄带信号的功率集中在中心频率附近，两者的功率谱密度和频谱密度图灶启有很大的差距。处理方法也有很大差距。
宽带信号与窄带信号是相对的，不满足窄带条件的信号就称为宽带信号。目前，窄带信号的定义也不尽相同。若信号带宽为B
，时宽为T，中心频率为f0，则窄带信号的定义有：
定义1：
B
<<
f
0，即相对带宽B/f
0<<
1
,一般B/f
0<
0.1。
定义2：
2v/c<<1/TB,其中v
是阵列与目标的相对径向运动速度，c
是信号在介质中的传播速度。
定义3：(
N－1)
d/c<<1/B，其中N
是阵元数目，d
是阵元间距。
定义4：2
Bτθ
+1≈1
,其中τθ是整个阵列以及延迟线的没租延时之和。
定义5：该信号空间协方差矩阵在没有噪声时的第二个特征值小于噪声功率。
“定义1”是对窄带信号的直观理解,同时也是窄带实信号可有效表示为其复解析形式的充分条件,在很多文献中均以该定义来区分信号是宽带信号还是窄带信号。
“定义2”是指在存在相对运动的系统中,在信号的持续时间T
内,相对于信号的距离分辨率,目标没有明显的位移,此时信号可视为窄带信号,否则信号就是宽带信号。
“定义3”是指在阵列信号处理中,如果信号带宽的倒数远远大于信号掠过阵列孔径的最大传播时间,就称为窄带信号,否则为宽带信号。
“定义4”和“定义5”是从阵列采样数据自相关矩阵特征值的角度来定义窄带信号的,窄带情况下,阵列采样数据自相关矩阵的大特征值个数等于信号个数。
可见，枯辩兆窄带信号定义的非一致性决定了宽带信号定义的非绝对性，不同的处理场合，应使用不同的定义确定信号是否是宽带的。
至于某些窄带信号的分析结果与算法不适合宽带信号，原因也是有差别的，不能给出统一的解释。例如，在DOA估计中，子空间方法中的MUSIC、ESPRIT
方法以及它们的改进算法，都是基于空间信号是窄带这一假设推导出的。因为当信号是窄带时，信号的延时近似等于信号的相移。当接收到的信号为宽带信号时，近似条件不满足，算法性能会严重下降。
虽然宽带信号没有官方或某一组织的严格定义,但对于超宽带信号,
FCC
(
Federal
Communica2tions
Commission)
给出了严格的定义.“超宽带”的定义是:信号的相对带宽(信号频谱的带宽与其中心频率之比)
大于等于20
%
,或者绝对带宽大于等于500MHz。该定义没有界定信号的时域波形特征。

⑤ 从经典谱估计到现代谱估计

谐波分析最早可追溯到古代对时间的研究。18世纪伯努利（Bernoulli）、欧拉（Euler）和拉格朗日（Lagrange）等人对波动方程及其正弦解进行了研究，19世纪初叶，傅立叶（Fourier）证明了在有限时间段上定义的任何函数都可以用正弦和余弦分量的无限谐波的总和来表示。1898年舒斯特（Schuster）以傅立叶分析为基础来拟合待分析信号，研究太阳黑子数的周期变化，并提出了周期图的概念。1930年维纳（Wiener）发表了经典性论文《广义谐波分析》，对平稳随机过程的自相关函数和功率谱密度作了精确的定义，证明了二者之间存在着傅氏变换（以下简称傅氏变换）的关系，从而为功率谱分析奠定了坚实的统计学基础。由于1934年辛钦（Khintchine）也独立地证明了自相关函数和功率谱之间的傅氏变换关系，即维纳-辛钦（Wiener-Khintchine）定理。根据这个定理（详见第一章），平稳离散随机信号x（n）的自相关函数r_xx（m）

r_xx（m）=E［x（m+n）x^*（n）］ （4-1）

与功率谱P_xx（e^jω）之间构成一傅氏变换对，即

地球物理信息处理基础

若x（n）还是各态遍历性的，则其自相关函数可由它的一个采样时间序渗销列用时间平均的方法求出，即

地球物理信息处理基础

在大多数应用中x（n）是实信号，于是上式可写成

地球物理信息处理基础

实际上，一般只能在时域观测到随机信号的有限个采样值（例如N个值），可表示为

x_N（n）=｛x（0），x（1），…，x（N-1）｝=｛x（n），n=0，1，…，N-1｝

其自相关函数只能由这N个采样数据进行估计，常用有偏估计

地球物理信息处理基础

这是一种渐近一致估计，称之为采样自相关函数。

用采样自相关函数的傅氏变换作为功率谱的估计，这种方法是布莱克曼（R.Blackman）和杜基（J.Tukey）在1958年提出来的，称为功率谱估计的自相关法（简称BT法）。此方法需要先求出有限个观测数据估计自相关函数，然后再根据式（4-2）计算出功率谱。在快速傅氏变换（FFT）算法提出之前，这是一种最流行的功率谱估计方法。

1965年库利（Cooley）和杜基（Tukey）完善了着名的FFT算法，使计算傅氏变换的速度提高了两个数量级，运算量显着降低，这样DFT变换很快在各领域，特别是在工程实践中得到了孝运广泛应用。由式（4-5）知，

为x（n）与x（-n）的卷积运算，因为

地球物理信息处理基础

若x（n）的傅氏变换为X（e^jω），则x（-n）的傅氏变换等于X^*（e^jω）。对式（4-5）两端取傅氏变换，得到

地球物理信息处理基础

这表明：通过对随机数据直接进行离散傅氏变换，然后取其幅值的平方，再对多样本进行此种运算并取平均值作为功率谱的估计，即舒斯特的周期图，这种谱估计受到了人们的普遍重视巧喊梁，因为它不需要计算自相关函数，而直接计算功率谱。

周期图和自相关法以及它们的改进方法称为功率谱估计的经典方法，周期图和自相关法是经典功率谱估计的两个基本方法。由于FFT的出现，周期图和自相关法往往被结合起来使用，其步骤如下：

（1）对x_N（n）补N个零，求

；

（2）由

作傅氏变换，得

，这时｜m｜≤M=N-1；

（3）对

加窗函数v（m），这时｜m｜≤M＜＜N-1，得

；

（4）利用

，求

的傅氏变换，即

地球物理信息处理基础

由周期图法得到的功率谱

，其估计方差并不随样本长度的增加而趋于零，

，出人意料的是，不管数据记录有多长，周期图和自相关法得到的估计都不是功率谱的良好估计。事实上，随着记录长度增加，这两种估计的随机起伏反而会更加严重！此外，它们存在着以下两个难以克服的固有缺点。

（1）频率分辨率（区分两个邻近频率分量的能力）不高。因为它们的频率分辨率（赫兹）反比于数据记录长度（秒）（即Δf=k/T_p=k/NT，k为常数，T_p=NT为数据的记录长度，T为采样周期），而实际应用中一般不可能获得很长的数据记录，即观察到的数据只能是有限个，而观察不到的数据被认为是0。这样，如果只有N个观测数据，而对于N以外的数据，信号仍有较强的相关性，那么估计出的功率谱就会出现很大的偏差。

（2）对于有限的观测数据，相当于将信号在时域内乘以矩形窗函数，因而在频域内则相当于使真正的功率谱与sinc函数进行卷积，由于sinc函数不同于δ 函数，它有主瓣和旁瓣，这样使卷积后的功率谱不同于真正的功率谱。sinc函数的主瓣不是无限窄的，引起功率谱向附近频域扩展，造成谱的模糊，降低谱的分辨率；同时，由于sinc函数的旁瓣存在，导致能量向旁瓣中“泄漏”（称之为旁瓣泄漏），即引起频谱间的干扰，信号强的功率谱旁瓣影响信号弱的功率谱检测，严重时，会使主瓣产生很大失真，检测不出弱信号，或者把旁瓣误认为是信号，造成假信号。为了对经典功率谱估计进行改进，可以采用各种不同的窗函数，但其结果都是以增加主瓣宽度来换取旁瓣的压低，因此功率谱分辨率低是经典功率谱估计的致命缺点。

为了克服以上缺点，人们曾做过长期努力，提出了平均、加窗平滑等办法，在一定程度上改善了经典功率谱估计的性能。实践证明，对于长数据记录来说，以傅氏变换为基础的经典功率谱估计方法，的确是比较实用的。但是，经典方法始终无法根本解决频率分辨率和功率谱估计稳定性之间的矛盾，特别是在数据记录很短的情况下，这一矛盾显得尤为突出。这就促进了现代功率谱估计方法研究的展开。

现代功率谱估计方法主要是以随机过程（Stochastic Process）的参数模型（Parameter Model）为基础的，称之为参数模型方法。虽然说现代功率谱估计技术的研究和应用主要起始于60年代，但实际上，时间序列模型在非工程领域早已被采用，如Yule在1927年、Walker在1931年都曾使用过自回归模型预测描述经济的时间序列的发展趋势，而Prony则早在1795年就曾采用指数模型去拟合在气体化学实验中获得的数据。在统计学和数值分析领域中，人们也曾采用过模型方法。

现代功率谱估计的提出主要是针对经典功率谱估计（周期图和自相关法）的分辨率和方差性能不好的问题而提出的。1967年Burg在地震学研究中受到线性预测滤波的启发，提出了最大熵谱估计方法，在提高分辨率方面作了最有意义的探索。1968年Parzen正式提出了自回归谱估计方法。1971年Van der Bos证明了一维最大熵谱估计与自回归谱估计等效。1972年出现的谱估计的Prony方法在数学上与自回归方法有某些类似。目前以自回归滑动平均模型为基础的谱估计已经比自回归模型谱估计具有更高的频率分辨率和更好的性能。1973年Pisarenko提出的谐波分解方法提供了可靠的频率估计方法。1981年Schmidt提出了谱估计的多信号分类（MUSIC）算法等。因此，现代功率谱分析主要有ARMA谱分析、最大似然、熵谱估计和特征分解四种方法。ARMA谱分析是一种建模方法，即通过平稳线性信号过程建立模型来估计功率谱密度；熵谱估计包括最大熵谱和最小交叉法；特征分解也叫特征构造法和子空间法，包括Pisarenko谐波分解法、Prony法、MUSIC法和ESPRIT法（用旋转不变技术估计参数方法）。

现代功率谱估计研究仍侧重于一维功率谱分析，而且大部分是建立在二阶矩（相关函数、方差、功率谱密度）基础上的。但由于功率谱密度是频率的实函数，缺少相位信息，因此，建立在高阶谱基础上的谱估计方法正引起人们的注意，特别是双谱估计和三谱估计的研究受到了高度的重视。其它如多维谱估计、多通道谱估计等的研究也正在发展中。人们希望这些新方法能更多地在提取信息、估计相位和描述非线性等方面获得应用。

⑥ 高分求一段汉译英，用软件直接翻译的勿扰。

With the development of science and technology, now in many areas of need to know the exact direction of useful signal, the frequency of parameters, the precise measurement of the signal to put a very high demand. Spectrum space technology in many areas has been widely used, and array signal processing plays a very important role.
Often, the signal space is measured up to one direction Wei Boda (DOA) estimates. In practical applications, often at the same time to signal the number of parameters of the estimates, such as signal DOA, pitch angle and frequency. This paper first introced for the two-dimensional space spectrum of the estimated MUSIC algorithm, ESPRIT algorithm for the basic idea and principle of algorithm, and then focus on the two algorithms based on the position - two-dimensional frequency spectrum estimation. Then this paper were used MATLAB language MUSIC algorithm and ESRPIT algorithm estimated two-dimensional and computer simulation results are analyzed and compared.

This paper first introced the concept of radio frequency identification technology, the rise and development. Radio frequency identification system with the basic structure and principles, described the antenna transceiver and the working principle of the process. Subsequent analysis of the transmitting power, receive power chip, chip power consumption and receiving antennas to receive the echo of the relationship between power and to examine in detail the chip to receive power and receiving antennas to receive the echo power as a range of changes, as well as the decision The size of these power factors, the actual design of radio frequency identification hardware to provide reference data

参考

⑦ 在阵列信号处理中，为什么噪声子空间与导向矢量是相互正交的

阵列信号处理经过多年来的发展，时域、空域上所蕴含的信息已得到充分的挖掘。近年来，新的研究工作在极化上得到开展。作为一种既能感知极化信息又能感知时空域信息的接收装置，电磁矢量传感器阵列接收到的数据具有高维结构。然而，传统的方法处理高维数据是将其转化为矩阵数据，这使得数据的高维信息未得到有效利用。

利用张量这一处理高维数据的专用工具，研究了非完全极化波信号的BTD（Block
term
decompositions）模型、基于Tucker分解的参数估计方法、基于CP（Canoical
Decomposition）分解的参数估计方法和基于BTD分解的DOA估计方法。本文围绕如...
展开
阵列信号处理经过多年来的发展，时域、空域上所蕴含的信息已得到充分的挖掘。近年来，新的研究工作在极化上得到开展。作为一种既能感知极化信息又能感知时空域信息的接收装置，电磁矢量传感器阵列接收到的数据具有高维结构。然而，传统的方法处理高维数据是将其转化为矩阵数据，这使得数据的高维信息未得到有效利用。
本文利用张量这一处理高维数据的专用工具，研究了非完全极化波信号的BTD（Block
term
decompositions）模型、基于Tucker分解的参数估计方法、基于CP（Canoical
Decomposition）分解的参数估计方法和基于BTD分解的DOA估计方法。本文围绕如何充分利用电磁矢量传感器阵列中蕴含的高维信息，
展开了以下几个部分工作：

1、改进了基于张量的Tucker分解的MUSIC算法。通过张量分解改善了估计的噪声子空间与真实信号导向矢量的之间的正交性，提升了对电磁矢量传感器阵列参数估计的分辨力。相对于矩阵方法对噪声子空间的估计，张量法利用到了数据的各个维度的信息，能从被噪声覆盖的数据中还原出更接近真实值的数据。仿真表明，利用张量分解估计的噪声子空间进行MUSIC参数估计能够提升来波信号的分辨力。这部分说明了张量各维度之间的整体性，即张量数据的高维特征，是一种值得利用的信息。

2、改进了基于张量的CP分解对阵列进行盲估计的方法。相比于Tucker分解的估计方法，CP分解不仅能反映数据各维度之间的整体性，更重要的是还可以盲估计出各维度上的组成成分。利用这种盲估计特性，可以进一步将对应组成成分具有的结构特点作为约束去影响CP分解的每一步迭代。即这种盲估计特性为充分利用阵列自身具有的结构特点创造了条件。另外，本文提出了基于参数化的加入结构约束的方法。即在每一步迭代中，利用盲估计出来的组成成分估算参数，并根据参数和信号模型生成满足模型特征的组成成分。仿真表明，加入了结构信息的约束之后，参数估计误差有所减小。
3、研究了基于
BTD分解的非完全极化波参数估计方法。对于非完全极化波的电磁矢量传感器阵列接收模型，同一方向的来波总是可以等效为两路不相干的信号被一个6×2矩阵导向。这时候，CP分解唯一性所需的条件不满足，本文首次采用了符合这种阵列模型的BTD分解。BTD分解可以盲估计出空域阵列导向矢量和极化部分的一个列空间。另外，根据电磁场中的坡印廷定理，本文研究了从这个列空间中求解方向向量的方法。最后，融合空域阵列导向矢量和方向向量的估计得到来波方向。最后与ESPRIT算法进行了仿真对比。结果表明，对于非完全极化波的接收模型，基于BTD分解的方法能够得到更精确的估计结果。

⑧ DOA估计算法

学号：20000300055

姓名：王铎澎

嵌牛导读：文章对DOA算法进行了简单的介绍。

嵌牛正文：https://blog.csdn.net/zhangziju/article/details/100730081?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160689878119725222413438%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=160689878119725222413438&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~_landing_v2~default-1-100730081.pc_first_rank_v2_rank_v28&utm_term=Musicsuanfa&spm=1018.2118.3001.4449

DOA估计算法

DOA（Direction Of Arrival）波达方向定位技术主要有ARMA谱分析、最大似然法、熵谱分析法和特征分解法，特征分解法主要有MUSIC算法、ESPRIT算法WSF算法等。

MUSIC (Multiple Signal Classification)算法，即多信号分类算法，由Schmidt等人于1979年提出。MUSIC算法是一种基于子空间分解的算法，它利用信号子空间和噪声子空间的正交性，构建空间谱函数，通过谱峰搜索，估计信号的参数。对于声源定位来说，需要估计信号的DOA。MUSIC算法对DOA的估计有很高的分辨率，且对麦克风阵列的形状没有特殊要求，因此应用十分广泛。

运用矩阵的定义，可得到更为简洁的表达式：

X = A S + N X=AS+NX=AS+N

式中

X = [ x 1 ( t ) , x 2 ( t ) , . . . x M ( t ) ] T X=[x_1(t),x_2(t),...x_M(t)]^TX=[x1​(t),x2​(t),...xM​(t)]T，

S = [ S 1 ( t ) , S 2 ( t ) , . . . S D ( t ) ] T S=[S_1(t),S_2(t),...S_D(t)]^TS=[S1​(t),S2​(t),...SD​(t)]T，

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^TA=[a(θ1​),a(θ2​),...a(θD​)]T，

N = [ n 1 ( t ) , n 2 ( t ) , . . . n M ( t ) ] T N=[n_1(t),n_2(t),...n_M(t)]^TN=[n1​(t),n2​(t),...nM​(t)]T。

X XX为阵元的输出，A AA为方向响应向量，S SS是入射信号，N NN表示阵列噪声。

其中 φ k = 2 π d λ s i n θ k \varphi_k=\frac{2\pi d}{\lambda}sin\theta_kφk​=λ2πd​sinθk​有

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[

1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD

\right]A=[a(θ1​),a(θ2​),...a(θD​)]T=⎣⎢⎢⎢⎡​1e−jφ1​⋮e−j(M−1)φ1​​1e−jφ2​⋮e−j(M−1)φ2​​⋯⋯⋱⋯​1e−jφD​⋮e−j(M−1)φD​​⎦⎥⎥⎥⎤​

对x m ( t ) x_m(t)xm​(t)进行N点采样，要处理的问题就变成了通过输出信号x m ( t ) x_m(t)xm​(t)的采样{ x m ( i ) = 1 , 2 , . . . , M } \{ x_m (i)=1,2,...,M\}{xm​(i)=1,2,...,M}估计信号源的波达方向角θ 1 , θ 2 . . . θ D \theta_1,\theta_2...\theta_Dθ1​,θ2​...θD​,由此可以很自然的将阵列信号看作是噪声干扰的若干空间谐波的叠加，从而将波达方向估计问题与谱估计联系起来。

对阵列输出X做相关处理，得到其协方差矩阵

R x = E [ X X H ] R_x=E[XX^H]Rx​=E[XXH]

其中H HH表示矩阵的共轭转置。

根据已假设信号与噪声互不相关、噪声为零均值白噪声，因此可得到：

R x = 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_x=E[(AS+N)(AS+N)^H] =AE[SS^H]A^H+E[NN^H]=AR_SA^H+R_NRx​=E[(AS+N)(AS+N)H]=AE[SSH]AH+E[NNH]=ARS​AH+RN​

其中R s = E [ S S H ] R_s=E[SS^H]Rs​=E[SSH]称为信号相关矩阵

R N = σ 2 I R_N=\sigma^2IRN​=σ2I是噪声相关阵

σ 2 \sigma^2σ2是噪声功率

I II是M × M M\times MM×M阶的单位矩阵

在实际应用中通常无法直接得到R x R_xRx​，能使用的只有样本的协方差矩阵：

R x ^ = 1 N ∑ i = 1 N X ( i ) X H ( i ) \hat{R_x}=\frac{1}{N} \sum_{i=1}^{N}X(i)X^H (i)Rx​^​=N1​∑i=1N​X(i)XH(i)，R x ^ \hat{R_x}Rx​^​是R x R_xRx​的最大似然估计。

当采样数N → ∞ N\to\inftyN→∞，他们是一致的，但实际情况将由于样本数有限而造成误差。根据矩阵特征分解的理论，可对阵列协方差矩阵进行特征分解，首先考虑理想情况，即无噪声的情况：R x = A R s A H R_x=AR_sA^HRx​=ARs​AH，对均匀线阵，矩阵A由

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[

1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD

\right]A=[a(θ1​),a(θ2​),...a(θD​)]T=⎣⎢⎢⎢⎡​1e−jφ1​⋮e−j(M−1)φ1​​1e−jφ2​⋮e−j(M−1)φ2​​⋯⋯⋱⋯​1e−jφD​⋮e−j(M−1)φD​​⎦⎥⎥⎥⎤​

所定义的范德蒙德矩阵，只要满足θ i ≠ θ j , i ≠ j \theta_i\neq \theta_j,i\neq jθi​​=θj​,i​=j，则他的各列相互独立。

若R s R_sRs​为非奇异矩阵R a n k ( R s ) = D Rank(R_s)=DRank(Rs​)=D，各信号源两两不相干，且M > D M>DM>D，则r a n d ( A R s A H ) = D rand(AR_sA^H)=Drand(ARs​AH)=D，

由于R x = E [ X X H ] R_x=E[XX^H]Rx​=E[XXH]，有：

R s H = R x R_s^H=R_xRsH​=Rx​

即R s R_sRs​为Hermite矩阵，它的特性是都是实数，又由于R s R_sRs​为正定的，因此A R s A … … H AR_sA……HARs​A……H为半正定的，它有D个正特征值和M − D M-DM−D个零特征值。

再考虑有噪声存在的情况

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx​=ARs​AH+σ2I

由于σ 2 > 0 \sigma^2>0σ2>0，R x R_xRx​为满秩阵，所以R x R_xRx​有M个正实特征值λ 1 , λ 2 . . . λ M \lambda_1,\lambda_2...\lambda_Mλ1​,λ2​...λM​

分别对应于M个特征向量v 1 , v 2 . . . v M v_1,v_2...v_Mv1​,v2​...vM​。又由于R x R_xRx​为Hermite矩阵，所以各特征向量是正交的，即：v i H v j = 0 , i ≠ j v_i^Hv_j=0,i\neq jviH​vj​=0,i​=j与信号有关的特征值只有D个,分别等于矩阵A R s A H AR_sA^HARs​AH的各特征值与σ 2 \sigma^2σ2之和，其余M − D M-DM−D个特征值为σ 2 \sigma^2σ2，即σ 2 \sigma^2σ2为R RR的最小特征值，它是M − D M-DM−D维的，对应的特征向量v i , i = 1 , 2 , . . . , M v_i,i=1,2,...,Mvi​,i=1,2,...,M中，也有D个是与信号有关的，另外M − D M-DM−D个是与噪声有关的，可利用特征分解的性质求出信号源的波达方向θ k \theta_kθk​。

MUSIC算法的原理及实现

通过对协方差矩阵的特征值分解，可得到如下结论：

将矩阵R x R_xRx​的特征值进行从小到大的排序，即λ 1 ≥ λ 2 ≥ . . . ≥ λ M > 0 \lambda_1 \geq \lambda_2\geq...\geq\lambda_M>0λ1​≥λ2​≥...≥λM​>0，其中D个较大的特征值对应于信号，M − D M-DM−D个较小的特征值对应于噪声。

矩阵R x R_xRx​的属于这些特征值的特征向量也分别对应于各个信号和噪声，因此可把R x R_xRx​的特征值（特征向量）划分为信号特征（特征向量）与噪声特征（特征向量）。

设λ i \lambda_iλi​为R x R_xRx​的第i ii个特征值，v i v_ivi​是与λ i \lambda_iλi​个相对应的特征向量，有：

R x v i = λ i v i R_xv_i=\lambda_iv_iRx​vi​=λi​vi​

再设λ i = σ 2 \lambda_i=\sigma^2λi​=σ2是R x R_xRx​的最小特征值R x v i = σ 2 v i i = D + 1 , D + 2... M R_xv_i=\sigma^2v_i i=D+1,D+2...MRx​vi​=σ2vi​i=D+1,D+2...M，

将R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx​=ARs​AH+σ2I代入可得σ 2 v i = ( A R s A H + σ 2 I ) v i \sigma^2v_i=(AR_sA^H+\sigma^2I)v_iσ2vi​=(ARs​AH+σ2I)vi​，

将其右边展开与左边比较得：

A R s A H v i = 0 AR_sA^Hv_i=0ARs​AHvi​=0

因A H A A^HAAHA是D ∗ D D*DD∗D维的满秩矩阵，( A H A ) − 1 (A^HA)^{-1}(AHA)−1存在；

而R s − 1 R_s^{-1}Rs−1​同样存在，则上式两边同乘以R s − 1 ( A H A ) − 1 A H R_s^{-1}(A^HA)^{-1}A^HRs−1​(AHA)−1AH，

有：

R s − 1 ( A H A ) − 1 A H A R s A H v i = 0 R_s^{-1}(A^HA)^{-1}A^HAR_sA^Hv_i=0Rs−1​(AHA)−1AHARs​AHvi​=0

于是有

A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,...,MAHvi​=0，i=D+1,D+2,...,M

上式表明：噪声特征值所对应的特征向量（称为噪声特征向量）v i v_ivi​，与矩阵A AA的列向量正交，而A AA的各列是与信号源的方向相对应的，这就是利用噪声特征向量求解信号源方向的出发点。

用各噪声特征向量为例，构造一个噪声矩阵E n E_nEn​:

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},...v_{M}]En​=[vD+1​,vD+2​,...vM​]

定义空间谱P m u ( θ ) P_{mu}(\theta)Pmu​(θ):

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) = 1 ∥ E n H a ( θ ) ∥ 2 P_{mu}(\theta)=\frac{1}{{a^H}(\theta)}E_nE_n^Ha(\theta)=\frac{1}{\Vert E_n^Ha(\theta)\Vert^2}Pmu​(θ)=aH(θ)1​En​EnH​a(θ)=∥EnH​a(θ)∥21​

该式中分母是信号向量和噪声矩阵的内积，当a ( θ ) a(\theta)a(θ)和E n E_nEn​的各列正交时，该分母为零，但由于噪声的存在，它实际上为一最小值，因此P m u ( θ ) P_{mu}(\theta)Pmu​(θ)有一尖峰值，由该式，使θ \thetaθ变化，通过寻找波峰来估计到达角。

MUSIC算法实现的步骤

1.根据N个接收信号矢量得到下面协方差矩阵的估计值：

R x = 1 N ∑ i = 1 N X ( i ) X H ( i ) R_x=\frac{1}{N} \sum_{i=1}^NX(i)X^H(i)Rx​=N1​∑i=1N​X(i)XH(i)

对上面得到的协方差矩阵进行特征分解

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx​=ARs​AH+σ2I

2.按特征值的大小排序 将与信号个数D DD相等的特征值和对应的特征向量看做信号部分空间，将剩下的M − D M-DM−D个特征值和特征向量看做噪声部分空间，得到噪声矩阵E n E_nEn​:

A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,...,MAHvi​=0，i=D+1,D+2,...,M

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},...v_{M}]En​=[vD+1​,vD+2​,...vM​]

3.使θ \thetaθ变化 ，按照式

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) P_{mu}(\theta)=\frac{1}{{a^H}(\theta)E_nE_n^Ha(\theta)}Pmu​(θ)=aH(θ)En​EnH​a(θ)1​

来计算谱函数，通过寻求峰值来得到波达方向的估计值。

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;            % 阵元间距

d=0:dd:(N-1)*dd;

A=exp(-1i*2*pi*d.'*sin(theta*derad));  %方向矢量

%%%%构建信号模型%%%%%

S=randn(M,K);            %信源信号，入射信号

X=A*S;                    %构造接收信号

X1=awgn(X,snr,'measured'); %将白色高斯噪声添加到信号中

% 计算协方差矩阵

Rxx=X1*X1'/K;

% 特征值分解

[EV,D]=eig(Rxx);                  %特征值分解

EVA=diag(D)';                      %将特征值矩阵对角线提取并转为一行

[EVA,I]=sort(EVA);                %将特征值排序 从小到大

EV=fliplr(EV(:,I));                % 对应特征矢量排序

% 遍历每个角度，计算空间谱

for iang = 1:361

    angle(iang)=(iang-181)/2;

    phim=derad*angle(iang);

    a=exp(-1i*2*pi*d*sin(phim)).';

    En=EV(:,M+1:N);                  % 取矩阵的第M+1到N列组成噪声子空间

    Pmusic(iang)=1/(a'*En*En'*a);

end

Pmusic=abs(Pmusic);

Pmmax=max(Pmusic)

Pmusic=10*log10(Pmusic/Pmmax);            % 归一化处理

h=plot(angle,Pmusic);

set(h,'Linewidth',2);

xlabel('入射角/(degree)');

ylabel('空间谱/(dB)');

set(gca, 'XTick',[-90:30:90]);

grid on;

实现结果