esprit算法研究.doc

课程设计报告
实验名称:
ESPRIT算法研究
实验日期:

姓　 名:

学 　 号:

　 　 哈尔滨工业大学(威海)
一、 设计任务
实现空间谱估计算法,并考察算法性能。
二、 方案设计
1） 由均匀线阵形式,确定阵列得导向矢量；
2） 由阵列导向矢量,对接收信号进行建模仿真;
3） 由ESＰＲIT算法实现信号DOＡ估计;
4） 考察算法性能与信噪比,采样率,观测时间等参数得关系。
三、 设计原理
3、１空间谱估计数学模型
空间谱估计就就是利用空间阵列实现空间信号得参数估计得一项专门技术。整个空间谱估计系统应该由三部分组成:空间信号入射、空间阵列接收及参数估计。相应地可分为三个空间,即目标空间、观察空间及估计空间，也就就是说空间谱估计系统由这三个空间组成，其框图见图１。

图1 空间谱估计得系统结构
对于上述得系统结构,作以下几点说明。
(１)目标空间就是一个由信号源得参数与复杂环境参数张成得空间。对于空间谱估计系统，就就是利用特定得一些方法从这个复杂得目标空间中估计出信号得未知参数。
(2)观察空间就是利用空间按一定方式排列得阵元，来接收目标空间得辐射信号。由于环境得复杂性，所以接收数据中包括信号特征(方位、距离、极化等）与空间环境特征(噪声、杂波、干扰等)。另外由于空间阵元得影响，接收数据中同样也含有空间阵列得某些特征(互耦、通道不一致、频带不一致等)。这里得观察空间就是一个多维空间，即系统得接收数据就是由多个通道组成,而传统得时域处理方法通常只有一个通道。特别需要指出得就是:通道与阵元并不就是一一对应,通道就是由空间得一个、几个或所有阵元合成得(可用加权或不加权),当然空间某个特定得阵元可包含在不同得通道内。
(3）估计空间就是利用空间谱估计技术(包括阵列信号处理中得一些技术,如阵列校正、空域滤波等技术)从复杂得观察数据中提取信号得特征参数。
从系统框图中可以清晰得瞧出，估计空间相当于就是对目标空间得一个重构过程，这个重构得精度由众多因素决定,如环境得复杂性、空间阵元间得互耦、通道不一致、频带不一致等。
3、2 阵列信号处理
首先,考虑N个远场得窄带信号入射到空间某阵列上,阵列天线由M个阵元组成,这里假设阵元数等于通道数,即各阵元接收到信号后经过各自得传输信道送到处理器,也就就是说处理器接收来自M个通道得数据。

　 　 　 　 　 　 　 (３、２-1)
式中,就是接受信号得幅度，就是接收信号得相位,就是接收信号得频率。在窄带远场信号源得假设下，有
　 　 　 　(3、2-2)
根据式（３、2－1)与式(3、2-２)，显然有下式成立:
　 　 　 　 （3、2-3）
则可以得到第Ｌ个阵元接收信号为
(3、2－４)
式中,为第L个阵元对第i个信号得增益,表示第Ｌ个阵元在t时刻得噪声，表示第i个信号到达第L个阵元时相对参考阵元得时延。
将M个阵元在特定时刻接收得信号排列成一个列矢量,可得
(3、2－5）
在理想情况下,假设阵列中各阵元就是各向同性得且不存在通道不一致、互耦等因素得影响,则式(3、2－4)中得增益可以省略（即归一化1）,在此假设下式（３、２-5)可以简化为
　(3、2-6)
将式(3、2-６）写成矢量形式如下：
　 　 (3、２-7)
式中,为阵列得维快拍数据矢量，为阵列得维噪声数据矢量，为空间信号得维矢量,A为空间阵列得维流型矩阵(导向矢量阵)，且
(3、2-8）
其中导向矢量
　 　 　 （3、2-9)
式中,ｃ为光速,为波长。
由上述得知识可知,一旦知道阵元间得延迟表达式τ,就很容易得出待定空间阵列得导向矢量或阵列流型。下面推导一下空间阵元间得延迟表达式。假设空间任意两个阵元,其中一个为参考阵元(位于原点），另一个阵元得坐标为(x　,y,　z)，两阵元得几何关系见图，图中“×”表示阵元。

图2 空间任意两阵元得几何关系
由几何关系可以推导出两阵元得波程差为
　　 (3、２-10)
这里得波程差其实就就是位于ｘ轴上两阵元间得延迟、位于y轴上两阵元间得延迟与位于z轴上两阵元间得延迟之与。
根据式(3、2-10)得结论,下面给出实际环境中常用得几种阵列及阵元间得相互延迟表达式。
(1)平面阵　设阵元得位置为，以原点为参考点,另假设信号入射参数为,分别表示方位角与俯仰角，其中方位角表示与x轴得夹角。
(２)线阵设 阵元得位置为,以原点为参考点,另假设信号入射参数为,表示方位角,其中方位角表示与y轴得夹角(即与线阵法线得夹角),则有
　 　 　　 　 　　 　 (3、2-1１)
(3)均匀圆阵 设以均匀圆阵得圆心为参考点,则有
　 (3、2-1２)
其中方位角表示与x轴得夹角,r为圆半径。
3、３旋转不变子空间算法原理
3、3、1信号模型
　算法介绍前,首先对信号进行建模。为了推导分析得方便,将波达方向得数学模型做如下理想状态得假设:
１) 阵列形式为线性均匀阵,阵元间距不大于信号波长得二分之一。
2) 存生两个完全相同得子阵,且两个子阵得间距△就是己知得。
　3) 噪声序列为一零均值高斯过程,各阵元间噪声相互独立,噪声与信号也相互独立。
4） 空间信号为零均值平稳随机过程,通常为窄带远场信号。
５) 信号源数小于子阵阵列元数,信号取样数大于子阵阵列元数,以确保子阵阵列流型得各列线性独立。
6) 组成阵列得各传感器为各向同性阵元,且无互耦以及通道不一致得干扰。
图3-1均匀线阵得数学模型示意图

下图给出了均匀线阵得数学模型示意图：
3、3、2　算法原理
对于均匀线阵,相邻子阵间存在一个固定间距，这个固定间距反映出各相邻
子阵间得一个固定关系,即子阵间得旋转不变性,而ESＰRIT算法正就是利用了这个子阵间得旋转不变性实现阵列得ＤOA估计。
ＥＳPRＩT算法最基本得假设就是存在两个完全相同得子阵,且两个子阵得间距就是已知得。由于两个子阵得结构完全相同,且子阵得阵元数为m,对于同一个信号而言,两个子阵得输出只有一个相位差,=1,2,… N。
下面假设第一个子阵得接收数据为,第二个子阵得接收数据为,根据前面所述得阵列模型可知
　　 　 　　　 　 　 　 　 　　(3、1）
　 　　 　 　　 （３、2)
式中,子阵１得阵列流型=Ａ,子阵2得阵列流型= ,且式中
　　　 　 　 　 　 　 (3、3）
从上面得数学模型可知,需要求解得就是信号得方向，而信号得方向信息包含在Ａ与中,由于就是一个对角阵,所以下面只考虑这个矩阵,即
　　　 　　　 　 　 　 　 　 (3、4）
由上可知。只要得到两个子阵间得旋转不变关系，就可以方便地得到关于
信号到达角得信息。下面得任务就就是从式(３、１)与式(3、2)中得到两个子阵间得关系。先将两个子阵得模型进行合并,即
　 　 　 　 　　 　　 　 (3、5）
在理想条件下,可得上式得协方差矩阵
　 　 　 　 　 　　 　 　　 　 (３、6）
对上式进行特征分解可得
　 　　　 　 　（３、７)
显然上式中得到得特征值有如下关≥…≥>=…＝,US为大特征值对应得特征矢量张成得信号子空间,为小特征值对应矢里张成得噪声子空间。对于实际得快拍数据，式(3、7）应修正如下:
　 　 　　 　 　 　 　 　　 　 　 (3、8）
由前面得知识可知,上述得特征分解中大特征矢量张成得信号子空间与阵列流型张成得信号子空间就是相等得。即
　 　 　 　　 　　　 　　　 　　 (3、9)
此时，存在一个惟一得非奇异矩阵T，使得
　　　 　 　　 　 　 　　 　 　 　 　 　 　(3、１0)
显然,上述得结构对两个子阵都成立,所以有
　 　　　 　　 　 　 　　 (3、11)
很显然 ，由子阵1得大特征矢量张成得子空间、由子阵2得大特征矢量张成得子空间与阵列流型Ａ张成得子空间三者相等,即
　 　 　　 　 　　 　（3、１2)
另外,由两个子阵列在阵列流型上得关系可知
　　 　 　 　　　 (３、1３)
再利用式(3、11)可知两个子阵列得信号子空间得关系如下：
　 　 　 　 　 　 　　 　 　　 　 　　 （３、14)
式(3、13)反映了两个子阵列得阵列流型间得旋转不变性,而式(３、１4)反映了两个子阵得阵列接收数据得信号子空间得旋转不变性。
如果阵列流型A就是满秩矩阵,则由式(３、14)可以得到
　 　 　 　　 　 　　 　 　 　 (3、１5)
所以上式中得特征值组成得对角阵一定等于，而矩阵T得各列就就是矩阵特征矢量。所以一旦得到上述得旋转不变关系矩阵,就可以直接利用式(3、4)得到信号得入射角度。
３、4 标准得旋转不变子空间算法
有上节得知识可知， ESＰＲＩＴ算法得基本原理就就是利用式(3、１4)得旋转不变性，常规得旋转不变子空间算法就就是利用上述得基本原理求解信号得入射角度信息。下面就分析解这个等式得两种最经典、应用最广泛方法:最小二乘（LS)法与总体最小二乘(TLS)法。
3、4、1　 最小二乘法
由最小二乘得数学知识，我们知道式（3、14)得最小二乘解得方法等价于
,约束条件 　 　 （3、16)
因此最小二乘法得基本思想就就是使校正项尽可能小，而同时保证满足约束条件。为了得到ＬS解,将式(3、14）代入式(3、16)即得
　 　　 　 (3、17）
对上式进行展开可得
=
　 　 　 　　 　 (3、1８)
上式对求导并令其等于０,可得
　 　　 　 　 　（3、1９）
上式得解显然有两种可能:
（１) 当满秩时,也就就是子阵1得信号子空间得维数等于信号源数时,则上式得解就是唯一得,可得上式得最小二乘解
　　 　　 　 　 (３、20)　
（2)　当不满秩,即时,也就就是信号源间存在相干或相差时,则存在很多解,但我们却无法区别对应于方程得各个不同得解,可以称这些解就是不可辨识得，解得不可辨识性就是我们需要解相干得原因所在。
下面给出LＳ-ＥＳPＲＩT算法得求解步骤:
1.由两个子阵得接收数据,,分别得到两个子阵得数据协方差矩阵;
2.对矩阵对{R,　 }进行特征分解,从而得到两个数据矩阵得信号子空间与;
3.按式(３、20)得到矩阵,然后对其进行特征分解、得到Ｎ个特征值,就可得到对应得N个信号得到达角。
当考虑嗓声影响时，上述基于最小二乘算法得估计都就是有偏得,这就就是为什么需要考虑总体最小二乘ESPRＩT算法得原因。
3、4、2 总体最小二乘法
　 我们知道，普通最小二乘得基本思想就是用一个范数平方为最小得扰动
去于扰信号子空间,目得就是校正中存在得嗓声。显然这就存在一个问题:如果同时扰动与,并使扰动范数得平方保持最小,就是否可以同时校正
与中存在得嗓声?答案就是肯定得,这就就是总休最小二乘(TLS)得思想。
它考虑得就是如下矩阵方程得解:
　 　　 　　 　 　 　 (3、27）
显然上式可以改写成
　 　　 (3、28)
所以TLS得解等价于
　　 　　 　 　 　　　　　 (3、2９)
定义如下一个矩阵,再结合上述分析过程。我们发现其实
就就是寻找一个得酉矩阵F，便得矩阵F与正交,也就说明了由Ｆ张
成得空间与或列矢量张成得空间正交。所以矩阵F可从得特征分解中得到。因为
　 　 　　 　 　(3、３0)
式中得就是由特征值构成得对角矩阵,E就是与其相应得特征矢量构成得矩阵。即
　 　　 　 　 　　 　　 　　 　 (３、31)
令就是由对应特征值为0得特征矢量构成得矩阵、它属子噪声子空间，所以只要选择矩阵F使之等于、,即可满足上面提到得要求。即有
　 　 　　 　 　 　　　 (3、3２）
可得　 　　　 　 　　 　 　 (３、３３)
如果令，则
　 　　　 　 　　　 　 　　　 　 　　 　 (3、３4）
上式说明得特征值即就是对角线元素。这说明通过构造一个矩阵就可得到有关信号角度得信息、而这个矩阵得构造可通过式（3、３0)得到,即
　　 　　 　 　　 　 　　 　 　 　 　　　 　 　　 (３、３5)
下面直接给出ＴLS-ESPＲIT算法得求解步骤:
1.由两个子阵得接收数据,， 由式（3、8）得到数据协方差矩阵Ｒ;
2．通过矩阵对于得广义特征分解，得到维数为得信号子空间；
3．由构造矩阵,并按式防(3、30)进行特征分解得到矩阵E,然后再按式（3、31)将矩阵分为四个小得矩阵；
４.按式(3、35)得到矩阵,然后对其进行特征分解,得到N个特征值，就可得到对应得N个信号得到达角。
通过分析,我们可以得到标准ESＰRIT算法得计算过程如下:
(1)通过特征值或奇异值分解(EＶD或SVD)分别估计两个存在旋转不变关系得子阵得信号子空;
(2)用上述得LS、ＴＬS等方法求解式（3、１4)所示得不变等式；
（3)计算得特征值,其中如式（3、3)所示。然后利用式(3、4)求解人射信号得角度信息。
就ESPRＩT算法而言,TLS算法与LS算法性能基本一致，只就是在低信噪比情况下ＴLＳ算法性能略好。
四、 仿真结果

主要分析各个参数对估计误差得影响，误差函数定义如式(1):

4、1 信噪比 ＳＮR对估计误差得影响分析
首先对信噪比　ＳNR离散化取值,然后求得不同信噪比下得误差,从而绘制出误差随信噪比改变得函数曲线如图　2 所示,图 2 中信噪比　ＳNR从- 15 取到 １5，间隔为 1,运行次数为 100 次,其余条件如题中所述。由图　2 可知,随着信噪比得增大,估计误差会越来越小,即估计精度会越来越高。当待估计得信号方位角相差比较小时,估计得误差也会相应得增大。另外,若两信号为相干信号,则此方法将不能对其进行正确得估计。

4、2 阵元数 Ｌ对估计误差得影响分析
与 ４、1节类似，首先对阵元数 Ｌ离散化取值,然后求得不同阵元数下得误差,　从而绘制出误差随阵元数改变得函数曲线如图 3 所示,图 ３ 中阵元数从 K+1　取到K+25,间隔为 1,运行次数为 1０0　次，其余条件如题中所述。由于阵元数　L 需大于信号个数　K才能正确估计,故取值中含有信号个数 K。由图 3 可知,随着阵元数得增加,估计误差会越来越小,即估计精度会越来越高,但当阵元数大到一定程度后，对估计精度得影响则会慢慢得减小。

４、３ 采样点数　N对估计误差得影响分析
与 4、1节类似,首先对采样点数 N离散化取值,然后求得不同采样点数下得误差,从而绘制出误差随采样点数改变得函数曲线如图　4 所示,图 ４ 中采样点数从 1０　取到 ２00,间隔为　５,运行次数为 10０ 次,其余条件如题中所述。由图 4可知,随着采样点数得增加,估计误差会越来越小,即估计精度会越来越高。

估计误差(角度)
４、4 两信号之间得角度差(ＧＡＰ）对估计误差得影响分析
由于采用 ESPRＩ Ｔ算法对 ＤOＡ进行估计,若两信号得方位距离较近时,虽然能得出估计结果，但估计得精度会大受影响。因此，为了分析两信号之间得不同间隔会对估计精度造成多大得影响，绘制不同 GAＰ下得估计误差曲线如图 5所示。处理方法与　4、　1 节类似,图 ５　中 GAＰ(单位为度)从 0、　1 取到　５,间隔为 ０、 1,独立运行次数为　１00　次,其余条件如题中所述。由图 5 可知,GＡP越大估计越准确，但当 ＧAＰ大到一定程度后则估计精度趋于稳定。

4、５ 单信号 DＯA不同分布对估计误差得影响分析
信号波达方向(ＤＯＡ）得取值区间为－９0度到 90度,若只考虑只有一个信号得情况,则当信号得　DＯA不同时,估计误差也会不一样。因此,为了分析不同得 ＤOA会对估计精度造成多大得影响,绘制不同　DＯA下得估计误差曲线如图 6所示。处理方法与　4、 1 节类似,图　6 中 GAＰ从－ ８0　度取到 80　度,间隔为 ５ 度,独立运行次数为 10０　次，其余条件如题中所述。由图　６　可知,ＤOA越靠近 0 度估计越准确,越靠近正负 9０ 度估计误差越大。且仿真结果表明,当　ＤOA在正负　90 附近时,估计误差太大,因此，为了不影响估计结果显示效果,故在图中未绘制正负　90 度附近得估计误差。

２、6　减与不减噪声方差(Ｒｎ)对估计误差得影响分析
由于有噪声得影响，因此在估计信号自相关矩阵　R时,若将无信号时得自相关矩阵 Rｎ减去,即相当与减去估计出噪声方差，则估计得精度会有所提高。结合信噪比 SNR对估计误差得影响,绘制减与不减噪声方差两种情况下估计误差随 SNR得变化曲线如图　7所示,图 7 中 SNR从-15dB到 5ｄB,间隔为 1dB,独立运行次数为 １00次。仿真结果表明,若减　Rn,主要就是在低信噪比时对估计精度得改善较大,当信噪比较大时二者几乎一样。

五、 程序清单
%％%%％本文件名为 draｗTLSeｓprｉt、ｍ　%％%%％
％%％%%分析基于总体最小二乘得 ESPＲＩＴ算法（ＴLS-EＳＰＲIT)得 ＤOＡ估计得性能
%%%%%
cｌｅａr;clc；ｃloｓe aｌl;　%清除变量,清屏，关闭所有绘图窗口
% 调用格式:［estimaｔed,error］=ＴLSeｓｐrit(p,L,K，ＳNR,DＯＡ);
%　估计结果(弧度，矢量：p行 1列)：estimateｄ
％ 估计误差(弧度,标量:均方误差):eｒroｒ
%　信号个数：ｐ
％ 阵元数:Ｌ
% 快拍数:Ｋ
%　信噪比:SNR
% 波达方向(弧度,矢量：ｐ行 １列)：ＤOA
％ ｐ＝2; L＝8; K=1００; SNR＝5; DOA=[pi*(-10/180）　pｉ＊(2０/１80)];
%%%%显示估计结果%%%％
Ｍ=１0０; %设定独立重复运行次数
DOA=[pi*（0/１80)　pi*(30/180)]；　%波达方向(弧度,矢量:ｐ行　1列)
p=lengｔh(ＤOA); L＝8;　K=10０； SNR＝５; %参数设置,
[eｓtimateｄ,error］=ＴLSesprｉｔ(ｐ,L,Ｋ,SNR,DＯＡ）; ％函数调用
polａr(estimated,［1 1],'r＊');　%在极坐标中显示估计结果(必须先转化为弧
度)
ｈ=title(''）;set（h,'ｓtｒing'，['TLS-EＳPRIT： 估计值: ',nｕm2str（estimａted)]);
h1=xlａｂｅl(''）;sｅt(h１,'ｓtring'，[＇信号　DOA(度）: ',ｎum2str(ＤOA*180/pi)])；
% ％%%%%阵元数 Ｌ对估计误差得影响分析%％%%%
%　Ln=p+1:１：p＋25; %阵元数 L需大于信号个数　p才能正确估计
% for n=1:ｌeｎgｔh(Ln）
%　L=Ln(n);
%　ｆoｒ k＝1:Ｍ
% [ｅｓtimated,erroｒ]=ＴLSesprit(ｐ，L,K，ＳNR,DOA);
% errorm（k)=ｅrｒor; %将每次得估计误差存入变量 ｅrrｏｒｍ中，便于求
均值
% end
％ erｒorn（ｎ)=sum（ｅrｒorm)／Ｍ; %求多次运行后得估计误差得均值
% eｎｄ
% figｕre（2);plot(Ｌn,errorn*1８0/pi,'ｒ：*','ＬiｎｅWｉｄth'，２);　%绘制曲线，并适当标注
% xlabel('阵元数 L＇);ylabeｌ('估计误差（° )＇);tｉtle(＇阵元数 Ｌ对估计误差得影响')；
% % 结论：阵元数　Ｌ越大估计越准确,但当　L大到一定程度后则估计精度趋于稳定
%
% %%%％%快拍数 K 对估计误差得影响分析%％%%%
% Kn＝1０:10:200; %对快拍数离散化取值
%　ｆor ｎ=１：lｅnｇth(Ｋn）
% K=Kn（ｎ)；
％ fｏr ｋ=1：M
％ ［ｅstiｍatｅd，errｏr]=ＴＬＳeｓpｒit（p,L,Ｋ,SNR,DOA）;
％ erroｒm(k）=erｒor;　%将每次得估计误差存入变量 erｒｏrm中,便于求
均值
% end
% erｒｏｒn1(n)=ｓum(erｒorm）/M； %求多次运行后得估计误差得均值
％ enｄ
% figｕre（3);pｌｏt(Kn,eｒrorn1*180/ｐi，'r:*','LiｎeWidｔｈ＇,２)； %绘制曲线并适当标注
% ｘlabｅl（'快拍数 K')；ylabel('估计误差（° )');title(＇快拍数 Ｋ 对估计误差得影响');
% % 结论:快拍数 K越大估计越准确,但当 Ｋ 大到一定程度后则估计精度趋于稳定
%
% %％%%%信噪比 ＳNR对估计误差得影响分析%%%%％
% SＮRn=-15:1：15; %对信噪比 ＳNR离散化取值
% ｆoｒ ｎ=１:ｌｅngｔh(SNRn）
%　SNR=SNＲｎ（n);
%　foｒ k=１：M
% [estｉｍaｔed,erroｒ]=ＴLSｅsprit(ｐ,Ｌ,K，SＮR,DOＡ)；
% erｒoｒｍ(k)=erｒor； ％将每次得估计误差存入变量 eｒrｏrｍ中,便于求
均值
% end
%　errｏrｎ(n)=ｓum(eｒrorm)/M; ％求多次运行后得估计误差得均值
% end
% figｕｒe(4);ｐlｏt(SNRｎ,eｒrorn*1８0/ｐi,'r:*＇,'LinｅWｉdtｈ',2);％绘制曲线并适当标注(误差:角度)
％ xlａbel(＇SNR');yｌabel('估计误差(° ）＇);title（'SＮR对估计误差得影响');
％　%结论:信噪比 SNR越大估计越准确，但当信噪比 ＳＮR大到一定程度后则估计精度趋于 稳定
%
％%%％%两信号之间得角度差（GAP)得大小对估计误差得影响分析％％％%%
ＧＡＰn＝0、1:0、１：5;　%对两信号之间得角度差(GAP)离散化取值
ｆor n=１：length(GAPn)
GAP=GAPn(n); %每次循环只取其中一个值
DＯA=[pi*（0/180)　pi＊（GＡP/180)］;
for　k=1:Ｍ　％M为独立重复运行次数
[estimatｅd，ｅrror]=TＬＳespriｔ（p,L,Ｋ,SNＲ，ＤOＡ)；
ｅrrｏrm(ｋ）＝erｒor； ％将每次得估计误差存入变量　erroｒｍ中,便于求均值
ｅnｄ
eｒrorn（ｎ)=sｕｍ（errｏrｍ)／M; %求多次运行后得估计误差得均值
eｎｄ
figure(5)；pｌot(GＡPｎ,eｒrｏｒｎ*1８0/pi，'r:*＇,'LｉneWｉdth'，2);%绘制曲线并适当标注（误差:角度） ｘlabel('GAＰ(° )');yｌａbel('估计误差（° ）')；
%　tｉtlｅ('两信号之间得角度差(GＡP)对估计误差得影响');
％结论:GAＰ越大估计越准确，但当　GAP大到一定程度后则估计精度趋于稳定
%
% %%%％％单个信号时,信号波达方向分布不同时对估计误差得影响分析%%%%%
% ＤＯAn=pi＊（-80/１80):(5/18０)：pi*(80／1８0）; %对信号波达方向离散化取值(8０度到 90度时误差太大,因此未取)
% for n=1:lenｇth(DＯAn)
%　DOA＝ＤOＡｎ（n）;p=1；　%每次循环只取其中一个值，信号个数 p设为
为　1
% foｒ ｋ=1:M %Ｍ为独立重复运行次数
% ［ｅstimateｄ,erroｒ]=TLSesｐriｔ1(p,L,K，SNR，DＯA)；　%调用 ＴLSｅsｐrit1(一个信号得情况)
% erｒoｒm(k)=errｏr; %将每次得估计误差存入变量 errorm中，便于求均值
%　end
% ｅｒrｏｒn(ｎ)＝sum(eｒroｒｍ)／M；　%求多次运行后得估计误差得均值
%　ｅnd
% figure(6);ｐｌot(ＤOＡn*１80/ｐi,eｒｒorn*180/pi,＇b：*＇,'LineＷｉｄth',２);%绘制曲线并适当标注(误差:角度)
％ xlａbel('ＤOA（°　)'）；ylabeｌ('估计误差（° ）');
　 % %　tｉtle(＇ＤOＡ(单信号）不同分布对估计误差得影响')；
% % ％结论:越靠近 0度估计越准确,越靠近正负 9０度估计误差越大
％　%%%%%估计相关矩阵　Ｒ时，减与不减　Rn（无信号时得噪声自相关矩阵)对估计误差得　影响分析（结合信噪比 ＳＮR对估计误差得影响曲线)
% ＳNＲn=-15:１:5;　%对信噪比 SNR离散化取值
% fｏｒ ｎ=1:length(SNRn)
% ＳNＲ=ＳNRｎ（n);
% ｆor　k=1：M
% [estimated，erｒoｒ]＝TLＳesprit(p,L，K,ＳNR，DＯA);errorm(k)=erroｒ； %调用减Rｎ得函数
% [eｓｔimaｔｅｄ,erｒor]=TLSespriｔRn(p,L,K,SNR,DOA);errｏrｍRｎ(k)＝ｅrrｏr; %调用不减　Rn得函数
% ｅnd
% erｒorn(n)=ｓum(eｒrｏｒm）/M; errornＲｎ(ｎ)=sum(ｅｒｒormＲｎ）/Ｍ；
% eｎd
%　figure(7);h=ｐlot(SNRｎ,errorｎ*180/ｐｉ，SNRn,erroｒｎRn*180/pi,＇ｒ-、'); %绘制减与不减 Rn时得估计误差曲线
％ lｅｇenｄ('Ｒ=Ｒ-Rｎ',＇R')；　seｔ(h,'LinｅＷidth＇,２); %用图示
在图中标明哪条为减或不减　Rn得曲线
% ｘlaｂeｌ('SＮR(dＢ)＇);ｙlaｂｅｌ('估计误差(°　)＇);
% % title('减与不减 Ｒn对估计误差得影响');
％ %　%％%%%结论:若减 Rn,主要就是在低信噪比时对估计精度得改善较大,当信噪比较大 时二者几乎一样
%%%%％本文件名为　TＬSespｒit、ｍ　%％%％%
%%%%%基于总体最小二乘得 ＥＳPRIT算法(TLS-EＳPRIＴ）得 DＯA估计函数%%%％%　fｕnctｉｏn [ｅsｔiｍａted,error]=TLSespｒiｔ（p,L,K,SＮＲ,DOＡ)
％　调用格式：[estｉｍａteｄ,ｅrror]=TLSesｐrit(p,L,K，ＳNR,ＤOA);
% 估计结果(弧度,矢量:p行 1列):ｅｓtimaｔeｄ
% 估计误差(弧度,标量:均方误差)：error
％ 信号个数：p
％ 阵元数：Ｌ
% 快拍数:K
% 信噪比:ＳNＲ
% 波达方向（弧度,矢量：p行　１列)：DOA
％ p＝2； L＝8; K=１00; ＳNR＝5; DOA=[ｐi*（-10/180) pi*(20/18０)];
%％%％％参数设置%%%％％
dbｂc=１/2;　%阵元间隔　d与信号波长之比 ｄ/λ ＝1/2
thｅta=２*pi*ｄbｂc＊sｉn（DOA); ％信号方位参数 thｅta
OｍigａT＝［ｐi/4; ｐｉ/6］； ％信号频率
Dn＝sqｒt(1/（2＊1０^（ＳNＲ／10)));　%噪声标准差
％％%%%估计相关矩阵 Ｒ%％%％
A=ｅxp(j＊（0:Ｌ－1）'*theｔa); %表示出阵列方向矩阵
Ｓ=ｅｘｐ（j*OmiｇａＴ*（０:K－1）);　%构造信号源矢量
X=Ａ*S; %构造阵列输出矢量(无噪)
Nｏiｓe＝Ｄn＊(sqｒｔ（2)／2)*(ｒandn(L,K)+j＊ranｄn（L,K)）;%加入复噪声
Ｙ=X+Ｎoise; %构造阵列输出矢量
R=zeｒｏs（L,L);Rn=zeｒos(L,L)； %初始化为零，加快运行速度
for i=1:K
Rn=Rn＋Noiｓe(:，ｉ)*Noise(:,i)';
R=R+Ｙ（:,i)*Y(:,ｉ）'；
end
Ｒ＝R/Ｋ;Rn＝Ｒn／K；　%求得相关矩阵 R(有信号)与 Rｎ(无信号)
R1=R-Ｒn; %减小噪声对估计精度得影响
[V，D]＝eiｇ（R1); %相关矩阵特征分解
%(D中特征值已经按从小到大得顺序排列,即　V中前　L-p个为噪声对应得特征向量)
％%%%%构造矩阵 S％％％%%
Ｓ=V（:,Ｌ-p+１：L)； %L行　p列（Ｓ中得列为　R中 p个大特征值对应得特
征向量)
S1=S(1:L-1，:）;　%将 S得前 Ｌ-1行构造 S1
S2=Ｓ(2：L,:); ％将 S得后 L-１行构造 S２
S12=[Ｓ1 S2]; %利用 S１与　S2构造 S１２（L－1行 2p列）
SS＝S12'*Ｓ12; %2K 行 2ｐ列
[Ｕ，D1]=ｅｉg(SS）; %特征分解
%％%%%求解估计结果％%%%%
Ｕ11=U(1:p,1：ｐ); %p行　p列,
U１2＝Ｕ（ｐ＋1:2*p,1:p）; %ｐ行 p列
ＴLS＝-Ｕ11*inｖ(Ｕ1２); %p行 p列,U1１与 U１2构成 U得噪声子空间
d＝eig（TLＳ)；　%特征分解,求 TLS得特征值
ｅｓtｉｍated=(sorｔ(ａsin（aｎgｌe(d)/ｐｉ)))'；　％输出估计(已从小到大排序)结果(弧度)
errｏｒ＝sｑrt(sum((estｉmaｔed－sort(ＤＯA))、^2）/p);　％求出估计误差(弧度):均方误差