《武汉工程大学学报》  2011年01期 71-74   出版日期：2011-01-30   ISSN:1674-2869   CN:42-1779/TQ

---

基于状态空间实现的EDS频率与时延联合估计


0引言指数衰减正弦信号（EDS: Exponentially Damped Sinusoid）是核磁共振成像、线性系统辨识、语音分析、暂态分析和生物信号分析等领域广泛应用的参数化模型［1］之一.在加性高斯白噪声条件下,指数衰减正弦信号的参数估计问题已经得到了较深入的研究.最大似然方法（ML）［2］、线性预测（如KT）［3］、子空间旋转不变［4］和子空间正交分解等一些有效的现代谱估计方法被应用到指数衰减正弦模型估计并取得了很好的性能.最近,文献［5］提出了一种非迭代前后向线性预测方法估计EDS信号模型的非线性参数,在特定的数据长度情况下,性能良好.文献［6］利用离散傅里叶变换(DFT)的峰值获得EDS信号模型的衰减因子和频率的估计值,此方法中的汉宁窗有效的减小了在衰减因子较大且非相干时的估计误差,但在衰减因子较小和相干的情况下,估计性能较差.文献［7］提出了一种迭代自适应方法(Iterative Adaptive Approach 简称IAA),这种方法不需要任何参数模型,也不需要采样是均匀这一假设,但此方法的计算复杂度较大.文献［8］提出了一种基于EDS信号模型的IMP自适应算法,此方法扩展到信号参数随时间变化的信号模型中.最近,基于双基元接收的指数衰减正弦信号的频率、衰减因子和时延［9］联合估计问题也引起了学者们的关注.文献［10］研究了一种在一般DDS(DDS: Damped and Delayed Sinusoidal)模型中扩展的数据模型—偏衰减时延正弦信号模型(PDDS),即考虑了多个指数衰信号具有同一个时延下的情形.文献［11］给出了一种基于DDS模型的子带处理方法,在文中有两种模型参数估计算法—基于子空间算法(DDSB)和基于傅里叶变换算法(DDSD).文献［11］给出了一种基于正弦信号模型的状态空间实现的频率和时延的联合参数估计方法,但未考虑指数衰减因子.根据文献［12］、［13］提出的正弦信号模型加以推广,提出一种基于状态空间实现的指数衰减正弦信号的频率、衰减因子和时延的联合估计方法,由状态过度矩阵的特征值可以获得其频率和衰减因子估计,而由观测矩阵及估计的频率及衰减因子可以获得时延估计,结果具有闭式解,不需要搜索计算,其方法具有对模型误差不敏感的优点,因而具有更好的数值稳定性.仿真结果表明在低信噪比情况下具有较高的估计精度.1信号模型及基本理论基于双基元的接收信号数据模型可以表示为：r1(n)=s(n)+q1(n)（1）
r2(n)=s(n-D)+q2(n),n=0,1,…,N-1（2）这里r1(n),r2(n)是两个接收的信号数据,n为采样的时间,s(n)=∑Pi=1αie(-di+jwi)n,信号幅度αi是未知的复常数,假定归一化的频率wi和di各不相同；q1(n)和q2(n)是任意的加性复高斯噪声.且假定正弦波数目P是先验已知或者已由某一检测方法获得P.参数D表示由两个阵元引起的相对时延,N为采样数,我们的目的是由给定的接收数据去估计时延D、各频率分量wi(i=1,2,…,P)以及相应的各衰减因子分量di(i=1,2,…,P).首先用接收的信号数据形成一个时间状态空间模型,即形成下面的两组向量,
X1(k)=［r1(k),r1(k+1),…,r1(k+M-1)］T（3）
X2(k)=［r2(k),r2(k+1),…,r2(k+M-1)］T（4）这里k=0,1,…,K-1=N-M+1,T表示转置算子,参数M表示每个向量的长度(P+1≤M≤N-P+1).将(1)和(2)代入(3)及(4)有：X1(k)=A(w)s(k)+Q1(k)（5）X2(k)=A(w)Δ(w,d,D)s(k)+Q2(k)（6）
这里
A(w)=［a1,a2,…,aP］
am=［1,e-dm+jwm,…,e(-dm+jwm)(M-1)］T
m=1,2,…,P
s(k)=
［α1e(-d1+jw1)k,α2e(-d2+jw2)k,…,αPe(-dP+jwP)k］T
Q1(k)=［q1(k),q1(k+1),…,q(k+M-1)］T
Q2(k)=［q2(k),q2(k+1),…,q(k+M-1)］T
Δ(w,d,D)=
diag{e-D(-d1+jw1,e-D(-d2+jw2),…,e-D(-dP+jwP)}第1期吴云韬,等：基于状态空间实现DDS频率与时延联合估算
武汉工程大学学报第33卷
从式（6）可以看出,频率和衰减因子信息包含在Vandermonde矩阵A(w)中,时延信息包含在对角矩阵Δ（w,D）中.设y(k)=X1(k)
X2(k)=
A(w)
A(w)Δ(w,d,D)s(k)+e(k)s(k)表示状态向量,因此可以形成如下的状态空间模型s(k+1)=Φs(k)（7）y(k)=B(w,D)s(k)+e(k)（8）这里Φ=diag［e-d1+jw1,…,e-dP+jwP］（9）称为状态过渡矩阵.B=A(w)
A(w)Δ(w,d,D)（10）称为观测矩阵.将由给定的数据去估计出系统矩阵B和Φ,如频率和衰减因子已知,那么可以利用Δ(w,d,D)的结构给出时延D的估计.2提出的算法2.1状态空间算法的实现及信号参数估计设L≥2是一给定的整数,定义一2ML×1向量
yL(k)=［yT(k),yT(k+1),…,yT(k+L-1)］T通过迭代状态空间模型有
yL(k)=ΩLs(k)+eL(k)k=1,2,…,K-L+1上式中eL(k)与yL(k)类似定义,ΩL=［BT(BΦ)T…(BΦL-1)T］T（11）这里2ML×P矩阵ΩL称为状态空间模型的可观测矩阵(通常定义L=P).定义yL(k)的采样协方差矩阵yL(L)=1K∑K-L+1k=L+1yL(k)yHL(k-L)设yL=∑2MLl=1λlulvHl为yL的奇异值分解,这里λl表示奇异值(按降序排列),ul和vl分别表示左、右奇异向量.因此信号子空间可以由矩阵=［u1,…,uP］的列矢量得到.则状态矩阵ΦT和观测矩阵BT能够用下面的形式估计得到T=11（12）T=#1L-12L（13）这里kl表示从k到l的块行(每一个块矩阵具有维数2M×P),(·)#表示矩阵的伪逆.设T的特征分解为T=EΛE-1,这里Λ=diag［λ1,…,λP］,从而频率可以由λm的相角得到(∠表示取相角),衰减因子可以由λm的实部取负对数得到,即m=∠λm（14）m=-ln［Re(λm)］,m=1,…,P（15）再根据B(w,D)的定义,Δ(w,d,D)的估计为(w,d,D)=#1MM+12M,从而时延D的估计由下式给出=∑Pm=1∠mm(,,D)-∑Pm=1(-m+m)（16）
式（16）中∠表示矩阵(,,D)的第m个对角元素的相角.2.2DDS模型的CRLB下界在基于双基元DDS模型中,设s=∑Pi=1Aie-din+j(win+pi)
sD=∑Pi=1Aie-di(n-D)+j［wi(n-D)+pi］
xi=e-din+j(win+pi)
xiD=e-di(n-D)+j［wi(n-D)+pi］
d=e-(d1+d2+…+dP)n+j［(w1-w2)-…-wP)n+(p1-p2-…-pP)］
dD=e-(d1+d2+…+dP)(n-D)+j［(w1-w2-…-wP)(n-D)+(p1-p2-…-pP)］这里Ai(i=1,2,…,P)是实振幅,pi(i=1,…,P)是相角.先给出此模型的对数似然函数是
lnf(r1,r2;A1,…,AP,p1,…,pP,
w1,…,wP,d1,…,dP,D)=
-Nln(2π)-Nlnσ2-
1σ2∑N-1n=0(r1-s)2+∑N-1n=0(r2-sD)2（17）这里σ2是噪声功率.然后写出Fisher信息矩阵,对似然函数中的每个参数求导可以得到Fisher信息矩阵中每个元素的表达式,即Iij=2lnfGiGj(i,j=1,…,4P+1)（18）其中,G=［A1,…,AP,p1,…,pP,w1,…,wP,d1,…,dP,D］,Gi表示向量G中的第i个元素,Gj表示向量G中的第j个元素.则频率、衰减因子和时延的CRLB下界可以由下式得到,令J=I-1(A1,…,AP,p1,…,pP,w1,…,
wP,d1,…,dP,D)（19）每个参数对应的CRLB界限就是上式中对角线上的元素数值.CRLB(wi)=J(2P+i,2P+i)=
1σ2∑N-1n=0［Air1n2(xi+x*i)-
A1…APn2(d+d*)+Air2(n-D)2(xiD+
x*iD)-A1…AP(n-D)2(d+d*)］CRLB(di)=J(3P+i,3P+i)=
1σ2∑N-1n=0［-Air1n2(xi+x*i)+
A1…APn2(d+d*)+4A2in2e-2din-
Air1(n-D)2(xiD+x*iD)+A1…AP(n-D)2
(dD+d*D)+4A2i(n-D)2e-2di(n-D)］CRLB(D)=J(4P+1,4P+1)=
1σ2∑N-1n=0［r2Ai(xiD+x*iD)(w2i-
d2i)+2jr2Aidiwi(xiD-x*iD)+4A2id2ie-2di(n-D)+
A1…AP(∑Pi=1d2i)(dD+d*D)+
2jA1…AP(∑Pi=1di)(w1-…-wP)(d*D-dD)-
A1…AP(w1-…-wP)2(d*D+dD)］3计算机仿真实验及性能比较在计算机仿真中,信号s(n)=α1e(-d1+jw1)n+α2e(-d1+jw2)n,α1=α2=1/2,w1=0.3π,w2=0.6π,d1=0.001,d2=0.002 和时延D=0.7 s.采样时间为1 s,独立实验次数为500.采样点数N=200,取L=2,M=25.假定加性噪声为白噪声,qi(t)(i=1,2)互不相关.图1分别给出了两个分量的频率估计的均方误差随信噪比的对比图,图2分别给出了两个分量的衰减因子估计的均方误差随信噪比的对比图,图3给出了时延估计的均方误差随信噪比的对比图,同时基于两阵元数据的ESPRIT算法的频率估计结果,相应的CRLB下界也一同给出.从图1、图2及图3中可以看出,当信噪比大于-5dB时,上述方法参数估计的均方误差接近于对应的CRLB下界,同时,该方法和ESPRIT方法具有相同的频率及衰减因子估计性能.图1频率估计的均方误差随信噪比的对比图
Fig.1Mean square error of frequency estimation versus SNR
注：σw为频率均方误差.图2衰减因子估计的均方误差随信噪比的对比图
Fig.2Mean square error of damped factor estimation versus SNR图3时延估计的均方误差随信噪比的对比图
Fig.3Mean square error of Delay estimation versus SNR
注：σd为衰减因子均方误差；σD为时延均方误差.4结语计算机仿真实验证明了本文的方法性能在低信噪比和较小的衰减因子情况下指数衰减时延正弦信号模型的频率，衰减因子和时延估计都能达到相应的CRLB下界.但在衰减因子较大的情况下，性能还有待提高.