esprit算法研究.doc

1、课程设计报告实验名称：ESPRIT 算法研究实验日期：姓 名：学 号：哈尔滨工业大学（威海）一、 设计任务实现空间谱估计算法，并考察算法性能。二、 方案设计1） 由均匀线阵形式，确定阵列的导向矢量；2） 由阵列导向矢量，对接收信号进行建模仿真；3） 由 ESPRIT 算法实现信号 DOA 估计；4） 考察算法性能与信噪比，采样率，观测时间等参数的关系。三、 设计原理3.1 空间谱估计数学模型空间谱估计就是利用空间阵列实现空间信号的参数估计的一项专门技术。整个空间谱估计系统应该由三部分组成：空间信号入射、空间阵列接收及参数估计。相应地可分为三个空间，即目标空间、观察空间及估计空间，也就是说空间谱

2、估计系统由这三个空间组成，其框图见图 1。图 1 空间谱估计的系统结构对于上述的系统结构，作以下几点说明。(1)目标空间是一个由信号源的参数与复杂环境参数张成的空间。对于空间谱估计系统，就是利用特定的一些方法从这个复杂的目标空间中估计出信号的未知参数。(2)观察空间是利用空间按一定方式排列的阵元，来接收目标空间的辐射信号。由于环境的复杂性，所以接收数据中包括信号特征(方位、距离、极化等)和空间环境特征(噪声、杂波、干扰等)。另外由于空间阵元的影响，接收数据中同样也含有空间阵列的某些特征(互耦、通道不一致、频带不一致等)。这里的观察空间是一个多维空间，即系统的接收数据是由多个通道组成，而传统的时

3、域处理方法通常只有一个通道。特别需要指出的是：通道与阵元并不是一一对应，通道是由空间的一个、几个或所有阵元合成的(可用加权或不加权)，当然空间某个特定的阵元可包含在不同的通道内。(3)估计空间是利用空间谱估计技术(包括阵列信号处理中的一些技术，如阵列校正、空域滤波等技术)从复杂的观察数据中提取信号的特征参数。从系统框图中可以清晰的看出，估计空间相当于是对目标空间的一个重构过程，这个重构的精度由众多因素决定，如环境的复杂性、空间阵元间的互耦、通道不一致、频带不一致等。3.2 阵列信号处理首先，考虑N个远场的窄带信号入射到空间某阵列上，阵列天线由M个阵元组成，这里假设阵元数等于通道数，即各阵元接收

4、到信号后经过各自的传输信道送到处理器，也就是说处理器接收来自M个通道的数据。 )(0()tjiietus(3.2-1)()(0ttjii式中， 是接受信号的幅度， 是接收信号的相位， 是接收信号的频率。在)(tui )(t 0窄带远场信号源的假设下，有(3.2-2)(ttuii根据式(3.2-1)和式(3.2-2)，显然有下式成立：(3.2-3)0)(jii etsts则可以得到第L个阵元接收信号为(3.2-4)Ni lliill tntsgtx1)()(Ml,21式中， 为第L个阵元对第i个信号的增益， 表示第 L个阵元在t时刻的噪声，lgl表示第i个信号到达第L个阵元时相对参考阵元的时延。

5、l将 M 个阵元在特定时刻接收的信号排列成一个列矢量，可得(3.2-5) )()()()()()( 212121 22 112121 02010 2020210 11 tnttstsegegetxt MNjMNjMj jjj NM 在理想情况下，假设阵列中各阵元是各向同性的且不存在通道不一致、互耦等因素的影响，则式(3.2-4)中的增益 可以省略 (即归一化 1)，在此假设下式(3.2-5)可以简化为li(3.2-6) )()()()()()( 212121 02010 2020210 11 tnttstseetxt MNjjj jjjMMMN 将式(3.2-6)写成矢量形式如下：(3.2-7

6、)()(tsAtx式中， 为阵列的 维快拍数据矢量， 为阵列的 维噪声数据矢量，)(tx1n1为空间信号的 维矢量 ，A为空间阵列的 维流型矩阵(导向矢量阵)，且)(tsNN(3.2-8)()(00201 Naa其中导向矢量(3.2-9)exp()()(02100Miii jjai,21式中 ，c为光速， 为波长。f20由上述的知识可知，一旦知道阵元间的延迟表达式 ，就很容易得出待定空间阵列的导向矢量或阵列流型。下面推导一下空间阵元间的延迟表达式。假设空间任意两个阵元，其中一个为参考阵元(位于原点)，另一个阵元的坐标为(x ,y, z)，两阵元的几何关系见图，图中“”表示阵元。图 2 空间任意

7、两阵元的几何关系由几何关系可以推导出两阵元的波程差为(3.2-10)sincosincos(1 zyx这里的波程差其实就是位于x轴上两阵元间的延迟、位于y轴上两阵元间的延迟和位于z轴上两阵元间的延迟之和。根据式(3.2-10)的结论，下面给出实际环境中常用的几种阵列及阵元间的相互延迟表达式。(1)平面阵 设阵元的位置为 ，以原点为参考点，另假设信),21)(,Mkyx号入射参数为 ，分别表示方位角与俯仰角，其中方位角表示与 x轴的,21)(,Nii夹角。(2)线阵设 阵元的位置为 ，以原点为参考点，另假设信号入射参),(k数为 ，表示方位角，其中方位角表示与y轴的夹角（即与线阵法线的夹角）),

8、21(i，则有(3.2-11)sin(1kkixc(3)均匀圆阵 设以均匀圆阵的圆心为参考点，则有(3.2-12)iiki Mcrcos)(2os其中方位角表示与x轴的夹角，r为圆半径。3.3 旋转不变子空间算法原理3.3.1 信号模型算法介绍前，首先对信号进行建模。为了推导分析的方便，将波达方向的数学模型做如下理想状态的假设：1) 阵列形式为线性均匀阵，阵元间距不大于信号波长的二分之一。2) 存生两个完全相同的子阵，且两个子阵的间距是己知的。3) 噪声序列为一零均值高斯过程，各阵元间噪声相互独立，噪声与信号也相互独立。4) 空间信号为零均值平稳随机过程，通常为窄带远场信号。5) 信号源数小于

9、子阵阵列元数，信号取样数大于子阵阵列元数，以确保子阵阵列流型的各列线性独立。6) 组成阵列的各传感器为各向同性阵元，且无互耦以及通道不一致的干扰。下图给出了均匀线阵的数学模型示意图:3.3.2 算法原理对于均匀线阵，相邻子阵间存在一个固定间距，这个固定间距反映出各相邻子阵间的一个固定关系，即子阵间的旋转不变性，而 ESPRIT 算法正是利用了这个子阵间的旋转不变性实现阵列的 DOA 估计。ESPRIT 算法最基本的假设是存在两个完全相同的子阵，且两个子阵的间距 是已知的。由于两个子阵的结构完全相同，且子阵的阵元数为 m，对于同一个信号而言，两个子阵的输出只有一个相位差 ， =1，2， N 。i

10、下面假设第一个子阵的接收数据为 ，第二个子阵的接收数据为 ，根1X2X图 3-1 均匀线阵的数学模型示意图1()Xtd 2()Xt3()t()MXt据前面所述的阵列模型可知(3.1)1 11()NXaSA(3.2)12 22()NjjeSN 式中，子阵 1 的阵列流型 =A，子阵 2 的阵列流型 = ，且式中1 2A1.Njjdiage（3.3）从上面的数学模型可知，需要求解的是信号的方向，而信号的方向信息包含在 A 和 中，由于 是一个对角阵，所以下面只考虑这个矩阵，即(2sin)/kk（3.4）由上可知。只要得到两个子阵间的旋转不变关系 ，就可以方便地得到关于信号到达角的信息。下面的任务就

11、是从式（3.1）和式（3.2）中得到两个子阵间的关系。先将两个子阵的模型进行合并，即12XASN（3.5）在理想条件下，可得上式的协方差矩阵(3.6) HHSNREXAR对上式进行特征分解可得21mHHiSsNeU（3.7）显然上式中得到的特征值有如下关 = ，U S 为大特1N12m征值对应的特征矢量张成的信号子空间， 为小特征值对应矢里张成的噪声子空间。对于实际的快拍数据，式(3.7)应修正如下：(3.8)HNHSUUR由前面的知识可知，上述的特征分解中大特征矢量张成的信号子空间与阵列流型张成的信号子空间是相等的。即(3.9)(AspanUsS此时，存在一个惟一的非奇异矩阵 T，使得S)(

12、（3.10）显然，上述的结构对两个子阵都成立，所以有TAUS21（3.11） 很显然 ，由子阵1的大特征矢量张成的子空间 、由子阵2的大特征矢量张成1SU的子空间 与阵列流型A张成的子空间三者相等，即2S21)(SSspanAspanUs（3.12） 另外，由两个子阵列在阵列流型上的关系可知12A（3.13）再利用式（3.11）可知两个子阵列的信号子空间的关系如下: 112SSUTU（3.14）式（3.13）反映了两个子阵列的阵列流型间的旋转不变性，而式（3.14）反映了两个子阵的阵列接收数据的信号子空间的旋转不变性。如果阵列流型A是满秩矩阵，则由式（3.14）可以得到1T（3.15）所以上式

13、中 的特征值组成的对角阵一定等于 ，而矩阵T的各列就是矩阵特征矢量。所以一旦得到上述的旋转不变关系矩阵 ，就可以直接利用式（3.4）得到信号的入射角度。3.4 标准的旋转不变子空间算法有上节的知识可知， ESPRIT算法的基本原理就是利用式（3.14）的旋转不变性，常规的旋转不变子空间算法就是利用上述的基本原理求解信号的入射角度信息。下面就分析解这个等式的两种最经典、应用最广泛方法:最小二乘（LS）法和总体最小二乘（TLS）法。3.4.1 最小二乘法由最小二乘的数学知识，我们知道式（3.14）的最小二乘解的方法等价于，约束条件 2minSU221SSSU（3.16）因此最小二乘法的基本思想就是

14、使校正项 尽可能小，而同时保证满足2S约束条件。为了得到LS解，将式（ 3.14）代入式（ 3.16）即得（3.17）212minin)(min SSSUUf 对上式进行展开可得=21)(SSf22121HHHSS SU（3.18）上式对 求导并令其等于0，可得02)( 11SHSUdf（3.19）上式的解显然有两种可能:(1) 当 满秩时，也就是子阵1的信号子空间的维数等于信号源数时，则SU上式的解是唯一的，可得上式的最小二乘解1212()()HLSSSSU（3.20） (2) 当 不满秩，即 时，也就是信号源间存在相干或相差时，1SUNrankS)(1则 存在很多解，但我们却无法区别对应于

15、方程的各个不同的解，可以称这些解是不可辨识的，解的不可辨识性是我们需要解相干的原因所在。下面给出 LS-ESPRIT 算法的求解步骤：1由两个子阵的接收数据 , ，分别得到两个子阵的数据协方差矩阵；1X22对矩阵对R, 进行特征分解，从而得到两个数据矩阵的信号子空间NR和 ；1SU23按式（3.20）得到矩阵 ，然后对其进行特征分解.得到 N 个特征值，就LS可得到对应的 N 个信号的到达角。当考虑嗓声影响时，上述基于最小二乘算法的估计都是有偏的，这就是为什么需要考虑总体最小二乘 ESPRIT 算法的原因。3.4.2 总体最小二乘法我们知道，普通最小二乘的基本思想是用一个范数平方为最小的扰动

16、2SU去于扰信号子空间 ，目的是校正 中存在的嗓声。显然这就存在一个问题：2SU2SU如果同时扰动 和 ，并使扰动范数的平方保持最小，是否可以同时校正1和 中存在的嗓声 ?答案是肯定的，这就是总休最小二乘（TLS）的思想。1SU2它考虑的是如下矩阵方程的解:(3.27)122()SSSUU显然上式可以改写成2121()()0SS z（3.28）所以 TLS 的解等价于2min()0Uz约 束 条 件 ：（3.29）定义如下一个矩阵 ，再结合上述分析过程。我们发现其实122|SS就是寻找一个 的酉矩阵 F，便得矩阵 F 与 正交，也就说明了由 F 张2Nm12SU成的空间与 或 列矢量张成的空间

17、正交。所以矩阵 F 可从 的特征1SU 12HSU分解中得到。因为12HHSE（3.30）式中的 是由特征值构成的对角矩阵，E 是与其相应的特征矢量构成的矩阵。即12E（3.31）令 是由对应特征值为 0 的特征矢量构成的矩阵.它属子噪声子空12NNE间，所以只要选择矩阵 F 使之等于 、，即可满足上面提到的要求。即有NE(3.32)1121212| 0SSSUUF可得 120ATF（3.33）如果令 ，则12FT（3.34）上式说明 的特征值即 是对角线元素。这说明通过构造一个矩阵 就可得到有关信号角度的信息.而这个矩阵的构造可通过式（3.30）得到，即12TLSE（3.35）下面直接给出

18、TLS-ESPRIT 算法的求解步骤：1由两个子阵的接收数据 ， , 由式(3.8)得到数据协方差矩阵 R；1X22通过矩阵对于 的广义特征分解，得到维数为 的信号子空,NR2MN间 ；SU3由 构造矩阵 ，并按式防（3.30）进行特征分解得到矩阵 E，然后12SU再按式(3.31) 将矩阵分为四个小的矩阵；4按式（3.35）得到矩阵 ，然后对其进行特征分解,得到 N 个特征值，TLS就可得到对应的 N 个信号的到达角。通过分析，我们可以得到标准 ESPRIT 算法的计算过程如下：(1)通过特征值或奇异值分解（EVD 或 SVD）分别估计两个存在旋转不变关系的子阵的信号子空；(2)用上述的 L

19、S、TLS 等方法求解式 （3.14）所示的不变等式；(3)计算 的特征值，其中 如式（3.3）所示。然后利用式(3.4)求1kT解人射信号的角度信息。就 ESPRIT 算法而言， TLS 算法与 LS 算法性能基本一致，只是在低信噪比情况下 TLS 算法性能略好。四、 仿真结果主要分析各个参数对估计误差的影响，误差函数定义如式(1)：4.1 信噪比 SNR对估计误差的影响分析首先对信噪比 SNR 离散化取值，然后求得不同信噪比下的误差，从而绘制出误差随信噪比改变的函数曲线如图 2 所示，图 2 中信噪比 SNR 从- 15 取到 15，间隔为 1，运行次数为 100 次，其余条件如题中所述。

20、由图 2 可知，随着信噪比的增大，估计误差会越来越小，即估计精度会越来越高。当待估计的信号方位角相差比较小时，估计的误差也会相应的增大。另外，若两信号为相干信号，则此方法将不能对其进行正确的估计。4.2 阵元数 L对估计误差的影响分析与 4.1 节类似，首先对阵元数 L 离散化取值，然后求得不同阵元数下的误差， 从而绘制出误差随阵元数改变的函数曲线如图 3 所示，图 3 中阵元数从 K+1 取到 K+25，间隔为 1，运行次数为 100 次，其余条件如题中所述。由于阵元数 L 需大 于信号个数 K 才能正确估计，故取值中含有信号个数 K。由图 3 可知，随着阵元数的增加，估计误差会越来越小，即

21、估计精度会越来越高，但当阵元数大到一定程度后，对估计精度的影响则会慢慢的减小。4.3 采样点数 N对估计误差的影响分析与 4.1 节类似，首先对采样点数 N 离散化取值，然后求得不同采样点数下的误差，从而绘制出误差随采样点数改变的函数曲线如图 4 所示，图 4 中采样点数从 10 取到 200，间隔为 5，运行次数为 100 次，其余条件如题中所述。由图 4 可知，随着采样点数的增加，估计误差会越来越小，即估计精度会越来越高。4.4 两信号之间的角度差（GAP）对估计误差的影响分析由于采用 ESPRI T 算法对 DOA 进行估计，若两信号的方位距离较近时，虽然能得出估计结果，但估计的精度会大

22、受影响。因此，为了分析两信号之间的不同间隔会对估计精度造成多大的影响，绘制不同 GAP 下的估计误差曲线如图 5 所示。处理方法与 4. 1 节类似，图 5 中 GAP（单位为度）从 0. 1 取到 5，间隔为 0. 1，独立运行次数为 100 次，其余条件如题中所述。由图 5 可知，GAP 越大估计越准确，但当 GAP 大到一定程度后则估计精度趋于稳定。估 计 误 差 （角 度）4.5 单信号 DOA不同分布对估计误差的影响分析信号波达方向（DOA）的取值区间为-90 度到 90 度，若只考虑只有一个信号的情况，则当信号的 DOA 不同时，估计误差也会不一样。因此，为了分析不同的 DOA 会

23、对估计精度造成多大的影响，绘制不同 DOA 下的估计误差曲线如图 6 所示。处理方法与 4. 1 节类似，图 6 中 GAP 从- 80 度取到 80 度，间隔为 5 度，独立运行次数为 100 次，其余条件如题中所述。由图 6 可知，DOA 越靠近 0 度估计越准确，越靠近正负 90 度估计误差越大。且仿真结果表明，当 DOA 在正负 90 附近时，估计误差太大，因此，为了不影响估计结果显示效果，故在图中未绘制正负 90 度附近的估计误差。2.6 减与不减噪声方差(Rn)对估计误差的影响分析由于有噪声的影响，因此在估计信号自相关矩阵 R 时，若将无信号时的自相关矩阵 Rn 减去，即相当与减去

24、估计出噪声方差，则估计的精度会有所提高。结合信噪比 SNR 对估计误差的影响，绘制减与不减噪声方差两种情况下估计误差随 SNR 的变化曲线如图 7 所示，图 7 中 SNR 从-15dB 到 5dB，间隔为 1dB，独立运行次数为 100 次。仿真结果表明，若减 Rn，主要是在低信噪比时对估计精度的改善较大，当信噪比较大时二者几乎一样。五、 程序清单%本文件名为 drawTLSesprit.m %分析基于总体最小二乘的 ESPRIT 算法（TLS-ESPRIT）的 DOA 估计的性能%clear;clc;close all; %清除变量，清屏，关闭所有绘图窗口% 调用格式:estimated,

25、error=TLSesprit(p,L,K,SNR,DOA);% 估计结果(弧度，矢量：p 行 1 列)：estimated% 估计误差(弧度，标量：均方误差)：error% 信号个数：p% 阵元数：L% 快拍数：K% 信噪比：SNR% 波达方向(弧度，矢量：p 行 1 列)：DOA% p=2; L=8; K=100; SNR=5; DOA=pi*(-10/180) pi*(20/180);%显示估计结果%M=100; %设定独立重复运行次数DOA=pi*(0/180) pi*(30/180); %波达方向(弧度，矢量：p 行 1 列)p=length(DOA); L=8; K=100; SN

26、R=5; %参数设置，estimated,error=TLSesprit(p,L,K,SNR,DOA); %函数调用polar(estimated,1 1,r*); %在极坐标中显示估计结果(必须先转化为弧度)h=title();set(h,string,TLS-ESPRIT: 估计值: ,num2str(estimated);h1=xlabel();set(h1,string,信号 DOA(度): ,num2str(DOA*180/pi);% %阵元数 L 对估计误差的影响分析% Ln=p+1:1:p+25; %阵元数 L 需大于信号个数 p 才能正确估计% for n=1:length(L

27、n)% L=Ln(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);% errorm(k)=error; %将每次的估计误差存入变量 errorm 中，便于求均值% end% errorn(n)=sum(errorm)/M; %求多次运行后的估计误差的均值% end% figure(2);plot(Ln,errorn*180/pi,r:*,LineWidth,2); %绘制曲线，并适当标注% xlabel(阵元数 L);ylabel(估计误差( );title(阵元数 L 对估计误差的影响);% % 结论：阵元数 L 越大估计越准

28、确，但当 L 大到一定程度后则估计精度趋于稳定% %快拍数 K 对估计误差的影响分析% Kn=10:10:200; %对快拍数离散化取值% for n=1:length(Kn)% K=Kn(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);% errorm(k)=error; %将每次的估计误差存入变量 errorm 中，便于求均值% end% errorn1(n)=sum(errorm)/M; %求多次运行后的估计误差的均值% end% figure(3);plot(Kn,errorn1*180/pi,r:*,LineWidth

29、,2); %绘制曲线并适当标注% xlabel(快拍数 K);ylabel(估计误差( );title(快拍数 K 对估计误差的影响);% % 结论：快拍数 K 越大估计越准确，但当 K 大到一定程度后则估计精度趋于稳定% %信噪比 SNR 对估计误差的影响分析% SNRn=-15:1:15; %对信噪比 SNR 离散化取值% for n=1:length(SNRn)% SNR=SNRn(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);% errorm(k)=error; %将每次的估计误差存入变量 errorm 中，便于求均值

30、% end% errorn(n)=sum(errorm)/M; %求多次运行后的估计误差的均值% end% figure(4);plot(SNRn,errorn*180/pi,r:*,LineWidth,2);%绘制曲线并适当标注(误差：角度)% xlabel(SNR);ylabel(估计误差( );title(SNR 对估计误差的影响);% %结论：信噪比 SNR 越大估计越准确，但当信噪比 SNR 大到一定程度后则估计精度趋于 稳定%两信号之间的角度差（GAP）的大小对估计误差的影响分析%GAPn=0.1:0.1:5; %对两信号之间的角度差（GAP）离散化取值for n=1:length

31、(GAPn)GAP=GAPn(n); %每次循环只取其中一个值DOA=pi*(0/180) pi*(GAP/180);for k=1:M %M 为独立重复运行次数estimated,error=TLSesprit(p,L,K,SNR,DOA);errorm(k)=error; %将每次的估计误差存入变量 errorm 中，便于求均值enderrorn(n)=sum(errorm)/M; %求多次运行后的估计误差的均值endfigure(5);plot(GAPn,errorn*180/pi,r:*,LineWidth,2);%绘制曲线并适当标注(误差：角度) xlabel(GAP( );ylab

32、el(估计误差( );% title(两信号之间的角度差（GAP）对估计误差的影响);%结论：GAP 越大估计越准确，但当 GAP 大到一定程度后则估计精度趋于稳定% %单个信号时，信号波达方向分布不同时对估计误差的影响分析% DOAn=pi*(-80/180):(5/180):pi*(80/180); %对信号波达方向离散化取值（80 度到 90 度时误差太大，因此未取）% for n=1:length(DOAn)% DOA=DOAn(n);p=1; %每次循环只取其中一个值,信号个数 p 设为为 1% for k=1:M %M 为独立重复运行次数% estimated,error=TLSe

33、sprit1(p,L,K,SNR,DOA); %调用 TLSesprit1（一个信号的情况）% errorm(k)=error; %将每次的估计误差存入变量 errorm 中，便于求均值% end% errorn(n)=sum(errorm)/M; %求多次运行后的估计误差的均值% end% figure(6);plot(DOAn*180/pi,errorn*180/pi,b:*,LineWidth,2);%绘制曲线并适当标注(误差：角度)% xlabel(DOA( );ylabel(估计误差( );% % title(DOA（单信号）不同分布对估计误差的影响);% % %结论：越靠近 0 度

34、估计越准确，越靠近正负 90 度估计误差越大% %估计相关矩阵 R 时，减与不减 Rn（无信号时的噪声自相关矩阵）对估计误差的 影响分析(结合信噪比 SNR 对估计误差的影响曲线)% SNRn=-15:1:5; %对信噪比 SNR 离散化取值% for n=1:length(SNRn)% SNR=SNRn(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);errorm(k)=error; %调用减 Rn 的函数% estimated,error=TLSespritRn(p,L,K,SNR,DOA);errormRn(k)=erro

35、r; %调用不减 Rn 的函数% end% errorn(n)=sum(errorm)/M; errornRn(n)=sum(errormRn)/M;% end% figure(7);h=plot(SNRn,errorn*180/pi,SNRn,errornRn*180/pi,r-.); %绘制减与不减 Rn 时的估计误差曲线% legend(R=R-Rn,R); set(h,LineWidth,2); %用图示在图中标明哪条为减或不减 Rn 的曲线% xlabel(SNR(dB);ylabel(估计误差( );% % title(减与不减 Rn 对估计误差的影响);% % %结论：若减 Rn

36、，主要是在低信噪比时对估计精度的改善较大，当信噪比较大 时二者几乎一样%本文件名为 TLSesprit.m %基于总体最小二乘的 ESPRIT 算法（TLS-ESPRIT）的 DOA 估计函数% function estimated,error=TLSesprit(p,L,K,SNR,DOA)% 调用格式:estimated,error=TLSesprit(p,L,K,SNR,DOA);% 估计结果(弧度，矢量：p 行 1 列)：estimated% 估计误差(弧度，标量：均方误差)：error% 信号个数：p% 阵元数：L% 快拍数：K% 信噪比：SNR% 波达方向(弧度，矢量：p 行 1

37、列)：DOA% p=2; L=8; K=100; SNR=5; DOA=pi*(-10/180) pi*(20/180);%参数设置%dbbc=1/2; %阵元间隔 d 与信号波长之比 d/ =1/2theta=2*pi*dbbc*sin(DOA); %信号方位参数 thetaOmigaT=pi/4; pi/6; %信号频率Dn=sqrt(1/(2*10(SNR/10); %噪声标准差%估计相关矩阵 R%A=exp(j*(0:L-1)*theta); %表示出阵列方向矩阵S=exp(j*OmigaT*(0:K-1); %构造信号源矢量X=A*S; %构造阵列输出矢量（无噪）Noise=Dn*(

38、sqrt(2)/2)*(randn(L,K)+j*randn(L,K);%加入复噪声Y=X+Noise; %构造阵列输出矢量R=zeros(L,L);Rn=zeros(L,L); %初始化为零，加快运行速度for i=1:KRn=Rn+Noise(:,i)*Noise(:,i);R=R+Y(:,i)*Y(:,i);endR=R/K;Rn=Rn/K; %求得相关矩阵 R(有信号)和 Rn（无信号）R1=R-Rn; %减小噪声对估计精度的影响V,D=eig(R1); %相关矩阵特征分解%（D 中特征值已经按从小到大的顺序排列，即 V 中前 L-p 个为噪声对应的特征向量）%构造矩阵 S%S=V(:

39、,L-p+1:L); %L 行 p 列（S 中的列为 R 中 p 个大特征值对应的特征向量）S1=S(1:L-1,:); %将 S 的前 L-1 行构造 S1S2=S(2:L,:); %将 S 的后 L-1 行构造 S2S12=S1 S2; %利用 S1 和 S2 构造 S12（L-1 行 2p 列）SS=S12*S12; %2K 行 2p 列U,D1=eig(SS); %特征分解%求解估计结果%U11=U(1:p,1:p); %p 行 p 列，U12=U(p+1:2*p,1:p); %p 行 p 列TLS=-U11*inv(U12); %p 行 p 列，U11 和 U12 构成 U 的噪声子空间d=eig(TLS); %特征分解,求 TLS 的特征值estimated=(sort(asin(angle(d)/pi); %输出估计(已从小到大排序)结果（弧度）error=sqrt(sum(estimated-sort(DOA).2)/p); %求出估计误差(弧度):均方误差