非参数谱估计算法与实现.doc

非参数谱估计的算法与实现 内容摘要 功率谱估计是随机信号分析中的一个重要内容 主要研究信号在频域中的 各种特征 目的是根据有限 观测数据在频域内提取被噪声 污染的有用信号 本设计针对非参数谱估计的算法与实现展开研究 首先讨论了一种最常 见的谱估计法 周期图法 由于周期图法是有偏估计 其方差过大 分辨率 较低 因此 需要研究基于周期图的各种改进方法 主要的改进途径包括改善 窗口形状 对数据进行平均和平滑等 为此 讨论了基于周期图的其他几种非 参数谱估计算法 即 BT 法 平均周期图法 包括 Bartlett 法和 Welch 法 及 其通过 MATLAB 的仿真实现 此外 信号来波方向估计是功率谱估计算法在空间领域的推广应用 具体 就是以多元天线阵结合现代数字信号处理为基础的新型测向技术 为此 本设 计还讨论比较了采用 Capon 法 MUSIC 算法实现波达方向估计的实现 关键词 周期图法 平均周期图法 BT谱估计 Capon法 Music算法 波达方向 Algorithms and Implementations of Non Parameter Spectral Estimations Abstract Power spectrum estimate is an important part of the random signal analysis which primarily studies the random signal characteristics in the frequency domain and aims at recovering the useful signal in the frequency domain contaminated with noise based on the limited observed data The design focuses on the algorithm and implementation of the non parameter spectral estimations First the most general spectral estimation periodogram method is discussed As the periodogram method is a biased estimate with the large variance and low resolution the improved methods based on the periodgram method need to be studied The main improvements include changing the shape of windows averaging and smoothing the data etc As a result several other non parameter spectral estimation algorithms based on the periodgram method such as the BT method the averaging periodgram method including the Bartlett method and Welch method and their realizations by the software of MATLAB are given Moreover the direction of arrival DOA estimation is to apply the power spectrum estimation algorithms into the space domain i e DOA is the new direction finding technique based on the multiple sensor array combined with the modern digital signal processing Accordingly the realizations of DOA estimation algorithms based on the method of Capon and MUSIC are discussed and compared in this design Key words Periodgram Averaging Periodgram BT Spectrum Capon method Music Algorithm Direction of arrival DOA 目 录 1 绪 论 1 1 1 引言 1 1 1 1 非参数谱估计 1 1 1 2 空间谱估计 1 1 2 本文的主要工作 2 2 几种常用的非参数谱估计算法与实现 2 2 1 引言 2 2 2 周期图法 3 2 2 1 周期图谱估计 3 2 2 2 周期图法的性质 3 2 2 3 周期图法在 MATLAB 中的仿真实现 4 2 3 BT 法 5 2 3 1 BT 法谱估计 5 2 3 2 BT 谱估计的性质 6 2 3 3 BT 谱估计在 MATLAB 中的仿真实现 7 2 4 平均周期图法 7 2 4 1 Bartlett 方法及在 MATLAB 中的仿真实现 8 2 4 2 Welch 方法及在 MATLAB 中的仿真实现 10 2 5 小结 13 3 空间谱估计的算法与实现 14 3 1 引言 14 3 2 全向天线的波束下倾 14 3 2 1 全向天线的波束下倾的基本原理 14 3 2 2 全向天线的波束下倾的仿真实现 15 3 3 天线的波达方向估计 15 3 3 1 天线的波达方向估计的基本原理 16 3 3 2 求波达方向的方法 17 3 4 天线阵的波束形成 20 3 5 小结 23 4 结 论 23 致谢 23 参考文献 24 附录 24 0 非参数谱估计的算法与实现 1 绪 论 1 1 引言 1 1 1 非参数谱估计 信号的频谱分析是研究信号特性的重要手段之一 对于确定性信号 可以 用傅里叶变换来考察其频谱性质 在通信系统中 最常见的信号往往不是确定 信号 而是具有某种统计特性的随机信号 对于随机信号 由于它一般不是周 期的 持续时间无限长 具有无限长能量的功率信号 不满足傅里叶变换条件 随机信号也不存在解析表达式 因此 对随机信号 我们不能像确定信号那样 进行频谱分析 然而 虽然随机信号的频谱不存在 但如果随机信号是平稳的 其相关函数就是确定的 那么相关函数的傅里叶变换就是它的功率谱密度函数 简称功率谱 通常可以通过求随机信号的功率谱来对这类信号进行频谱分析 功率谱是频率的函数 其物理单位是W Hz 反映了随机信号各频率成份功率 能量的分布情况 可以揭示信号中隐含的周期性及相邻谱峰等有用信息 功率 谱估计的应用极其广泛 例如 在语音信号识别 雷达杂波分析 地震勘探信 号处理 水声信号处理 系统辨识中非线性系统识别 物理光学中透镜干涉 流体力学的内波分析 太阳黑子活动周期研究等许多领域 发挥了重要作用 然而 实际应用中的平稳随机信号通常是有限长的 只能根据有限长信号估计 原信号的真实功率谱 这就是功率谱估计问题 功率谱估计分为非参数模型谱估计和参数模型谱估计 本文主要研究非参 数模型谱估计 非参数模型谱估计又称经典谱估计 主要方法有 周期图法 平均周期图法 BT法及改进的周期图法 如Bartlett法和Welch法 1 1 2 空间谱估计 电磁波到达方向 DOA 估计在无线电通信 雷达 声纳和导航等领域有着广 泛的应用 由于空间信号存在多源可能 加上电磁波传播的多径效应 将会给 波束法测向带来较大测量误差 并受干扰信号的影响严重 利用阵列天线进行 测向 其空间分辨率因阵列孔径尺寸的选择受经典瑞利限的限制 在瑞利限以 内的空间目标则不能分辨 为了突破瑞利限约束 提高角度分辨力 新的测向技术不断涌现 在波束 形成技术 零点技术和时域谱估计技术的基础上发展起来的空间谱估计技术是 基于阵列信号处理技术的一种超分辨测向技术 用于提高在处理带宽内空问信 号角度的估计精度和分辨力 空间谱分析采用了类似时域谱估计中的非线性处理 从而产生了一些特殊 1 的算法 如多重信号分类法 MUSIC 旋转不变技术的参数估计法 ESPRIT 最 小内积法 JVW 投影矩阵法和矩阵分解法等 这些方法都巧妙地利用接收信号 的相关矩阵的特征结构 其中 MUSIC算法将阵列输出数据的自相关矩阵进行特 征分解 然后利用这两个子空间的正交性来估计信号的参数 这一算法的提出 开创了空间谱估计算法研究的新时代 成为空间谱估计理论体系中的标志性算 法 这种算法也是本文会研究到的 空间谱估计通过阵列信号处理 估计出阵 列信号的空间谱 即信号源在空间的分布情况 采用该技术不仅可以提高测向 精度和分辨力 而且加快了信号处理的运算速度 1 2 本文的主要工作 第一章 扼要介绍了谱估计 非参数谱估计以及空间谱估计等概念 第二章 首先研究周期图法的算法与实现 用MATLAB对其仿真 分析仿真 结果 总结周期图法的优缺点 再引出其他几种对周期图法的改进方法 如 BT法 平均周期图法 Bartlett法以及Welch法 从各种估计法的算法着手 结 合MATLAB仿真 分析仿真结果 与周期图法作比较 总结各种算法对周期图法 作出了哪些改进 第三章 首先引入什么是空间谱估计 简介其主要应用 采用Capon算法 和MUSIC算法仿真实现了阵列侧向以及波束形成 第四章 总结了本课题研究过程中的主要工作 理论结论 2 几种常用的非参数谱估计算法与实现 2 1 引言 谱估计的非参数方法完全依赖于PSD的定义式 1 k ti ekr 2 2 1 1 lim N i ti N ety N E 来实现功率谱的估计 这些方法构成了PSD估计的 经典手段 本章主要讨论 几种非参数谱估计方法 这些方法的性质以及这些方法在MATLAB中的仿真 这里首先介绍一种最常用的谱估计算法 即由 2 式直接导出的周期图法 对于足够长的数据长度 周期图法能够获得满意的分辨率 但是不是可靠的谱 估计算法 这是因为它的方差很大 而且并不随着数据长度的增加而降低 由于周期图法的估计方差很大 促使人们不断地研究新的改进方法 试图 以降低分辨率为代价来获得较小的估计方差 目前已经提出了一些改进的方法 我们仅讨论其中最常见的几种 可以证明 在大量数据样本的条件下 所有这 些方法在性质和性能方面或多或少是等价的 2 2 2 周期图法 2 2 1 周期图谱估计 一种信号功率谱密度估计方法 它的特点是 为得到功率谱估值 先取 信号序列的离散傅里叶变换 然后取其幅频特性的平方并除以序列长度 由于序列的离散傅里叶变换具有周期性 因而这种功率谱也具有周期性 常 称为周期图 Schuster于1899年首先提出周期图法 也称直接法 取平稳随机信号Y t 的有限个观察值 y O y 1 y T 1 对功率谱进行估计 周期图方 法依赖于PSD的定义 2 式 略去只有信号样本信息 的条件下不能 N i ty 1 实现的期望和极限运算 我们得到 周期图 2 1 1 N i ti p ety N 3 用来确定时间序列中可能存在 隐周期 是周期图谱估计子 3 式最早的 应用之一 这也许可以认为这种方法的由来 2 2 2 周期图法的性质 分析 的统计特性的重要性在于 证明了周期图作为 PSD 估计的不 p 可靠性 此外 还提供了一些如何改进周期图以便获得更好的谱估计的内部机 理 我们用两个部分来分析周期图方法的性质 偏差分析和方差分析 通常用偏差和方差这两种基本度量来刻划估计子的性能 其主要原因是一 个估计的总平方误差就是它的偏差的平方与方差之和 为了说明这个问题 设 a 为任意一个带估计量 是啊的一个估计 估计的均方差 MSE 为 a MSE 2 2 aaEaEaEaaE aaEaEaEaaEaEaERe 22 4 2 var abiasa 通过分别考虑 MSE 的偏差和方差分量 我们可以获得误差来源的深入了解 并 得到减少误差的方法 3 一 估计的均值 5 WE 2 1 其中是窗函数的傅里叶变换 当时 W n N E 是无偏估计 当 N 为有限值时 是有偏估计 偏差为 6 dWbias 2 1 二 估计的方差 7 EdDD N 2 00 2 1 var 其中 是矩形窗 的傅里叶变换 可见 不 0 D 1 0 nd1 Nn 是的一致估计 随着 N 增大 谱估计的起伏增大 时 N 2 var 三 估计的分辨率 除了通过偏差和方差来分析周期图法的性质 我们还可以通过分辨率来 分析 数据窗为长度为 的矩形窗时 增大时 分辨N N s 2 89 0Re N 率提高 但会使的起伏加剧 可见方差与分辨率是一对矛盾 y 周期图法应用比较广泛 主要是由于它与序列的频谱有对应关系 可以采 用FFT快速算法来计算 但是 这种方法需要对无限长的平稳序列进行截断 相 当于对其加矩形窗 使之成为有限长数据 同时 这也意味着对自相关函数加 三角窗 使功率谱与窗函数卷积 从而产生频谱泄漏 容易使弱信号的主瓣被 强信号的旁瓣所淹没 造成频谱的模糊和失真 使得谱分辨率较低 2 2 3 周期图法在 MATLAB 中的仿真实现 仿真实例仿真实例 1 1 用周期图法对信号 进行功率谱估计 其中为白 213 0 2sin 2 2 02sin 20 nwnnnx nw 噪声 对信号的采样率为 1Hz 满足取样定理 程序代码见附录 nx s f n 0 N 1 当N 256和N 1024时 MATLAB仿真得到的功率谱估计图分别为 4 00 10 20 30 40 50 60 70 80 91 40 20 0 20 40 三 三 三 Hz三 三 三 三 三 dB三 256三 三 三 三 00 10 20 30 40 50 60 70 80 91 40 20 0 20 40 三 三 三 Hz三 三 三 三 三 dB三 1024三 三 三 三 图 1 周期图法得出的功率谱估计 随着采样点数增加 该估计是渐进无偏的 从图1中可以看出 采用周期图 法估计得到的功率谱很不平滑 也就是估计的协方差比较大 而且采用增加采 样点的方法也不能使周期图变得更加平滑 这是周期图法的缺点 周期图法的优点是能应用离散傅里叶变换的快速算法来进行估值 这 种方法适用于长信号序列的情况 在有足够的序列长度时 应用改进的周期 图法 可以得到较好的功率谱估值 因而应用很广 2 3 BT 法 BT法即Blackman Turkey法 这种方法是Blackman与Turkey两人于1958年提 出的 故以他们的名字命名 这一节 我们研究BT法 并与周期图法作比较 2 3 1 BT 法谱估计 正如我们已经看到的 即使有很长的样本长度 方差仍然很大时周期图是这 种谱估计子存在的主要问题 周期图谱估计较差的统计质量可以直观地解释为 是由这样两个方面的原因同时引起的 即中末端滞后的估计 c Nk kr 精度低 以及大量 即便是少量 的协方差估计误差在中累积相加 c 的定义公式 相关图 c c 1 1 N Nk ki ekr 5 8 是基于相关的PSD定义 1 直接导出的相关图谱估计算法 通过在的定义式 c 8 式中截短求和的方法可以减少这两方面的影响 这个思想就导出了BT估计子 1 1 M Mk ki BT ekrkw 9 其中 是一个偶函数 且 NM kw kwkw 10 w 同时随平滑衰减到零 由于在式 9 中对样本协方 Mkkw 0 kwk kw 差序列的滞后进行加权 故被称为滞后窗 因为由这种方法求出的功率谱是通过自相关函数间接得到的 所以称为间 接法 又称为自相关法或BT法 当较小时 式 9 的计算量不是很大 因此 M 该方法是在FFT问世之前常用的谱估计方法 当时 BT法与周期图法估计的功率谱是一样的 当 时 BTNM NM 法的偏差大于周期图法 在窗函数满足一定条件下是渐进无偏估计 方差小于 周期图的方差 分辨率比周期图法低 与窗函数选择有关 令为的DTFT W kw 10 W 1 1 M MK kiki ekekw 利用DTFT的性质 我们可以写成 ki BT ekrk 和 0 0 1 1 0 0 MwMw 的积的DTFT 0 0 1 1 0 0 NrNr 11 kwDTFTkrDTFT 记 可得DTFT 0 0 1 1 0 0 NrNr P 12 dWW PPBT 2 1 由于极大多数常用的窗函数的在处都有一个相对较窄的主峰 由 12 式 W0 可得BT谱估计子 9 式就是一个局部加权的平均周期图 6 2 3 2 BT 谱估计的性质 由于 自然也要求 选择适当的窗函数就能够使估计普 0 0 BT 拥有这样的性质 如果滞后窗式正半定的 即 则加窗的协方差序列 kw 0 W 也是正半定的 这意味着对所有的 krkw 13 0 BT 即BT谱估计具有非负性 证明证明 既然对于是实对称的 那么其 是一个实的 kw0 kDTFT W 偶函数 进一步是一个正半定序列 则对于所有 由于定义 kw 0 W 由 12 式可知 即意味着 0 P 0 W 0 BT 应该注意的是 有些滞后窗 如三角窗 并不满足 13 式中的假设 因此 这些窗函数的使用可能导致估计谱出现负值 2 3 3 BT 谱估计在 MATLAB 中的仿真实现 仿真实例仿真实例 2 2 用BT法对信号 进行功率谱估计 其中为白 213 0 2sin 2 2 02sin 20 nwnnnx nw 噪声 对信号的采样率为 1Hz 满足取样定理 n 0 N 1 程序代码见 nx s f 附录 当N 512时 MATLAB仿真得到的功率谱估计图为 00 10 20 30 40 50 60 70 80 91 5 0 5 10 15 20 25 30 35 三 三 三 Hz三 三 三 三 三 dB三 512三 三 三 三 三 三 图2 BT法得出的功率谱估计图 通过实验仿真可以直观地看出功率谱估计中BT法和周期图法所得到的结果 是一致的 其特点是离散性大 曲线粗糙 方差较大 但是分辨率较高 7 2 4 平均周期图法 周期图是信号功率谱的一个有偏估值 而且 当信号序列的长度增大 到无穷时 估值的方差不趋于零 因此 随着所取的信号序列长度的不同 所得到的周期图也不同 这种现象称为随机起伏 由于随机起伏大 使用周 期图不能得到比较稳定的估值 一些学者对此作了改进 为了减小随机起伏 M S Bartlett 特提出平均周期图法 即先把信号 序列分为若干段 对每段分别计算其周期图 然后取各个周期图的平均作为 功率谱的估值 平均周期图可以减小随机起伏 但是 如果信号序列不是足 够长 由于每段序列长度变短 功率谱估值对不同频率成分的分辨能力也随 之下降 另一种改进方法是将周期图与一个适当的频域窗函数相褶积 从而 对周期图产生平滑作用 以减小随机起伏 加窗处理的结果虽然可以使随机 起伏减小 但也会使周期图的分辨能力下降 P O Welch 提出一种把加窗处理与平均处理结合起来的方法 先把分段 的数据乘以窗函数 进行加窗处理 分别计算其周期图 然后进行平均 Welch 方法是较常用的一种计算方法 为了得到较好的功率谱估值 加窗和 平均处理均应兼顾减小随机起伏和保证有足够的谱分辨率两个方面 2 4 1 Bartlett 方法及在 MATLAB 中的仿真实现 一 Bartlett方法 Bartlett方法的基本思想很简单 为了减少周期图的起伏 把长度为的观N 测样本分割为长度为的个子样本 然后对从这些子样本得到的周MMNL 期图对于每一个取平均 将Bartlett方法表示成如下的数学形式 令 14 1tMjytyj Lj Mt 1 1 表示第个子样本的观测值 并令j 15 2 1 1 M i ti jj ety M 表示对应的周期图 Bartlett谱估计由下式给出 16 L j jB L 1 1 由于Bartlett方法是在长度为的数据段上计算的 故其分辨率应该约为 MM 1 因此 与与原来周期图方法比较 Bartlett方法的分辨率降低了倍 作为对分L 辨率降低的回报 可以预料 Bartlett方法的方差减小了 事实上 可以发现 Bartlett方法把周期图的方差减小了倍 在选择 或 时 分辨率和方差LML 之间的折中是很明显的 二 Bartlett方法的性质 8 正如我们所知道的 16 式中的可以写成 j 17 1 1 M mk ki j j ekr 其中是第个子样本的协方差序列 把 17 式代入 16 式得 jrj 18 ki M Mk L j j B ekr L 1 11 1 我们发现在形式上与使用矩形窗的BT估计很相似 子样本协方差在 B krj 上的平均是 ACS 的一个估计 但是 在 18 式中 ACS 的估计并没有有j kr 效地利用可以得到的数据滞后的乘积 特别是当接近时 ktyty k1 M 实际上 对于 只有大约可以滞后乘积的被用来形成 18 式中 1 MkM 1 ACS 的估计 可以预料 这些滞后的方差要大于BT估计中对应的滞后 kr 类似地 方差要比的方差大 此外 由于Bartlett方法使用了一个 B BT 固定的矩形滞后窗 所以其分辨率 泄漏折中的灵活性也不如BT法 由于上述 原因 我们总结为 在形式上 14 16 式定义的Bartlett估计与使用长度为的固定矩形滞M 后窗的BT估计类似 但其方差略大于后者 由此和以前讨论的BT谱估计的性质 可知 与基本的周期图方法比较 Bartlett估计降低了分辨率 减少了方差 其 倍数均为 矩形窗的主瓣要比极大多数其他滞后窗的主瓣窄 所以 MNL 在BT类谱估计中 Bartlett估计可望得到最少的谱峰兼并 即有最大的分辨率 但有最明显的泄漏 三 Bartlett估计的MATLAB仿真 仿真实例仿真实例 3 3 用Bartlett估计对信号进行 tntttx 1102sin2602sin 功率谱估计 其中为高斯白噪声 对信号的采样频率为 500Hz 满 tn tx s f 足抽样定理 N 500 源程序代码见附录 得到的功率谱估计图如下 Bartlett估计对周期图法德改进的思想是将信号分段进行估计 并将这些估 计结果进行平均而减少估计的方差 使得估计功率谱图变得平滑 如上例中将 501点的信号分为三段 分别作周期图法估计 然后加以平均 功率谱估计结果 如图3所示 将信号分段得到的估计值显然平滑了许多 9 050100150200250300350400450500 60 40 20 0 三 三 三 Hz三 三 三 三 三 dB三 三 三 三 三 三 三 三 500三 三 三 三 050100150200250300350400450500 10 20 30 40 50 三 三 三 Hz三 三 三 三 三 dB三 三 500三 三 三 三 三 三 三 三 三 三 三 三 图3 Bartlett法频率谱估计图 2 4 2 Welch 方法及在 MATLAB 中的仿真实现 一 Welch谱估计 Welch方法是Bartlett方法经过两个方面的改进得到的 首先 Welch方法 中的数据段允许重叠 其次 在计算在周期图前 对每一个数据段加窗 为了 用数学方式描述Welch方法 令 19 tKjytyj 1 Sj Mt 1 1 表示第个数据段 在 19 式中 是第个序列起始点 若 则j kj1 jMK 序列间没有重叠 但是都是邻接的 我们得到用于Bartlett方法中的样本分割 得到个数据子样本 但是 在Welch方法中的推荐值为MNLS K 在这种情况下 得到个数据段 前后数据段之间重叠50 2 MK NMS 2 计算与对应的加窗周期图 tyj 20 2 1 1 M i ti jj etyt MP 其中表示时间窗的功率 P t 21 2 1 1 M i t M P 10 通过对 20 式中的加窗周期图取平均得到PSD的Welch估计 22 S j jW S 1 1 要解释上述对Bartlett方法的改进 从而导出Welch方法的理由很简单 通 过允许数据段之间的重叠 就有更多的周期图参与 22 式中的平均 故可以降 低PSD估计的方差 在周期图计算中引入窗函数 就能够更有效地控制PSD估 计中的偏差 分辨率性质 除此之外 时间窗还可用来对位于没个子样本末端的 数据给予较小的加权 这样 尽管子样本间是重叠的 但可以使相邻的子样本 间的相关性减弱 通过 22 式中平均 这种去相关的重要作用能够更有效地降 低方差 二 Welch谱估计的MATLAB仿真 仿真实例仿真实例 4 4 用Welch谱估计对信号进 tntttx 1102sin2602sin 行功率谱估计 其中为高斯白噪声 对信号的采样频率为 500Hz tn tx s f 满足抽样定理 N 500 源程序代码见附录 得到的功率谱估计图如下 050100150200250300350400450500 60 40 20 0 三 三 三 Hz三 三 三 三 三 dB三 三 三 三 三 三 三 三 500三 三 三 三 050100150200250300350400450500 80 60 40 20 0 三 三 三 Hz三 三 三 三 三 dB三 三 500三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 图4 采用样本混叠增加分段数得到的功率谱 11 050100150200250300350400450500 80 60 40 20 0 三 三 三 Hz三 三 三 三 三 dB三 三 三 三 三 三 三 三 500三 三 三 三 050100150200250300350400450500 80 60 40 20 0 三 三 三 Hz三 三 三 三 三 dB三 三 三 welch三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 图5 Welch 平均修正周期图法得到的功率谱估计图 增加分段数可以进一步降低估计的方差 然而若每段中的数据点数太少 就会使得估计的频率分辨率下降很多 在样本信号总点数一定的条件下 我们 可以采用使分段相互重叠的方法来增加分段数 而保持每段中信号点数不变 这样就在保证频率分辨率的前提下进一步降低估计的协方差 如f仿真实例4所 示 将xn分为5段 而且每段之间有重叠 对每段分别求FFT后求平均值 得出 功率估计值 功率谱估计结果如图4所示 与图3相比 显然估计值又平滑了许 多 另外 采用加窗法也可以降低估计的协方差 这也是窗函数应用的一个方 面 即在计算周期图法之前 对数据分段并加非矩形窗 如海明窗 汉宁窗或 凯瑟窗等 然后再采用分段长度一半的混叠率 就能够大大降低估计方差 这种方法称为Welch平均修正周期图法 简称Welch 法 仿真实例4是采用海明 窗的Welch 法 谱估计结果如图5所示 与图3 图4相比 显然Welch 法得出 的功率谱估计显然平滑性较好 三 不同窗函数的Welch谱估计的MATLAB仿真 仿真实例仿真实例 5 5 用Welch方法对信号 进行功率谱估计 其中为白 213 0 2sin 2 2 02sin 20 nwnnnx nw 噪声 对信号的采样率为 1Hz 满足取样定理 n 0 N 1 程序代码见 nx s f 附录 当N 1000时 窗函数采用矩形窗 Hanning窗 hamming窗 MATLAB仿真得 到的功率谱估计图分别为 12 00 10 20 30 40 50 60 70 80 91 50 0 50 三 三 三 Hz三 三 三 三 三 dB三 三 三 三 00 10 20 30 40 50 60 70 80 91 0 20 40 三 三 三 Hz三 三 三 三 三 dB三 hanning三 00 10 20 30 40 50 60 70 80 91 0 20 40 三 三 三 Hz三 三 三 三 三 dB三 hamming三 图6 Welch法加不同的窗函数得出的功率谱估计图 仿真结果分析 由图6可以看出 使用不同的窗函数谱估计的质量是不一样 的 矩形窗的主瓣较窄 分辨率较好 但方差较大 噪声水平较高 而Haning 窗和Hamming窗的主瓣较宽 分辨率较低 但方差较小 噪声水平较低 因此 在进行谱分析时选择何种窗函数 要视具体情况而定 如果强调高分辨率 能 精确读出主瓣频率 而不关心幅度的精度 例如测量震动物体的自震频率 可 以选用主瓣宽度比较窄的矩形窗 对受到强干扰的窄带信号 若干扰靠近信号 则可选用旁瓣幅度较小的窗函数 若离开通带较远 则可选用渐近线衰减速度 比较快的窗函数 总之 要针对不同的信号和不同的处理目的来选择合适的窗 函数 这样才能得到良好的效果 2 5 小结 本章先研究周期图谱估计 从均值 方差 分辨率讨论周期图法的性能 估计的均值得出当时 是无偏估计 当 N 为有限值时 是有偏估计 N 估计的方差得出不是的一致估计 随着 N 增大 谱估计的起伏增 大 增大时 分辨率提高 但会使的起伏加剧 可见方差与分辨率是N y 一对矛盾 为此 提出了改进周期图法的几种途径 即改善窗口形状 平均和 平滑 BT 法中涉及到两次加窗 一次是用于截取数据 叫数据窗 一次是用于 截取相关函数 叫延迟窗 延迟窗常用三角窗 相关法的两次加窗将造成泄漏 效应 影响到谱峰分辨率和方差 平均就是将截取的数据段再分成 L 个小段 分别计算功率谱后取功率谱的平均 这种方法使估计的方差减少 但偏差加大 13 分辨率下降 平滑是用一个适当的窗函数与算出的功率谱进行卷积 使谱线平 滑 这种方法得出的谱估计是无偏的 方差也小 但分辨率下降 现在比较常 用的改进方法是 Welch 法 又叫加权交叠平均法 简记为 WOSA 法 这种方法以 加窗 加权 求取平滑 以分段重叠求得平均 因此集平均与平滑的优点于一 体 同时也不可避免带有两者的缺点 因此归根到底是一种折中 Welch 法得 出的是无偏谱估计 当段数 L 增大时方差减小 趋于一致估计 分段平均减小 了由数据样本本身的随机性带来的方差 段数越多方差越小 但分辨率下降 非参数谱估计法的基础是 FFT 它们出现较早 运算量较小 物理概念明 确 便于工程实现 仍是目前常用的谱估计方法 有的已被固化到 FFT 仪和谱 分析仪中 但有一些难以克服的缺点 谱的分辨率较低 它正比于 N 是 数据长度 加窗的坏影响不可避免 较宽的主瓣降低分辨力 较大的旁瓣有 可能掩盖真实谱中较弱的成分 或是产生假的峰值 而没有一个窗函数能使谱 估计在方差 偏差和分辨率各方面同时改善的 使用窗函数只是一种折中的技 巧 不是解决问题的根本办法 不是真实谱的一致估计 且 N 增大时谱曲线起 伏加剧 基于上述原因 非参数谱估计还有继续研究的必要性 3 空间谱估计的算法与实现 3 1 引言 本章将讨论使用一个元无源阵元组成的阵列对个辐射源进行定位的问mn 题 这里 从辐射源来的能量可能是声音 电磁波等等 接收传感器将该能量 转换为电信号 接收阵元可能为电磁天线 水声传感器或者地震检波器 这类 问题在雷达或声纳系统 通信 天体物理学 生物医学 地震水下监视以及其 他许多领域具有广泛的应用 这一问题的基本内容包括确定能量在空间 空气 水下或地下 的分布 这里能量高度集中的点就是辐射源 这类问题称为空间 谱估计 之所以称为谱估计还因为辐射源定位问题和前面讨论的谱估计有着非 常紧密的联系 事实上我们将看到 前面所述的方法几乎都可以用作推导辐射 源定位问题的解 3 2 全向天线的波束下倾 3 2 1 全向天线的波束下倾的基本原理 广播用的发射机天线建在天线塔上 希望覆盖广大的地域 在覆盖区的边 沿 由于距离远 地球表面弯曲 电磁波场强衰减很快 因此为了有效地覆盖 既定的区域 保证区域内场强不低于特定值 通常采用波束下倾的方法 将轴 向排列的半波振子天线 通过调节天线的轴向间距 馈电的相位 使得轴向天 线阵列的方向图实现波束下倾 轴向排列的半振子天线结构 阵元之间的距离是 垂直轴线与电磁波辐d 射方向的夹角是 相邻阵元馈入信号的相位差是 相邻阵元发出电磁波到 14 达同一地点的程差是 相邻阵元发出的电磁波到达同一地点的相位差 cosd 参看公式 23 2 cos 1 12 0 d eeeEE njjj 23 3 2 2 全向天线的波束下倾的仿真实现 根据以上讨论 编写出绘制轴向排列天线阵列的方向图程序 源程序代码 见附录 仿真结果如下图 2 4 6 8 30 210 60 240 90 270 120 300 150 330 1800 图7 轴向阵列波束下倾仿真图 3 3 天线的波达方向估计 天线阵列信号处理作为信号处理的一个重要分支 在通信 雷达 地震勘 测 射电天文等领域获得了广泛应用和迅速发展 它包括 信源定位 确定阵列到信源的仰角和方位角 甚至距离 若信源位于近 场 信源分离 确定各个信源发射的信号波形 各个信源从不同方向到达阵 15 列 即使它们在时域和频域上是叠加的 信道估计 确定信源与阵列之间的传输信道的参数 多径参数 阵列信号处理时改善蜂窝和个人通信服务系统质量 覆盖范围和容量的一 个强有力的工具 建立在波达方向估计 波束形成基础上的智能天线的应用 抑制了干扰信号 改善了欲接收信号的信噪比 降低了数字通信系统的误码率 将接收天线阵列用于公众网的反向连接 客户到基站 时 多个接收天线能够 收集更多的信号能量 若天线在空间足够分离或极化各异 则多个天线能够提 供很好的分集接收 并抑制多径传输引起的衰落 这些好处可以扩大基站的覆 盖范围 改善通信质量 3 3 1 天线的波达方向估计的基本原理 全向天线不仅利用率不高 而且对各种信号不加区别地接收 降低了通信 质量 定点无线通信采用定向天线 大幅度地改善了通信质量 面对众多移动 用户的公众通信网基站和专用移动通信网 采用天线指向即波束可变的天线 智能天线 可以使移动通信的质量得到很大的改善 为使天线的波束指向可控 甚至形状可控 采用阵列天线是合适的 在距 离通信源足够远的空间里 可以将到达的电磁波视为平面波 对于等距离直线 阵天线 由于调制在载波上的基带信号码元宽度与波速的乘积远大于天线阵的 尺寸 因此多个天线阵元上的信号的幅度可视为不变 而它们的载波的相位差 则取决于其相互位置 尺寸 波长和到达方向 无线接收的无线信号中有许多成分 其中我们关心的是信号 天线阵列S 各个阵元接收的电磁波信号因为阵元排列位置的不同带来相位差 经过特定参 数的加权控制器处理后 进一步改变了各个阵元输出信号的相位和幅度 处w 理的目标是使得阵元输出信号的和天线输出中的成分具有最大输出 用YS 信号作为基准信号 反馈控制单元的功能就是将输出信号与基准信号的SYS 差值 即误差信号 作为调节控制加权控制器参数的依据 反馈控制是 w 使减小 中的成分加大 也就是说 天线阵列接收方向图指向了信号的 YSS 方向 以下通过两个简单的例子介绍智能天线中波达方向估计和波束形式的原理 等距离直线阵智能天线 我们把天线阵元沿轴排列 从到 若有一平x01 M 面波以角 入射线与轴的夹角 和角 入射线与轴的夹角 入射到阵列 z x 上 第号阵元上产生的信号为 它与号阵元的相位差为 K K x0 cossin 2 dKAK 24 式中 和分别是入射波的波长和阵元的间距 亦称阵因子 计入阵因子 d K A 的影响 第号阵元的输出时 即 为了使天线阵的输出满足需要 K KKx A K u 在每个阵元上 用加权因子进行控制 这样第号阵元上输出的信号为 K wK 即 若到达天线阵的信号是个 则天线阵的输出是个信 KKK xAw KKu wNN 号在个阵元上的输出的叠加 将问题简化为平面的二维问题 Mxy1sin 并用解析式表达如下 16 nxnxnxnX N 21 25 1 1 11 cos12exp cos2exp 1 d Mj d j A 26 其中 为第一个信号的入射角 1 N AAAA 21 27 M WWWW 21 28 UWAXWty HHH 29 其中 H M uuuU 110 求多个信号到达的方向 波达方向 的方法有多种 下面讨论其中两种方 法 Capon法和MUSIC算法 及其仿真验证结果 3 3 2 求波达方向的方法 1 Capon 法 Capon法亦称最小方差无畸变响应 MVDR Minimum Variance Distortion Less Response 天线阵列中的阵元数决定了阵列方向图设计中的自由度数 Capon法将阵列中可控的自由度用来形成期望的波束形状 达到对有用信号进行 提升和对无用信号进行抑制的目的 并将其优化问题表达为 WRWnyE UU H ww minmin 2 30 其约束条件为 可以证明上式的解为 1 0 AW H ARA AR W UU H UU cap 1 1 代入式 29 可以得到相应到的功率为 17 ARA P UU H cap 1 1 31 2 MUSIC法 MUSIC法亦称多重信号分类 Multiple Signal Classification 当入射 信号数小于阵元数时 的矩阵有与进入阵列的信号数目相等的NMMM UU R 非零特征值及个为零的特征值 NM 是相应的噪声特征矢量 1 Nn vV 2 N v M v n V UU R 因为与的正交性 分母很小 峰值很大 这样可以得出MUSIC法的空A n V 间谱为 AA AA P H H music 32 其中 称为噪声子空间的正交投影估计 H nnV V 3 仿真 以下是对七单元线性天线阵在四信号输入情况下了仿真 仿真实例仿真实例 6 6 设信号1从方向入射 信号2从方向入射 信号3从4 3 6 方向入射 信号4从方向入射 下面是用两种方法 分别用两种坐标 直4 3 角坐标和极坐标 求波达方向估计 程序代码见附录 仿真结果如下 4999999999999999800000000000000 9999999999999999600000000000000 14999999999999999000000000000000 30 210 60 240 90 270 120 300 150 330 1800 图8 Capon法作出的波达方向估计 极坐标 18 020406080100120140160180 10 0 10 5 10 10 10 15 10 20 10 25 10 30 10 35 图9 MUSIC法作出的波达方向估计 直角坐标 500000 1000000 1500000 30 210 60 240 90 270 120 300 150 330 1800 图10 MUSIC法作出的波达方向估计 极坐标 19 020406080100120140160180 10 15 10 10 10 5 10 0 10 5 10 10 图11 Capon法作出的波达方向估计 直角坐标 从对以上四图的分析可以得到如下结论 两种方法估计得都比较准确 在 处有尖锐的方向图线 MUSIC法方向图像幅度更大 4 3 3 4 6 3 4 天线阵的波束形成 1 原理 我们以等距离圆阵为例来讨论天线阵的波束形成 我们取圆阵的圆心为 在阵列圆上任取一点 把天线阵元顺序定为OB 到 若有一平面波以角入射到阵列上 第号阵元上产生的信号为01 M K 它与到达阵元中心的波前的相位差是 K x kK r jAcos2exp 与 分别是入射波的波长与阵列圆的半径 亦称阵因子 为了使天线阵的 r K A 输出满足需要 在每个阵元上加上加权因子控制 这样第号阵元上输出的 K K 信号为 若到达天线阵的信号是个 天线阵输出的是个信号的 KKK xA NN 个阵元上的输出的叠加 用解析表达式如下 M cos 2exp 1 kK r jA 33 H MKKK AAAA 21 H N nxnxnxnX 110 20 34 式中 是阵元以为基准顺时针画出的角度 K KOB N AAAA 21 35 H M wwwW 110 36 UWAXWty HH 37 为了求得多个信号到达的方向 波达方向 可以采用上述的 Capon MUSIC两种方法 波束形成可以采用下面的方法 当有多 个信号输入时 其中1个信号时我们关心的 个信号是N1 N 需要抑制的 方程组 38 描述了上述需求的约束条件 四个信号输入 第一个信号是 我们关心的 其余信号是需要抑制的 T HU W0 0 0 1 38 根据信号波达方向的知识及约束条件求解方程组 38 可以 AXUU 得到 110 M wwwW 39 代入式 37 可得到阵列输出的方向特性 2 仿真 下面是四个输入信号八单元圆阵列天线阵的波达方向估计和波束形成 仿真实例仿真实例 7 设信号1从方向入射 信号2从方向入射 信号3从3 2 3 方向入射 信号4从方向入射 下面是用两种算法求波达方向估计2 3 6 8 5 及波束形成 源程序代码见附录 仿真结果如下 21 050100150200250300350400 10 0 10 5 10 10 10 15 10 20 10 25 10 30 10 35 图12 Capon算法波达方向估计 050100150200250300350400 10 15 10 10 10 5 10 0 10 5 10 10 图13 MUSIC算法波达方向估计 22 0 5 1 1 5 30 210 60 240 90 270 120 300 150 330 1800 图14 波束形成 3 5 小结 本章在对非参数谱估计的理解基础上 讨论了其在空间领域进行推广 具 体就是了采用Capon算法和MUSIC算法求天线阵的波达方向 讨论了天线阵的 波束形成 及其MATLAB的仿真实现 4 结 论 本设计的主要任务是研究非参数谱估计的算法与实现 另外对功率谱估计 算法在空间领域的推广应用尤其是阵列侧向问题进行了讨论 在经典谱估计中 无论是周期图法还是其改进的方法 都存在着频率分辨 率低 方差性能不好的问题 原因是谱估计时需要对数据加窗截断 用有限个 数据或其自相关函数来估计无限个数据的功率谱 这其实是假定了窗以外的数 据或自相关函数全为零 这种假定是不符合实际的 正是由于这些不符合实际 的假设造成了经典谱估计分辨率较差 另外 经典谱估计的功率谱定义中既无 求均值运算又无求极限运算 因而使得谱估计的方差性能较差 当数据很短时 这个问题更为突出 如何选取最佳窗函数 提高频谱分辨率 如何在短数据情 况下提高信号谱估计质量 还需要进一步研究 Capon法和MUSIC算法可以用于空间谱估计测向 充分利用了天线各个阵元从空间电 磁场接收到的全部信息 而传统测向方法仅仅利用了其中一部分信息 因而具有多方面的 优越性 1 可以实现对同一一信道中存在的多个信号同时测向 2 超分辨率测向 即对处 23 于天线阵固有波束宽度以内的两个以上方向的来波同时测向 3 阵元位置可以随意放置 各个阵元的方向特性也可以是任意的 4 具有优良的测向灵敏度和准确度 5 可以较好地 克服测向