中图分类号:TN911.23
文献标志码:A
文章编号:1671-5497(2016)03-0923-06
Standard normalized correlation-based MUSIC method for harmonic signal frequency estimation in α Noise
0 引 言噪声背景下的谐波恢复是信号处理领域的核心问题之一, 传统的谐波恢复方法一般都是假设背景噪声为高斯噪声。然而在自然界中, 能够真正符合高斯分布的随机过程毕竟是少数。在大量的信号处理实践中, 其背景噪声通常为非高斯噪声[1, 2, 3, 4], 例如由放电引起的大气噪声, 水下生物引起的水下声学噪声, 数字电路的开关噪声, 无线通信中的信号窜扰等噪声。它们比高斯噪声出现大幅度的数据突变的频率更高。反映在时域上, 它们显现出大量的显著尖峰脉冲特性; 反映在概率密度上, 它们具有比高斯分布更加厚重的拖尾现象[5, 6]。
研究结果[6]表明, 这种非高斯冲击噪声可以用 α稳定分布来表征, 一般称为 α稳定分布噪声(简称为 α噪声)。 α稳定分布是唯一一种满足广义中心极限定理的广义高斯分布, 高斯分布是 α稳定分布中, 当 α=2的一个特例, 因此 α稳定分布比高斯分布具有更广泛的代表意义。 α稳定分布的一个重要性质是它不存在 p>α(当 α<2时)的统计量[6], 当 α<1时, 甚至连均值的运算也不可能。因此, 在 α噪声存在的场合, 现有的各种基于二阶矩、高阶矩和高阶累积量的信号处理方法均失效[7]。目前, 任何涉及 α稳定分布的信号基本上都是采用 p<α<2的分数低阶矩[7]或由其派生的统计工具进行处理。但采用分数低阶矩类方法解决 α噪声背景下的谐波恢复问题时, 由于 α稳定分布过程不存在有限的方差, 因此, 传统信号处理常用的方法, 如极大似然法、最小二乘法等, 均无法使用[6]; 分数低阶矩类方法由于其本身固有的非线性, 无法得到该问题的闭式解[8]。针对上述问题, 文献[8]将 α噪声背景下的谐波信号用一个 α白噪声驱动的高阶AR模型表示, 然后利用分数低阶矩中, 当且仅当其阶次 p=1时具有线性的性质, 采用GCC(Generalized covariation coefficient)方法巧妙地解决了这一问题。但该方法仅适于阶次 p=1的情况, 当 p≠1时, 其谐波信号的频率估计精度显著下降, 因而该方法仅适用于 α>1的稳定分布噪声的情况。
针对分数低阶矩类方法所存在的非线性问题, 文献[9, 10]指出, 若 α稳定分布过程 y(n)可以用一个极点均在单位圆内的AR模型描述, 则适当定义的归一化相关函数估计(Normalized correlation estimate)存在, 因而可以用于 α稳定分布信号处理。然而这种归一化相关函数虽然具有线性特性, 但其AR模型极点均在单位圆内这一约束条件却使得它无法用来解决谐波信号频率估计问题。本文对标准归一化相关函数的定义域进行了扩展, 证明了对于 α稳定分布过程( 0<α≤2), 其标准归一化相关函数存在。在此基础上提出了一种 α噪声背景下, 基于标准归一化相关函数的谐波信号频率估计的MUSIC方法。
1 含噪谐波信号的数学模型本文主要考虑这样一类噪声中的谐波信号频率估计问题:
x(n)=∑i=1qαiexp(j(ωin+φi))(1)
式中: x(n)为待测无噪谐波信号; αi、ωi(i=1,…,q)分别为谐波信号幅值和频率, 均为实常数; φi为在 [-π,π]内均匀分布的独立同分布随机相位; q为谐波个数, 可以是未知的; j为虚部标志。
x(n)的含噪观测信号 y(n)为:
y(n)=x(n)+e(n),n=1,…,N(2)
式中: e(n)为独立同分布Sα S稳定分布随机噪声, 且与 x(n)相互独立; N为观测序列长度。
式(2)中的 SαS稳定分布过程的定义为[6]:
若存在参数 0<α≤2,γ≥0,-1≤β≤1和实数 a使得随机变量 x的特征函数具有式(3)的形式:
φ(u)=exp{jau-γ|u|α[1+jβsgn(u)ω(u,α)]}(3)
式中:
ω(u,α)={tan(πα/2),α≠1(2/π)log|u|,α=1(4)
sgn(u)={1,u>00,u=0-1,u<0(5)
则随机变量 x服从 α稳定分布, 记为: x~Sα(γ,β,a)。
式(3)中, 若有 β=a=0,则称随机变量 x服从Sα S稳定分布。
显然, 谐波信号 x(n)的 q个谐波频率可由特征方程(6)的根( |zi|=1)决定[11], 方程(6)的表达式如下:
∏i=1p(z-zi)(z-z*i)=∑i=02pbiz2p-i=0(6)
因而式(1)所示的谐波信号所对应的差分方程为:
x(n)+∑i=12qbix(n-i)=0(7)
将式(2)代入式(7)可得:
y(n)+∑i=12qbiy(n-i)=e(n)+∑i=12qbie(n-i)(8)
由式(8)可见: α白噪声背景下的谐波过程 y(n)为一个特殊的ARMA( 2q,2q)过程, 该ARMA模型的极点均在单位圆上, 且其AR参数与MA参数相同。
2 本文方法式(8)中, 由于 e(n)为 α稳定分布随机过程, 所以 y(n)不存在有限的方差。若采用分数低阶矩类方法, 由于该方法本身所固有的非线性( 0<p<α,p≠1), 所以无法得到含噪谐波信号的封闭表达式[8]。
文献[9, 10]指出, 若 α稳定分布过程 y(n)可以用一个极点均在单位圆内的AR模型描述, 则有如下归一化相关函数估计存在:
^ryy(m)=∑n=1Ny(n)y(n+m)∑n=1Ny2(n)(9)
但该归一化相关函数的定义要求 α稳定分布过程 y(n)所满足的AR模型的极点均在单位圆内, 这一条件过于严格。由于谐波信号模型的零极点均在单位圆上, 因而无法使用这种归一化相关函数来解决谐波信号的频率估计问题。
为了解决这一问题, 本文将文献[9, 10]所提出的归一化相关函数的定义域进行扩展, 提出如下对于 α稳定分布过程的标准归一化相关函数存在的定理。
定理 设 y(n)为服从 Sα(σ,β,a)稳定分布的随机过程, 且 0<α≤2,则有如下标准归一化相关函数存在:
ryy(m)=Ryy(m)Ryy(0)=E[y(n)y(n+m)]E[y2(n)](10)
证明 由标准归一化相关函数的定义, 有:
Ryy(m)=E[y(k)y(k+m)]E[y2(k)]=limN→∞∑k=1N[y(k)y(k+m)]limN→∞∑k=1N[y2(k)](11)
因为:
∑k=1N[y(k)-y(k+m)]2≥0(12)
将式(12)展开, 有:
∑k=1N[y(k)-y(k+m)]2=∑k=1Ny2(k)-2∑k=1Ny(k)y(k+m)+∑k=1Ny2(k+m)(13)
由于 y(n)为服从 Sα(σ,β,a)稳定分布, 因此y(k)具有平稳遍历性, 所以当N很大时, 有:
∑k=1Ny2(k)=∑k=1Ny2(k+m)(14)
则式(12)为:
∑k=1N[y(k)-y(k+m)]2=2∑k=1Ny2(k)-2∑k=1Ny(k)y(k+m)≥0(15)
所以有:
∑k=1Ny2(k)≥∑k=1Ny(k)y(k+m)(16)
将式(16)代入式(10), 有:
ryy(m)=E[y(k)y(k+m)]E[y2(k)]=limN→∞∑k=1N[y(k)y(k+m)]limN→∞∑k=1N[y2(k)]≤limN→∞∑k=1N[y2(k)]∑k=1N[y2(k)]=1(17)
这说明式(10)所示的标准归一化相关函数有界, 即对于 y(n)~Sα(γ,β,a)的稳定随机过程, 式(10)所示的标准归一化相关函数存在。
证毕。
定理将文献[9, 10]所提出的归一化相关函数的定义域扩展至任何服从 Sα(σ,β,a)稳定分布的随机过程。注意到在定理中, α稳定分布过程的特征指数 α的值域为 0<α≤2,当 α=2时, α稳定分布退化为高斯分布, 因此, 标准归一化相关函数不仅对 0<α<2稳定分布信号有效, 对高斯分布信号也是适用的。
特别值得指出的是, 与目前在 α稳定分布信号处理中广泛使用的分数低阶矩类方法相比, 本文所提出的标准归一化相关函数的最大优点在于它是一个线性函数, 这不仅大大简化了 α稳定分布信号处理的分析和计算过程, 而且使得现有的各种在高斯背景噪声下基于二阶矩的信号处理方法均可以很容易地应用于 α稳定分布信号处理, 这一点对于 α稳定分布噪声背景下的信号处理具有重要意义。
在实际应用中, 标准归一化相关函数可以使用式(18)估计:
^ryy(m)=∑k=1N[y(k)y(k+m)]∑k=1N[y2(k)](18)
根据定理, 对式(2)所示的含噪谐波信号取标准归一化相关函数, 得:
ryy(m)=∑i=1qejωim+σ2δ(m)(19)
式中: σ2为Sα S稳定分布白噪声方差。
将式(19)代入如下 k×k(k>q)维相关函数矩阵, 有:
R=[1ryy(1)…ryy(k-1)ryy(1)1…ryy(k-2)︙︙︙ryy(k-1)ryy(k-2)…1](20)
可得:
R=∑i=1qFiFHi+σ2I=FFH+σ2I(21)
式中: Fi=[1ejωi…ej(k-1)ωi]T为一个 k×1维复列向量; F=[F1F2…Fq]为 k×q维复矩阵。
由于矩阵 F完全是由信号的 q个谐波构成, 因此称矩阵 F为谐波矢量阵。
不难证明, 式(20)中的矩阵 R的秩为 q。对矩阵 R进行特征值分解, 并将特征值递减排列, 得:
R=UΣUH(22)
式中: U为其特征矢量阵; Σ为其特征值矩阵, Σ=diag[σ21,σ22,…,σ2k]。
由于矩阵 R的秩为 q,且其按特征值递减排列, 所以有:
σ21≥σ22≥…≥σ2q>σ2=…=σ2(23)
即矩阵 R有 q个较大特征值, 而另外 k-q个特征值明显较小, 因此, 当谐波个数未知时, 可以据此作为谐波个数的估计。
将特征矢量矩阵 U分块, U=[U1 ︙ U2]T, 其中 U1由 U的前 q个较大特征值所对应的特征矢量构成, U2由 U的后 k-q个较小特征值所对应的特征矢量构成, 并依此对应地将特征值矩阵 Σ分块成 Σ1和 Σ2,则式(22)可写成:
R=[U1U2][Σ1Σ2][U1U2]H(24)
用 U右乘式(24)两端, 并进行整理, 有:
RU2=σ2U2(25)
另一方面, 用 U2右乘式(21)两端, 得:
RU2=FFHU2+σ2U2(26)
比较式(25)和式(26), 明显有:
FFHU2=0(27)
由于谐波矢量阵 F为高矩阵, 因此可得:
FHU2=0(28)
这说明 q个谐波信号分量均为由式(28)所构成的特征方程的根。因此令:
P(jω)=1|FU2|2=1FU2UH2FH(29)
然后通过搜索 P(jω)的峰值即可得到 q个谐波信号频率的估值。这就是本文所提出的基于标准归一化相关函数的谐波信号频率估计的MUSIC方法。
3 仿真实验取含噪观测信号模型为:
y(k)=s(k)+n(k)=cos(2πf1k)+cos(2πf2k)+n(k)(30)
式中: f1、f2均为归一化数字频率,f1=0.26rad/s,f2=0.27rad/s;n(k)为附加噪声。
数据长度 N=512,归一化相关函数矩阵 R的维数为50× 50。
仿真实验1 为了考察本文方法在 α噪声背景下的谐波信号频率估计性能, 取 n(k)~Sα(γ,0, 0)的独立同分布的稳定分布噪声, α=1.1,n(k)按文献[12]方法产生; 其信噪比(SNR)按文献[13]所提出的广义信噪比(GSNR)计算, 即:
GSNR=10log(1γN∑k=1N|s(k)|2)(31)
取GSNR=0 dB, 使用本文方法所得到的一个典型实验结果如图1所示。
由图1可见, 在 α噪声背景下, 本文所提方法的谐波信号频率估计的精度及分辨率比较高。
仿真实验2 为了考察本文所提方法在各种不同GSNR条件下的谐波信号频率估计性能, 采用与实验1相同的实验条件及实验方法, 但GSNR从0 dB到25 dB变化, 变化间隔为5 dB, 每个实验点均进行50次Monte Carlo实验, 所得到的谐波信号频率估计的均方根误差(Root-mean-square error, RMSE)随GSNR变化曲线如图2(a)所示。
由图2(a)可以看出, 本文所提方法对 α噪声具有很强的抑制能力, 并且随着GSNR的增加, 其谐波信号频率估计的RMSE逐渐降低。
仿真实验3 为了考察本文所提方法在各种不同特征指数的 α背景噪声条件下的谐波信号频率估计性能, 采用与实验1相同的实验条件及实验方法, 但 α噪声的特征指数从0.5到2变化, 变化间隔为0.3, 每个实验点均进行50次Monte Carlo实验, 所得到的谐波信号频率估计的均方根误差(RMSE)随特征指数 α的变化曲线如图2(b)所示。
由图2(b)可见, 谐波信号频率估计的RMSE随 α噪声的特征指数 α的增加变化不大, 这说明标准归一化相关函数对不同特征指数的 α噪声具有很强的适应能力。特别是在该实验中, 当特征指数 α=2时, α稳定分布噪声已经退化为常见的高斯噪声, 因而对高斯噪声背景下的谐波信号处理依然有效, 说明本文方法具有良好的鲁棒性。
4 结 论(1)对于特征指数 0<α≤2的任何稳定分布过程, 标准归一化相关函数存在, 从而为 α稳定分布信号处理提供了一个更加有效、方便的信号处理工具。
(2)标准归一化相关函数是一个线性统计算子, 这不仅大大简化了 α稳定分布过程信号处理的分析和计算过程, 而且使得目前许多在高斯噪声背景下基于二阶矩的优化算法, 如最小二乘法、极大似然法等, 都可以应用到 α噪声背景下的信号处理中。
(3)本文方法不仅具有优良的谐波频率估计性能, 而且对于包括高斯噪声在内的 α稳定分布噪声均具有很强的鲁棒性。
The authors have declared that no competing interests exist.