AR双谱估计的Matlab实现.doc

Journal of Shaoguan University · Natural Science 第 33 卷 第 10 期 Vol.33 No.10 AR 双谱估计的 Matlab 实现 何 莹 （韶关学院 物理与机电工程学院， 广东 韶关 512005） 摘 要 ：根据 AR 双谱估计的基本原理 ，对其 Matlab 实现方法进 行 研 究 ，给 出 了 AR 模 型双 谱 估 计的 Matlab 实 现 ，最 后 给出一个仿真测试实例，仿真结果表明了该方法的有效性. 关 键 词 ：双谱；AR 模型；Matlab 中图 分 类 号 ：TN911 文献 标 识 码 ：A 文章 编 号 ：1007-5348（2012）10-0040-04 双谱估计在通信 、雷 达 、机 械 学 、生 物 医 学 、电 磁 学 、地 球 物 理 学 、流 体 动 力 学 、故 障诊断等领域具有广 泛的应用［1］.双谱提供了比功率谱更好的特性［2-3］，它可以检验信号偏离高斯分布的程度 ，能保留信号的相位 信息，并能自动抑制加性高斯噪声的影响等 ，是分析非高斯信号、非因果、非最小相位系统的主要工具. Matlab 是一种以矩阵运算为基础的交互式的可视化软件 系 统 环 境 ，集 数 值 分 析 、矩 阵 运 算 、信 号 处 理 和图形显示于一体 ，为科学研究 、工程设计及众多学科领域提供了一种简捷 、高效的编程工具.用 MATLAB 对信号进行双谱估计和分析，其实现方法非常方便、简单. 1 双谱基本理论 设{x(n)}为零均值、实平稳随机信号，其三阶累积量 c3x(m1,m2)定义为： c3x(m1,m2)=E［x(n)x(n+m1)x(n+m2)］ 其双谱定义为三阶累积量的二维傅里叶变换 ，即有： （1） ∞ ∞ Bx(ω1, ω2)= Σ Σ c3x(m1,m2)exp［-j(ω1m1+ω2m2］ m1 =-∞ m2=-∞ 其中，|ω1|燮π, |ω2|燮π, |ω1+ω2|燮π［1］. 在实际应用中由于可供处理的观测信号一般都是具有有限长度的 ，因 此 无 法 精 确 应 用 式(2)求 解 随 机 过程的双谱， 而只能对双谱进行估计.AR 模型法在观测数据较短的情况下就可获得高分辨率的双谱估计 ， 是双谱估计的常用方法之一. 设零均值采样数据{x(n)}（n=1,2,…,M）由以下 p 阶 AR 模型表示： p （2） x(n)+Σakx(n-k)=w(n) （3） k=1 式中 w(n)为零均值、具有有限的非零累积量 γ3e 的高阶白噪声.对（3）式求三阶自相关，可得如下的正则方程［4］： p c3x(m,n)+Σakc3x(m-k,n)=γ3eδ(m,n) k=1 收 稿 日 期 ：2012-09-03 作 者 简 介 ：何莹(1974-)，女，江西崇义人，韶关???院物理与机电工程学院讲师，硕士，主要从事随机信号分析与处理方面的研究. 第 10 期 何莹：AR 双谱估计的 Matlab 实现 ·41· m0, 坌n 或 (m,n)=(0,0) （4） 上式可用矩阵表示为： Ra=b a=［1,a1,…,ap］T （5） b=［γ3e, 0, …, 0］T 坌 坌 坌 坌 坌 坌c3x (0,0) c3x (-1, 0) c3x (0, 0) … c3x (p-1, p) … … … … c3x (-p,0) 坌 坌 坌 坌 坌c3x (1,0) c3x (-p+1,0) 坌 坌 坌 坌 坌 坌 坌 坌 坌 坌 坌 坌 坌 坌 坌 坌 … … c3x (0,p) 坌 坌c3x (p,p) 坌 坌 坌 坌 坌 矩 阵 R 是 把(4)式的各个方程按 m(m=0,…,p)从 小 到 大 、而 在 m 相 同 时 n(n=0,…,m)由小到大排列而形成的. 由于 R 不是方阵，因此方程组（4）只能在最小二乘意义下用伪逆运算求解. 其中三阶累积量 c3x(m1,m2)可由下式估计： k2 1 c3x (m1,m2)= M Σx(t)x(t+m1)x(t+m2) （6） t=k1 式中：k1=max(1, 1-m1, 1-m2)，k2=min(M, M-m1, M-m2).对应（3）式的{x(n)}的双谱为： B(ω1, ω2)=γ3eH(ω1)H(ω2)H*(ω1+ω2). 2 AR 双谱的 Matlab 实现 （7） 在 Matlab 环境中要进行双谱估计必须根据算法编写相应的函数文件.在实际应用中 ，AR 模型的阶次 p 一般是未知的 ， 阶次的选择是必须要解决的问 题.AR 模型阶次确定主要有奇异值分解法和双谱互相关法 等，具体可参考文献［1］，在此不作讨论. 要获得序列{x(n)}的 AR 双谱估计，可按以下步骤实