毕业论文_非参数谱估计算法与实现

非参数谱估计的算法与实现内容摘要功率谱估计是随机信号分析中的一个重要内容，主要研究信号在频域中的各种特征，目的是根据有限观测数据在频域内提取被噪声污染的有用信号。本设计针对非参数谱估计的算法与实现展开研究。首先讨论了一种最常见的谱估计法周期图法。由于周期图法是有偏估计，其方差过大，分辨率较低，因此，需要研究基于周期图的各种改进方法。主要的改进途径包括改善窗口形状、对数据进行平均和平滑等。为此，讨论了基于周期图的其他几种非参数谱估计算法（即BT法、平均周期图法，包括BARTLETT法和WELCH法）及其通过MATLAB的仿真实现。此外，信号来波方向估计是功率谱估计算法在空间领域的推广应用，具体就是以多元天线阵结合现代数字信号处理为基础的新型测向技术。为此，本设计还讨论比较了采用CAPON法、MUSIC算法实现波达方向估计的实现。关键词周期图法平均周期图法BT谱估计CAPON法MUSIC算法波达方向ALGORITHMSANDIMPLEMENTATIONSOFNONPARAMETERSPECTRALESTIMATIONSABSTRACTPOWERSPECTRUMESTIMATEISANIMPORTANTPARTOFTHERANDOMSIGNALANALYSIS,WHICHPRIMARILYSTUDIESTHERANDOMSIGNALCHARACTERISTICSINTHEFREQUENCYDOMAINANDAIMSATRECOVERINGTHEUSEFULSIGNALINTHEFREQUENCYDOMAINCONTAMINATEDWITHNOISEBASEDONTHELIMITEDOBSERVEDDATATHEDESIGNFOCUSESONTHEALGORITHMANDIMPLEMENTATIONOFTHENONPARAMETERSPECTRALESTIMATIONSFIRST,THEMOSTGENERALSPECTRALESTIMATIONPERIODOGRAMMETHODISDISCUSSEDASTHEPERIODOGRAMMETHODISABIASEDESTIMATEWITHTHELARGEVARIANCEANDLOWRESOLUTION,THEIMPROVEDMETHODSBASEDONTHEPERIODGRAMMETHODNEEDTOBESTUDIEDTHEMAINIMPROVEMENTSINCLUDECHANGINGTHESHAPEOFWINDOWS,AVERAGINGANDSMOOTHINGTHEDATA,ETCASARESULT,SEVERALOTHERNONPARAMETERSPECTRALESTIMATIONALGORITHMSBASEDONTHEPERIODGRAMMETHODSUCHASTHEBTMETHOD,THEAVERAGINGPERIODGRAMMETHOD,INCLUDINGTHEBARTLETTMETHODANDWELCHMETHODANDTHEIRREALIZATIONSBYTHESOFTWAREOFMATLABAREGIVENMOREOVER,THEDIRECTIONOFARRIVALDOAESTIMATIONISTOAPPLYTHEPOWERSPECTRUMESTIMATIONALGORITHMSINTOTHESPACEDOMAIN,IE,DOAISTHENEWDIRECTIONFINDINGTECHNIQUEBASEDONTHEMULTIPLESENSORARRAYCOMBINEDWITHTHEMODERNDIGITALSIGNALPROCESSINGACCORDINGLY,THEREALIZATIONSOFDOAESTIMATIONALGORITHMSBASEDONTHEMETHODOFCAPONANDMUSICAREDISCUSSEDANDCOMPAREDINTHISDESIGNKEYWORDSPERIODGRAMAVERAGINGPERIODGRAMBTSPECTRUMCAPONMETHODMUSICALGORITHMDIRECTIONOFARRIVALDOA目录1绪论111引言1111非参数谱估计1112空间谱估计112本文的主要工作22几种常用的非参数谱估计算法与实现221引言222周期图法3221周期图谱估计3222周期图法的性质3223周期图法在MATLAB中的仿真实现423BT法5231BT法谱估计5232BT谱估计的性质6233BT谱估计在MATLAB中的仿真实现724平均周期图法7241BARTLETT方法及在MATLAB中的仿真实现8242WELCH方法及在MATLAB中的仿真实现1025小结133空间谱估计的算法与实现1431引言1432全向天线的波束下倾14321全向天线的波束下倾的基本原理14322全向天线的波束下倾的仿真实现1533天线的波达方向估计15331天线的波达方向估计的基本原理16332求波达方向的方法1734天线阵的波束形成2035小结234结论23致谢23参考文献24附录24非参数谱估计的算法与实现1绪论11引言111非参数谱估计信号的频谱分析是研究信号特性的重要手段之一。对于确定性信号，可以用傅里叶变换来考察其频谱性质。在通信系统中，最常见的信号往往不是确定信号，而是具有某种统计特性的随机信号。对于随机信号，由于它一般不是周期的，持续时间无限长，具有无限长能量的功率信号，不满足傅里叶变换条件，随机信号也不存在解析表达式，因此，对随机信号，我们不能像确定信号那样进行频谱分析。然而，虽然随机信号的频谱不存在，但如果随机信号是平稳的，其相关函数就是确定的，那么相关函数的傅里叶变换就是它的功率谱密度函数，简称功率谱。通常可以通过求随机信号的功率谱来对这类信号进行频谱分析。功率谱是频率的函数，其物理单位是W/HZ,反映了随机信号各频率成份功率能量的分布情况，可以揭示信号中隐含的周期性及相邻谱峰等有用信息。功率谱估计的应用极其广泛，例如，在语音信号识别、雷达杂波分析、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域，发挥了重要作用。然而，实际应用中的平稳随机信号通常是有限长的，只能根据有限长信号估计原信号的真实功率谱，这就是功率谱估计问题。功率谱估计分为非参数模型谱估计和参数模型谱估计，本文主要研究非参数模型谱估计非参数模型谱估计又称经典谱估计，主要方法有周期图法，平均周期图法，BT法及改进的周期图法，如BARTLETT法和WELCH法。112空间谱估计电磁波到达方向DOA估计在无线电通信、雷达、声纳和导航等领域有着广泛的应用。由于空间信号存在多源可能，加上电磁波传播的多径效应，将会给波束法测向带来较大测量误差，并受干扰信号的影响严重。利用阵列天线进行测向，其空间分辨率因阵列孔径尺寸的选择受经典瑞利限的限制，在瑞利限以内的空间目标则不能分辨。为了突破瑞利限约束，提高角度分辨力，新的测向技术不断涌现。在波束形成技术、零点技术和时域谱估计技术的基础上发展起来的空间谱估计技术是基于阵列信号处理技术的一种超分辨测向技术，用于提高在处理带宽内空问信号角度的估计精度和分辨力。空间谱分析采用了类似时域谱估计中的非线性处理，从而产生了一些特殊的算法，如多重信号分类法MUSIC、旋转不变技术的参数估计法ESPRIT、最小内积法JVW、投影矩阵法和矩阵分解法等。这些方法都巧妙地利用接收信号的相关矩阵的特征结构。其中，MUSIC算法将阵列输出数据的自相关矩阵进行特征分解，然后利用这两个子空间的正交性来估计信号的参数。这一算法的提出开创了空间谱估计算法研究的新时代，成为空间谱估计理论体系中的标志性算法。这种算法也是本文会研究到的。空间谱估计通过阵列信号处理，估计出阵列信号的空间谱，即信号源在空间的分布情况。采用该技术不仅可以提高测向精度和分辨力，而且加快了信号处理的运算速度。12本文的主要工作第一章，扼要介绍了谱估计，非参数谱估计以及空间谱估计等概念。第二章，首先研究周期图法的算法与实现，用MATLAB对其仿真，分析仿真结果，总结周期图法的优缺点，再引出其他几种对周期图法的改进方法，如，BT法，平均周期图法BARTLETT法以及WELCH法。从各种估计法的算法着手，结合MATLAB仿真，分析仿真结果，与周期图法作比较，总结各种算法对周期图法作出了哪些改进。第三章，首先引入什么是空间谱估计，简介其主要应用。采用CAPON算法和MUSIC算法仿真实现了阵列侧向以及波束形成。第四章，总结了本课题研究过程中的主要工作、理论结论。2几种常用的非参数谱估计算法与实现21引言谱估计的非参数方法完全依赖于PSD的定义式KTIER121LIMNITINYE2来实现功率谱的估计，这些方法构成了PSD估计的“经典手段”。本章主要讨论几种非参数谱估计方法、这些方法的性质以及这些方法在MATLAB中的仿真。这里首先介绍一种最常用的谱估计算法，即由2式直接导出的周期图法。对于足够长的数据长度，周期图法能够获得满意的分辨率，但是不是可靠的谱估计算法，这是因为它的方差很大，而且并不随着数据长度的增加而降低。由于周期图法的估计方差很大，促使人们不断地研究新的改进方法，试图以降低分辨率为代价来获得较小的估计方差。目前已经提出了一些改进的方法，我们仅讨论其中最常见的几种。可以证明，在大量数据样本的条件下，所有这些方法在性质和性能方面或多或少是等价的。22周期图法221周期图谱估计一种信号功率谱密度估计方法。它的特点是为得到功率谱估值，先取信号序列的离散傅里叶变换，然后取其幅频特性的平方并除以序列长度，由于序列的离散傅里叶变换具有周期性，因而这种功率谱也具有周期性，常称为周期图。SCHUSTER于1899年首先提出周期图法，也称直接法，取平稳随机信号YT的有限个观察值YO，Y1，YT1对功率谱进行估计，周期图方法依赖于PSD的定义2式。略去只有信号样本信息NITY1的条件下不能实现的期望和极限运算，我们得到21NITIPEY（周期图）3用来确定时间序列中可能存在“隐周期”，是周期图谱估计子3式最早的应用之一，这也许可以认为这种方法的由来。222周期图法的性质分析P的统计特性的重要性在于，证明了周期图作为PSD估计的不可靠性，此外，还提供了一些如何改进周期图以便获得更好的谱估计的内部机理。我们用两个部分来分析周期图方法的性质偏差分析和方差分析。通常用偏差和方差这两种基本度量来刻划估计子的性能，其主要原因是一个估计的总平方误差就是它的偏差的平方与方差之和。为了说明这个问题，设A为任意一个带估计量，A是啊的一个估计。估计的均方差（MSE）为MSE22AEEAEAARE222VRBIS4通过分别考虑MSE的偏差和方差分量，我们可以获得误差来源的深入了解，并得到减少误差的方法。（一）估计的均值WE215其中W是窗函数N的傅里叶变换。当N时，E,是无偏估计，当N为有限值时，是有偏估计，偏差为DWBIAS216（二）估计的方差EDDN20021VAR7其中，0D是矩形窗0ND，1的傅里叶变换，可见，不是的一致估计。随着N增大，谱估计的起伏增大，N时，2VAR。（三）估计的分辨率除了通过偏差和方差来分析周期图法的性质，我们还可以通过分辨率来分析。数据窗为长度为N的矩形窗时，NS2890RE，增大时，分辨率提高，但会使Y的起伏加剧，可见方差与分辨率是一对矛盾。周期图法应用比较广泛，主要是由于它与序列的频谱有对应关系，可以采用FFT快速算法来计算。但是，这种方法需要对无限长的平稳序列进行截断，相当于对其加矩形窗，使之成为有限长数据。同时，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，从而产生频谱泄漏，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，使得谱分辨率较低。223周期图法在MATLAB中的仿真实现仿真实例1用周期图法对信号2130SIN20SINNWX进行功率谱估计，其中NW为白噪声。对信号X的采样率为F1HZ，满足取样定理。程序代码见附录。N0,N1，当N256和N1024时，MATLAB仿真得到的功率谱估计图分别为00102030405060708091402002040三三三HZ三三三三三DB三256三三三三00102030405060708091402002040三三三HZ三三三三三DB三1024三三三三图1周期图法得出的功率谱估计随着采样点数增加，该估计是渐进无偏的。从图1中可以看出，采用周期图法估计得到的功率谱很不平滑，也就是估计的协方差比较大，而且采用增加采样点的方法也不能使周期图变得更加平滑，这是周期图法的缺点。周期图法的优点是能应用离散傅里叶变换的快速算法来进行估值。这种方法适用于长信号序列的情况，在有足够的序列长度时，应用改进的周期图法，可以得到较好的功率谱估值，因而应用很广。23BT法BT法即BLACKMANTURKEY法，这种方法是BLACKMAN与TURKEY两人于1958年提出的，故以他们的名字命名。这一节，我们研究BT法，并与周期图法作比较。231BT法谱估计正如我们已经看到的，即使有很长的样本长度，方差仍然很大时周期图是这种谱估计子存在的主要问题。周期图谱估计较差的统计质量可以直观地解释为是由这样两个方面的原因同时引起的，即C中末端滞后NK的KR估计精度低，以及大量（即便是少量）的协方差估计误差在C中累积相加。C的定义公式C1NKKIER（相关图）8是基于相关的PSD定义1直接导出的相关图谱估计算法。通过在C的定义式8式中截短求和的方法可以减少这两方面的影响。这个思想就导出了BT估计子，1MKKIBTERW9其中N,K是一个偶函数（KW），且10，MKW0，同时W随K平滑衰减到零。由于在式9中对样本协方差序列的滞后进行加权，故被称为滞后窗。因为由这种方法求出的功率谱是通过自相关函数间接得到的，所以称为间接法，又称为自相关法或BT法。当M较小时，式9的计算量不是很大，因此，该方法是在FFT问世之前常用的谱估计方法。当NM时，BT法与周期图法估计的功率谱是一样的；当NM时，BT法的偏差大于周期图法，在窗函数满足一定条件下是渐进无偏估计；方差小于周期图的方差；分辨率比周期图法低，与窗函数选择有关。令W为KW的DTFT,1MKKIKIEE10利用DTFT的性质，我们可以写成KIBTER,01,0MW和,NRR的积的DTFTKDTFKTF11记DTF0,1,0NRRP，可得DWWPBT212由于极大多数常用的窗函数的在0处都有一个相对较窄的主峰。由12式可得BT谱估计子9式就是一个局部加权的平均周期图。232BT谱估计的性质由于0，自然也要求0BT。选择适当的窗函数就能够使估计普拥有这样的性质。如果滞后窗KW式正半定的（即0W），则加窗的协方差序列RK也是正半定的，这意味着对所有的，0BT13即BT谱估计具有非负性。证明既然KW对于是实对称的，那么其DTFW是一个实的偶函数。进一步是一个正半定序列，则对于所有，0。由于定义0P，由12式可知，0W即意味着BT。应该注意的是，有些滞后窗，如三角窗，并不满足13式中的假设，因此，这些窗函数的使用可能导致估计谱出现负值。233BT谱估计在MATLAB中的仿真实现仿真实例2用BT法对信号2130SIN20SINNWX进行功率谱估计，其中NW为白噪声。对信号X的采样率为SF1HZ，满足取样定理，N0,N1。程序代码见附录。当N512时，MATLAB仿真得到的功率谱估计图为00102030405060708091505101520253035三三三HZ三三三三三DB三512三三三三三三图2BT法得出的功率谱估计图通过实验仿真可以直观地看出功率谱估计中BT法和周期图法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。24平均周期图法周期图是信号功率谱的一个有偏估值；而且，当信号序列的长度增大到无穷时，估值的方差不趋于零。因此，随着所取的信号序列长度的不同，所得到的周期图也不同，这种现象称为随机起伏。由于随机起伏大，使用周期图不能得到比较稳定的估值。一些学者对此作了改进。为了减小随机起伏，MSBARTLETT特提出平均周期图法，即先把信号序列分为若干段，对每段分别计算其周期图，然后取各个周期图的平均作为功率谱的估值。平均周期图可以减小随机起伏，但是，如果信号序列不是足够长，由于每段序列长度变短，功率谱估值对不同频率成分的分辨能力也随之下降。另一种改进方法是将周期图与一个适当的频域窗函数相褶积，从而对周期图产生平滑作用，以减小随机起伏。加窗处理的结果虽然可以使随机起伏减小，但也会使周期图的分辨能力下降。POWELCH提出一种把加窗处理与平均处理结合起来的方法。先把分段的数据乘以窗函数进行加窗处理，分别计算其周期图，然后进行平均。WELCH方法是较常用的一种计算方法。为了得到较好的功率谱估值，加窗和平均处理均应兼顾减小随机起伏和保证有足够的谱分辨率两个方面。241BARTLETT方法及在MATLAB中的仿真实现（一）BARTLETT方法BARTLETT方法的基本思想很简单为了减少周期图的起伏，把长度为N的观测样本分割为长度为M的NL/个子样本，然后对从这些子样本得到的周期图对于每一个取平均。将BARTLETT方法表示成如下的数学形式。令,1TJYTJLJM,114表示第J个子样本的观测值，并令21MITIJJEY15表示对应的周期图。BARTLETT谱估计由下式给出LJB116由于BARTLETT方法是在长度为M的数据段上计算的，故其分辨率应该约为M/1。因此，与与原来周期图方法比较，BARTLETT方法的分辨率降低了L倍。作为对分辨率降低的回报，可以预料，BARTLETT方法的方差减小了。事实上，可以发现BARTLETT方法把周期图的方差减小了L倍。在选择M（或L）时，分辨率和方差之间的折中是很明显的。（二）BARTLETT方法的性质正如我们所知道的，16式中的J可以写成1MMKKIJJER17其中JR是第J个子样本的协方差序列。把17式代入16式得KIMKLJBER1118我们发现B在形式上与使用矩形窗的BT估计很相似。子样本协方差KRJ在J上的平均是ACSKR的一个估计。但是，在18式中ACS的估计并没有有效地利用可以得到的数据滞后的乘积KTY，特别是当K接近1M时。实际上，对于1MK，只有大约可以滞后乘积的/1被用来形成18式中ACS的估计。可以预料，这些滞后的方差要大于BT估计中对应的KR滞后，类似地，B方差要比BT的方差大。此外，由于BARTLETT方法使用了一个固定的矩形滞后窗，所以其分辨率泄漏折中的灵活性也不如BT法。由于上述原因，我们总结为在形式上，1416式定义的BARTLETT估计与使用长度为M的固定矩形滞后窗的BT估计类似，但其方差略大于后者。由此和以前讨论的BT谱估计的性质可知，与基本的周期图方法比较，BARTLETT估计降低了分辨率，减少了方差（其倍数均为MNL/）。矩形窗的主瓣要比极大多数其他滞后窗的主瓣窄。所以，在BT类谱估计中，BARTLETT估计可望得到最少的谱峰兼并（即有最大的分辨率），但有最明显的泄漏。（三）BARTLETT估计的MATLAB仿真仿真实例3用BARTLETT估计对信号TNTTX102SI602SIN进行功率谱估计，其中TN为高斯白噪声。对信号的采样频率为SF500HZ，满足抽样定理。N500，源程序代码见附录。得到的功率谱估计图如下BARTLETT估计对周期图法德改进的思想是将信号分段进行估计，并将这些估计结果进行平均而减少估计的方差，使得估计功率谱图变得平滑。如上例中将501点的信号分为三段，分别作周期图法估计，然后加以平均。功率谱估计结果如图3所示，将信号分段得到的估计值显然平滑了许多。0501001502002503003504004505006040200三三三HZ三三三三三DB三三三三三三三三500三三三三0501001502002503003504004505001020304050三三三HZ三三三三三DB三三500三三三三三三三三三三三三图3BARTLETT法频率谱估计图242WELCH方法及在MATLAB中的仿真实现一WELCH谱估计WELCH方法是BARTLETT方法经过两个方面的改进得到的。首先，WELCH方法中的数据段允许重叠。其次，在计算在周期图前，对每一个数据段加窗。为了用数学方式描述WELCH方法，令TKJYTJ1，SJM,119表示第J个数据段。在19式中，K是第J个序列起始点。若MK，则序列间没有重叠（但是都是邻接的），我们得到用于BARTLETT方法中的样本分割（得到MNLS/个数据子样本）。但是，在WELCH方法中的推荐值为2/，在这种情况下，得到NS/2个数据段（前后数据段之间重叠50）。计算与TYJ对应的加窗周期图21MITIJJEYTP20其中P表示时间窗T的功率21MIT21通过对20式中的加窗周期图取平均得到PSD的WELCH估计SJW122要解释上述对BARTLETT方法的改进，从而导出WELCH方法的理由很简单。通过允许数据段之间的重叠，就有更多的周期图参与22式中的平均，故可以降低PSD估计的方差。在周期图计算中引入窗函数，就能够更有效地控制PSD估计中的偏差/分辨率性质。除此之外，时间窗还可用来对位于没个子样本末端的数据给予较小的加权，这样，尽管子样本间是重叠的，但可以使相邻的子样本间的相关性减弱。通过22式中平均。这种去相关的重要作用能够更有效地降低方差。（二）WELCH谱估计的MATLAB仿真仿真实例4用WELCH谱估计对信号TNTTX102SI602SIN进行功率谱估计，其中TN为高斯白噪声。对信号的采样频率为SF500HZ，满足抽样定理。N500，源程序代码见附录。得到的功率谱估计图如下0501001502002503003504004505006040200三三三HZ三三三三三DB三三三三三三三三500三三三三050100150200250300350400450500806040200三三三HZ三三三三三DB三三500三三三三三三三三三三三三三三三三三三三三图4采用样本混叠增加分段数得到的功率谱050100150200250300350400450500806040200三三三HZ三三三三三DB三三三三三三三三500三三三三050100150200250300350400450500806040200三三三HZ三三三三三DB三三三WELCH三三三三三三三三三三三三三三三三图5WELCH平均修正周期图法得到的功率谱估计图增加分段数可以进一步降低估计的方差，然而若每段中的数据点数太少，就会使得估计的频率分辨率下降很多。在样本信号总点数一定的条件下，我们可以采用使分段相互重叠的方法来增加分段数，而保持每段中信号点数不变，这样就在保证频率分辨率的前提下进一步降低估计的协方差。如F仿真实例4所示，将XN分为5段，而且每段之间有重叠，对每段分别求FFT后求平均值，得出功率估计值。功率谱估计结果如图4所示，与图3相比，显然估计值又平滑了许多。另外，采用加窗法也可以降低估计的协方差，这也是窗函数应用的一个方面。即在计算周期图法之前，对数据分段并加非矩形窗（如海明窗、汉宁窗或凯瑟窗等），然后再采用分段长度一半的混叠率，就能够大大降低估计方差。这种方法称为WELCH平均修正周期图法，简称WELCH法。仿真实例4是采用海明窗的WELCH法。谱估计结果如图5所示，与图3，图4相比，显然WELCH法得出的功率谱估计显然平滑性较好。（三）不同窗函数的WELCH谱估计的MATLAB仿真仿真实例5用WELCH方法对信号2130SIN20SINNWX进行功率谱估计，其中NW为白噪声。对信号X的采样率为F1HZ，满足取样定理，N0,N1。程序代码见附录。当N1000时，窗函数采用矩形窗、HANNING窗、HAMMING窗，MATLAB仿真得到的功率谱估计图分别为0010203040506070809150050三三三HZ三三三三三DB三三三三0010203040506070809102040三三三HZ三三三三三DB三HANNING三0010203040506070809102040三三三HZ三三三三三DB三HAMMING三图6WELCH法加不同的窗函数得出的功率谱估计图仿真结果分析由图6可以看出，使用不同的窗函数谱估计的质量是不一样的，矩形窗的主瓣较窄，分辨率较好，但方差较大，噪声水平较高；而HANING窗和HAMMING窗的主瓣较宽，分辨率较低，但方差较小，噪声水平较低。因此，在进行谱分析时选择何种窗函数，要视具体情况而定。如果强调高分辨率，能精确读出主瓣频率，而不关心幅度的精度，例如测量震动物体的自震频率，可以选用主瓣宽度比较窄的矩形窗；对受到强干扰的窄带信号，若干扰靠近信号，则可选用旁瓣幅度较小的窗函数，若离开通带较远，则可选用渐近线衰减速度比较快的窗函数。总之，要针对不同的信号和不同的处理目的来选择合适的窗函数，这样才能得到良好的效果。25小结本章先研究周期图谱估计，从均值，方差，分辨率讨论周期图法的性能。估计的均值得出当N时，是无偏估计，当N为有限值时，是有偏估计。估计的方差得出不是的一致估计。随着N增大，谱估计的起伏增大。N增大时，分辨率提高，但会使Y的起伏加剧，可见方差与分辨率是一对矛盾。为此，提出了改进周期图法的几种途径，即改善窗口形状、平均和平滑。BT法中涉及到两次加窗一次是用于截取数据，叫数据窗；一次是用于截取相关函数，叫延迟窗，延迟窗常用三角窗。相关法的两次加窗将造成泄漏效应，影响到谱峰分辨率和方差。平均就是将截取的数据段再分成L个小段，分别计算功率谱后取功率谱的平均，这种方法使估计的方差减少，但偏差加大，分辨率下降。平滑是用一个适当的窗函数与算出的功率谱进行卷积，使谱线平滑。这种方法得出的谱估计是无偏的，方差也小，但分辨率下降。现在比较常用的改进方法是WELCH法，又叫加权交叠平均法，简记为WOSA法，这种方法以加窗（加权）求取平滑，以分段重叠求得平均，因此集平均与平滑的优点于一体，同时也不可避免带有两者的缺点，因此归根到底是一种折中。WELCH法得出的是无偏谱估计，当段数L增大时方差减小，趋于一致估计，分段平均减小了由数据样本本身的随机性带来的方差，段数越多方差越小，但分辨率下降。非参数谱估计法的基础是FFT，它们出现较早，运算量较小，物理概念明确，便于工程实现，仍是目前常用的谱估计方法，有的已被固化到FFT仪和谱分析仪中。但有一些难以克服的缺点谱的分辨率较低，它正比于（N是数据长度）。加窗的坏影响不可避免。较宽的主瓣降低分辨力，较大的旁瓣有可能掩盖真实谱中较弱的成分，或是产生假的峰值。而没有一个窗函数能使谱估计在方差、偏差和分辨率各方面同时改善的，使用窗函数只是一种折中的技巧，不是解决问题的根本办法。不是真实谱的一致估计，且N增大时谱曲线起伏加剧。基于上述原因，非参数谱估计还有继续研究的必要性。3空间谱估计的算法与实现31引言本章将讨论使用一个M元无源阵元组成的阵列对N个辐射源进行定位的问题。这里，从辐射源来的能量可能是声音、电磁波等等，接收传感器将该能量转换为电信号。接收阵元可能为电磁天线、水声传感器或者地震检波器。这类问题在雷达或声纳系统、通信、天体物理学、生物医学、地震水下监视以及其他许多领域具有广泛的应用。这一问题的基本内容包括确定能量在空间（空气、水下或地下）的分布，这里能量高度集中的点就是辐射源。这类问题称为空间谱估计。之所以称为谱估计还因为辐射源定位问题和前面讨论的谱估计有着非常紧密的联系。事实上我们将看到，前面所述的方法几乎都可以用作推导辐射源定位问题的解。32全向天线的波束下倾321全向天线的波束下倾的基本原理广播用的发射机天线建在天线塔上，希望覆盖广大的地域。在覆盖区的边沿，由于距离远，地球表面弯曲，电磁波场强衰减很快，因此为了有效地覆盖既定的区域，保证区域内场强不低于特定值，通常采用波束下倾的方法，将轴向排列的半波振子天线，通过调节天线的轴向间距、馈电的相位，使得轴向天线阵列的方向图实现波束下倾。轴向排列的半振子天线结构，阵元之间的距离是D，垂直轴线与电磁波辐射方向的夹角是，相邻阵元馈入信号的相位差是，相邻阵元发出电磁波到达同一地点的程差是COSD，相邻阵元发出的电磁波到达同一地点的相位差参看公式（23）2COS1120DEEENJJJ（23）322全向天线的波束下倾的仿真实现根据以上讨论，编写出绘制轴向排列天线阵列的方向图程序。源程序代码见附录。仿真结果如下图24683021060240902701203001503301800图7轴向阵列波束下倾仿真图33天线的波达方向估计天线阵列信号处理作为信号处理的一个重要分支，在通信、雷达、地震勘测、射电天文等领域获得了广泛应用和迅速发展。它包括信源定位确定阵列到信源的仰角和方位角，甚至距离（若信源位于近场）信源分离确定各个信源发射的信号波形。各个信源从不同方向到达阵列，即使它们在时域和频域上是叠加的。信道估计确定信源与阵列之间的传输信道的参数（多径参数）。阵列信号处理时改善蜂窝和个人通信服务系统质量、覆盖范围和容量的一个强有力的工具。建立在波达方向估计、波束形成基础上的智能天线的应用，抑制了干扰信号，改善了欲接收信号的信噪比，降低了数字通信系统的误码率。将接收天线阵列用于公众网的反向连接（客户到基站）时，多个接收天线能够收集更多的信号能量，若天线在空间足够分离或极化各异，则多个天线能够提供很好的分集接收，并抑制多径传输引起的衰落。这些好处可以扩大基站的覆盖范围，改善通信质量。331天线的波达方向估计的基本原理全向天线不仅利用率不高，而且对各种信号不加区别地接收，降低了通信质量。定点无线通信采用定向天线，大幅度地改善了通信质量。面对众多移动用户的公众通信网基站和专用移动通信网，采用天线指向即波束可变的天线（智能天线），可以使移动通信的质量得到很大的改善。为使天线的波束指向可控，甚至形状可控，采用阵列天线是合适的。在距离通信源足够远的空间里，可以将到达的电磁波视为平面波。对于等距离直线阵天线，由于调制在载波上的基带信号码元宽度与波速的乘积远大于天线阵的尺寸，因此多个天线阵元上的信号的幅度可视为不变，而它们的载波的相位差则取决于其相互位置、尺寸、波长和到达方向。无线接收的无线信号中有许多成分，其中我们关心的是S信号。天线阵列各个阵元接收的电磁波信号因为阵元排列位置的不同带来相位差。经过特定参数的加权控制器W处理后，进一步改变了各个阵元输出信号的相位和幅度。处理的目标是使得阵元输出信号的和天线输出Y中的S成分具有最大输出。用S信号作为基准信号，反馈控制单元的功能就是将输出信号Y与基准信号S的差值（即误差信号），作为调节控制加权控制器W参数的依据。反馈控制是使减小，Y中的S成分加大，也就是说，天线阵列接收方向图指向了信号的方向。以下通过两个简单的例子介绍智能天线中波达方向估计和波束形式的原理。等距离直线阵智能天线，我们把天线阵元沿X轴排列，从0到1M。若有一平面波以角（入射线与Z轴的夹角）和角（入射线与轴的夹角）入射到阵列上，第K号阵元上产生的信号为KX，它与0号阵元的相位差为COSIN2DA（24）式中，和D分别是入射波的波长和阵元的间距，KA亦称阵因子。计入阵因子的影响，第K号阵元的输出时KXA，即U。为了使天线阵的输出满足需要，在每个阵元上，用加权因子W进行控制。这样第号阵元上输出的信号为KXAW，即KUW。若到达天线阵的信号是N个，则天线阵的输出是N个信号在M个阵元上的输出的叠加。将问题简化为XY平面的二维问题（1SIN），并用解析式表达如下NXNXXN,21（25）111COS2EXPSDMJJA（26）其中，1为第一个信号的入射角。NA,21（27）MW（28）UAXTYHH（29）其中，HMUU110,。求多个信号到达的方向（波达方向）的方法有多种。下面讨论其中两种方法（CAPON法和MUSIC算法）及其仿真验证结果。332求波达方向的方法1CAPON法CAPON法亦称最小方差无畸变响应MVDR（MINIMUMVARIANCEDISTORTIONLESSRESPONSE）。天线阵列中的阵元数决定了阵列方向图设计中的自由度数。CAPON法将阵列中可控的自由度用来形成期望的波束形状，达到对有用信号进行提升和对无用信号进行抑制的目的，并将其优化问题表达为WRNYEUHWWMII2（30）其约束条件为10AWH可以证明上式的解为ARUHCAP1代入式（29），可以得到相应到的功率为ARPUHCAP1（31）2MUSIC法MUSIC法亦称多重信号分类（MULTIPLESIGNALCLASSIFICATION）。当入射信号数N小于阵元数M时，的UR矩阵有与进入阵列的信号数目相等的非零特征值及个为零的特征值。1NNVV2MV（NV是U相应的噪声特征矢量）因为A与的正交性，分母很小，峰值很大，这样可以得出MUSIC法的空间谱为APHMUSIC（32）其中，HNV称为噪声子空间的正交投影估计。3仿真以下是对七单元线性天线阵在四信号输入情况下了仿真仿真实例6设信号1从4/方向入射，信号2从3/方向入射，信号3从6/方向入射，信号4从/3方向入射。下面是用两种方法，分别用两种坐标（直角坐标和极坐标）求波达方向估计。程序代码见附录。仿真结果如下49999999999999998000000000000009999999999999999600000000000000149999999999999990000000000000003021060240902701203001503301800图8CAPON法作出的波达方向估计（极坐标）020406080100120140160180100105101010151020102510301035图9MUSIC法作出的波达方向估计（直角坐标）500000100000015000003021060240902701203001503301800图10MUSIC法作出的波达方向估计（极坐标）020406080100120140160180101510101051001051010图11CAPON法作出的波达方向估计（直角坐标）从对以上四图的分析可以得到如下结论两种方法估计得都比较准确。在4/3,/,6处有尖锐的方向图线。MUSIC法方向图像幅度更大。34天线阵的波束形成1原理我们以等距离圆阵为例来讨论天线阵的波束形成。我们取圆阵的圆心为O，在阵列圆上任取一点B，把天线阵元顺序定为0到1M。若有一平面波以角入射到阵列上，第K号阵元上产生的信号为KX，它与到达阵元中心的波前的相位差是KKRJACOS2EXP；与R分别是入射波的波长与阵列圆的半径，亦称阵因子。为了使天线阵的输出满足需要，在每个阵元上加上加权因子K控制。这样第号阵元上输出的信号为KXA。若到达天线阵的信号是N个，天线阵输出的是N个信号的M个阵元上的输出的叠加。用解析表达式如下COS2EXP1KKRJA（33）HMKA,NNXNXX110（34）式中，K是阵元以OB为基准顺时针画出的角度。NA,21（35）HMWW10（36）UWAXTYH（37）为了求得多个信号到达的方向（波达方向），可以采用上述的CAPON、MUSIC两种方法）。波束形成可以采用下面的方法当有多（N）个信号输入时，其中1个信号时我们关心的，1N个信号是需要抑制的。方程组（38）描述了上述需求的约束条件（四个信号输入，第一个信号是我们关心的，其余信号是需要抑制的）。THUW0,1（38）根据信号波达方向的知识AX及约束条件求解方程组（38），可以得到110,MW（39）代入式（37）可得到阵列输出的方向特性。2仿真下面是四个输入信号八单元圆阵列天线阵的波达方向估计和波束形成。仿真实例7设信号1从3/2方向入射，信号2从3/方向入射，信号3从2/3方向入射，信号4从685方向入射。下面是用两种算法求波达方向估计及波束形成。源程序代码见附录。仿真结果如下050100150200250300350400100105101010151020102510301035图12CAPON算法波达方向估计050100150200250300350400101510101051001051010图13MUSIC算法波达方向估计051153021060240902701203001503301800图14波束形成35小结本章在对非参数谱估计的理解基础上，讨论了其在空间领域进行推广，具体就是了采用CAPON算法和MUSIC算法求天线阵的波达方向，讨论了天线阵的波束形成，及其MATLAB的仿真实现。4结论本设计的主要任务是研究非参数谱估计的算法与实现，另外对功率谱估计算法在空间领域的推广应用尤其是阵列侧向问题进行了讨论。在经典谱估计中，无论是周期图法还是其改进的方法，都存在着频率分辨率低、方差性能不好的问题，原因是谱估计时需要对数据加窗截断，用有限个数据或其自相关函数来估计无限个数据的功率谱，这其实是假定了窗以外的数据或自相关函数全为零，这种假定是不符合实际的，正是由于这些不符合实际的假设造成了经典谱估计分辨率较差。另外，经典谱估计的功率谱定义中既无求均值运算又无求极限运算，因而使得谱估计的方差性能较差，当数据很短时，这个问题更为突出。如何选取最佳窗函数、提高频谱分辨率，如何在短数据情况下提高信号谱估计质量，还需要进一步研究。CAPON法和MUSIC算法可以用于空间谱估计测向，充分利用了天线各个阵元从空间电磁场接收到的全部信息，而传统测向方法仅仅利用了其中一部分信息，因而具有多方面的优越性1可以实现对同一一信道中存在的多个信号同时测向；2超分辨率测向，即对处于天线阵固有波束宽度以内的两个以上方向的来波同时测向；3阵元位置可以随意放置，各个阵元的方向特性也可以是任意的；4具有优良的测向灵敏度和准确度；5可以较好地克服测向模糊问题。参考文献1PETRESTOICA,RANDOLPHLMOSESINTRODUCTIONTOSPECTRALANALYSIS,PRENTICEHALL,INC19972魏平，唐斌，陆明泉，唐海燕，谈觅舟，肖先赐谱估计原理J电子科技大学研究生教材，19993胡广书数字信号处理一理论、算法与实现M北京清华大学出版社，20034王晓峰，王炳和周期图及其改进方法中谱分析率的MATLAB分析J武警工程学院学报，2003675775姚武川，姚天任经典谱估计方法的MATLAB分析J华中理工大学学报，2000，28445476刘刚，吕新华，攸阳阵列信号处理中基于MUSIC算法的空间谱估计J微计算机信息，2006，2243302303附录附录直接法CLEARALLCLCN256FS1N256N0N1XNSQRT20SIN2PI02NSQRT2SIN2PI0213NRANDNSIZENNFFT256PXX10LOG10ABSFFTXN,NFFT2/NF0LENGTHPXX1FS/LENGTHPXXSUBPLOT211PLOTF,PXXXLABEL频率（HZ）YLABEL功率谱（DB）TITLE256点功率谱GRIDN1024FS1N1024N0N1XNSQRT20SIN2PI02NSQRT2SIN2PI0213NRANDNSIZENNFFT1024PXX10LOG10ABSFFTXN,NFFT2/NF0LENGTHPXX1FS/LENGTHPXXSUBPLOT212PLOTF,PXXXLABEL频率（HZ）YLABEL功率谱（DB）TITLE1024点功率谱GRID附录间接法，512点功率谱估计CLEARALLCLCFS1N512N0N1XNSQRT20SIN2PI02NSQRT2SIN2PI0213NRANDNSIZENNFFT512PXX10LOG10FFTXCORRXN,UNBIASED,NFFTF0LENGTHPXX1FS/LENGTHPXXPLOTF,PXXXLABEL频率（HZ）YLABEL功率谱（DB）TITLE512点功率谱估计GRID附录BARTLETTCLEARALLCLCN500FS500T01/FS1F01FSXNSIN2PI60T2SIN2PI110TRANDN1,LENGTHTPXX1ABS1/FSFFTXN2PXX2ABSFFTXN11672ABSFFTXN1683342ABSFFTXN3355012/3SUBPLOT211PLOTF,10LOG10PXX1XLABEL频率（HZ）YLABEL功率谱（DB）TITLE周期图法得到的500点功率谱SUBPLOT212PLOT03FS,10LOG10PXX2XLABEL频率（HZ）YLABEL功率谱（DB）TITLE将500点信号分三段得到的功率谱附录WELCHCLEARALLCLCN500FS500T01/FS1F01FSXNSIN2PI60T2SIN2PI110TRANDN1,LENGTHTPXX1ABS1/FSFFTXN2PXX3ABS3/FSFFTXN11672ABS3/FSFFTXN832492ABS3/FSFFTXN1683342ABS3/FSFFTXN2504162ABS3/FSFFTXN3355012/5SUBPLOT211PLOTF,10LOG10PXX1XLABEL频