系统估计方法及程序和记录媒体、系统估计装置的制作方法

专利名称：系统估计方法及程序和记录媒体、系统估计装置的制作方法
技术领域：
本发明涉及系统估计方法及程序和记录媒体、系统估计装置，特别是使用根据H∞评价基准开发的超高速H∞滤波器的高速H∞滤波器算法，同时实现状态估计的稳健化和忘却系数的最佳化的系统估计方法及程序和记录媒体、系统估计装置。
背景技术：
通常，所谓系统估计是根据输入输出数据估计系统的输入输出关系的数学模型(传递函数、或脉冲响应等)的参数。作为代表性的应用例，有国际通信的回波消除器、数据通信的自动均衡器、音响系统的回声消除器和声场重放以及汽车等的主动噪声控制等。详细情况参照非专利文献1.1993年电子情报通信学会的「数字信号处理手册」等。
基本原理图8表示系统估计用的结构图的例子(未知系统也可以用IIR(Infinite ImpulseResponse；无限脉冲响应)滤波器表达)这种系统具备未知系统1、自适应滤波器2。又，自适应滤波器2具有FIR数字式滤波器3、自适应算法4。
下面对同定未知系统1的输出误差方式的一个例子进行说明。其中，uk是未知系统1的输入，dk是作为所希望的信号的系统的输出，d^k是滤波器的输出。(还有，「^」是估计值的意思，是附在字符正上方的符号，但是由于输入不方便，只是记载于字符的右上方。下同。)作为未知系统的参数，通常使用脉冲响应，因此自适应滤波器利用自适应算法对FIR数字式滤波器3的系数进行调节，使图的评价误差ek＝dk-d^k为最小。
又，向来，为了估计系统的参数(状态)，广泛使用以误差协方差矩阵的更新式(Riccati微分方程式)为依据的卡尔门(Kalman)滤波器。详细情况示于非专利文献2、即S.HaykinAdaptive filter theory，Prentice-Hall(1996)等。
下面对卡尔门滤波器的基本原理进行说明。
如下式所示，以状态空间模型表示的线性系统
xk+1＝ρ-1/2xk，yk＝Hkxk+vk(1)的状态xk的最小方差估计值x^k|k用状态的误差协方差矩阵∑^k|k-1如下所示得到。
x^k|k=x^k|k-1+Kk(yk-Hkx^k|k-1)]]>x^k+1|k=ρ-12x^k|k---(2)]]>Kk=Σ^k|k-1HkT(ρ+HkΣ^k|k-1HkT)-1---(3)]]>Σ^k|k=Σ^k|k-1-KkHkΣ^k|k-1]]>Σ^k+1|k=Σ^k|k/ρ---(4)]]>但是，x^0|-1=0,Σ^0|-1=ϵ0I,ϵ0>0---(5)]]>其中，xk为状态矢量或只是状态；未知，这是估计的对象。
yk为观测信号，是滤波器的输入，已知。
Hk为观测矩阵；已知。
vk为观测噪声；未知。
ρ为忘却系数；通常由尝试错误决定。
Kk为滤波器增益；从矩阵∑^k|k-1得出。
∑^k|k对应于x^k|k的误差的协方差矩阵；利用Riccati微分方程式得出。
∑^k+1|k对应于x^k+1|k的误差的协方差矩阵；利用Riccati微分方程式得出。
∑^1|0对应于初始状态的协方差矩阵；本来是未知的，但是为了方便，使用ε0I。
又，本发明人已经提出了利用高速H∞滤波器的系统同定算法(参照专利文献1)。这是为了系统的同定新决定H∞评价基准，开发以此为依据的超高速H∞滤波器的高速算法，同时以该高速H∞滤波器算法为依据的高速时变系统同定方法。高速H∞滤波器算法能够跟踪每单位时间步骤(step)在计算量O(N)急剧变化的时变系统。在上限值的极限与高速卡尔门(Kalman)滤波器算法完全一致。利用这样的系统同定能够实现时不变以及时变系统的高速实时时间同定估计。
还有，作为系统估计的领域中通常知道的方法，应该参照例如非专利文献2、3。
在回波消除器的应用例在国际电话等长途电话线路中，由于信号的放大等理由，采用4线式线路。另一方面，用户线由于距离比较短，采用2线式线路。
图9是通信系统与回波的说明图。在2线式线路与4线式线路的连接部，如图所示引入混合变压器，进行了阻抗匹配。该阻抗匹配如果是完全的，则来自说话人B的信号(声音)只到达说话人A。但是通常实现完全匹配很难，接受信号的一部分泄漏到4线式线路，发生在放大之后再度返回接收者(说话人A)的现象。这就是回波现象(echo)。回波现象随着传输距离的加长(延迟时间的加长)而影响变大，导致通话质量明显变差(在脉冲传输中，即使是近距离对通话质量变坏的影响也很大)。
图10表示回波消除器的原理图。
其中，如图所示导入回波消除器(echo canceller)，用能够直接观察的接受信号与回波逐次估计回波通道的脉冲响应，通过从实际回波减去利用其得到的拟似回波消除回波，谋求将其去除。
回波通道的脉冲响应的估计的进行要使残留回波ek的均方根误差为最小。这时妨碍回波通道估计的要素是线路噪声和来自说话人A的信号(声音)。通常两个说话人同时开始说话时中断脉冲响应的估计。又，混合变压器的脉冲响应程度为50ms左右，因此取样周期采用125微秒时回波通道的脉冲响应的次数实际上为400左右。
非专利文献11993电子情报通信学会“数字信号处理手册”非专利文献2S.HaykinAdaptive filter theory，Prentice-Hall(1996)非专利文献3B.Hassibi，A.H.Sayed，and T.Kailath“Indenfinite-Quadratic Estimation andControl”SIAM(1996)专利文献1日本特开2002-135171号公报发明内容但是，式(1)～(5)那样的加入已有的忘却系数ρ的卡尔门(Kalman)滤波器中，忘却系数ρ的值必须根据尝试错误决定，非常麻烦。而且忘却系数ρ的值是否果真是最佳值也没有判定手段。
又，在卡尔门滤波器使用的误差协方差矩阵，本来与非零的任意矢量的2次形式通常是正(以下称为“正定”)，但是在利用计算机以单精度进行计算的情况下，其2次形式为负(以下称为“负定”)，可知数值上是不稳定的。而且，计算量为O(N2)(或者O(N3))，因此状态矢量xk的幂N大的情况下，每1步骤的运算次数急剧增大，对实时处理不适合。
本发明鉴于上述存在问题，其目的在于确立能够在理论上最合适地决定忘却系数的估计方法，同时研究出其数值上稳定的估计算法和高速度的算法。又，本发明的目的在于提供能够适用于通信系统和音响系统的回波消除器、声场的重放或噪声控制等的系统估计方法。
本发明为了解决上述存在问题，应用新考虑的H∞最佳化方法，导出能够最佳地决定忘却系数的状态估计算法。而且取代应该是正定的误差的协方差矩阵，通过更新其因数矩阵，研究出数值上稳定的估计算法和高速算法。
本发明的第1种解决手段提供一种系统估计方法、使计算机执行各步骤用的系统估计程序、以及记录该程序的计算机可读记录媒体，是对下式表示的状态空间模型，xk+1＝Fkxk+Gkwkyk＝Hkxk+vkzk＝Hkxk其中，xk为状态矢量或只是状态wk是系统噪声vk是观测噪声yk是观测信号zk是输出信号Fk是系统的动力学Gk是驱动矩阵作为评价基准，从用忘却系数ρ加权的干扰决定，而且将向滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γf的项小的估计算法中，同时进行状态估计的可靠性和忘却系数ρ的最佳化用的系统估计方法、程序以及记录该程序的计算机可读记录媒体，包含下述步骤，即处理部从存储部或输入部输入包含上限值γf、作为滤波器的输入的观测信号yk、观测矩阵Hk的值的步骤、处理部根据上述上限值γf决定关于状态空间模型的忘却系数ρ的步骤、处理部从存储部读取包含初始值或某一时刻的观测矩阵Hk的值，用所述忘却系数ρ执行下式表示的超高速H∞滤波的步骤、
x^k|k＝Fk-1x^k-1|k-1+Ks，k(yk-HkFk-1x^k-1|k-1)其中，x^k|k用观测信号yo～yk的时刻k的状态xk的估计值Fk系统的动力学Ks，k滤波器增益处理部存储关于超高速H∞滤波器的求得的值的存储步骤、处理部根据观测矩阵Hi与滤波器增益Ks，i，计算以所述上限值γf及所述忘却系数ρ为依据的存在条件的步骤、以及处理部使γf小下去，反复执行所述超高速H∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部的步骤。
又，本发明的第2种解决手段提供一种系统估计装置，是对下式表示的状态空间模型，xk+1＝Fkxk+Gkwkyk＝Hkxk+vkzk＝Hkxk其中，xk为状态矢量或只是状态wk是系统噪声vk是观测噪声yk是观测信号zk是输出信号Fk是系统的动力学Gk是驱动矩阵作为评价基准，从用忘却系数ρ加权的干扰决定，而且将向滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γf的项小的估计算法中，同时进行状态估计的可靠性和忘却系数ρ的最佳化用的系统估计装置，具备执行估计算法的处理部、以及利用所述处理部进行读出和/或写入，存储与状态空间模型相关的各观测值、设定值、估计值的存储部，所述处理部从存储部或输入部输入包含上限值γf、作为滤波器的输入的观测信号yk、观测矩阵Hk的值、所述处理部根据上述上限值γf决定关于状态空间模型的忘却系数ρ、
所述处理部从存储部读取包含初始值或某一时刻的观测矩阵Hk的值，用所述忘却系数ρ执行下式表示的的超高速H∞滤波、x^k|k＝Fk-1x^k-1|k-1+Ks、k(yk-HkFk-1x^k-1|k-1)其中，x^k|k用观测信号yo～yk的时刻k的状态xk的估计值Fk系统的动力学Ks，k滤波器增益所述处理部存储关于超高速H∞滤波器的求得的值、所述处理部根据观测矩阵Hi与滤波器增益Ks，i，计算以所述上限值γf及以所述忘却系数ρ为依据的存在条件、而且所述处理部使γf小下去，反复执行所述超高速H∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部。
本发明的估计方法能够最佳地决定忘却系数，而且算法即使是单精度也能够稳定工作，因此能够以低成本实现高性能。通常，一般的民间通信设备等往往从成本和速度方面以单精度进行计算。因此本发明作为实用性的状态估计算法在各种工业领域都能够带来其效果。


图1是本发明实施形态的硬件的结构图。
图2是H∞滤波器的可靠性和忘却系数最佳化的流程图。
图3是图2中的H∞滤波器(S105)的算法的流程图。
图4是定理2的平方根阵列算法的说明图。
图5是定理3的数值上稳定的高速算法的流程图。
图6表示脉冲响应{hi}I＝023的值。
图7是定理3的数值上稳定的高速算法得到的脉冲响应的估计结果。
图8是用于系统估计的结构图。
图9是关于通信系统与回波的说明图。
图10表示回波消除器的原理图。
最佳实施方式下面对本发明的实施形态进行说明。
1.记号的说明首先对本发明实施形态中使用的主要记号及其已知或未知进行说明。
xk状态矢量或只是状态，未知，这是估计的对象。
x0初始状态，未知。
Wk系统噪声，未知。
vk观测噪声，未知。
yk观测信号；作为滤波器的输入，已知。
zk输出信号；未知。
Fk系统的动力学；已知。
Gk驱动矩阵；执行时已知。
Hk观测矩阵；已知。
x^k|k用观测信号yo～yk的时刻k的状态xk的估计值；由滤波器方程式提供。
x^k+1|k用观测信号yo～yk的时刻k+1的状态xk+1的估计值；由滤波器方程式提供。
x^0|0状态的初始估计值；本来未知，但是为了方便使用0。
∑^k|k对应于x^k|k的误差的协方差矩阵；由Riccati微分方程式提供。
∑^k+1|k对应于x^k+1|k的误差的协方差矩阵；由Riccati微分方程式提供。
∑^1|0对应于初始状态的协方差矩阵；本来未知，但是为了方便使用ε0I。
Ks，k滤波器增益；从矩阵∑^k|k-1得到。
ρ忘却系数；在定理1～3的情况下，如果γ决定，则根据ρ＝1-χ(γf)可自动决定ρ。
ef，i滤波器误差Re，k辅助变量还有，记号上所附加的“^”“ˇ”是估计值的意思。又，“～”、“-”、“∪”等是为了方便附加的记号。这些记号，为了输入方便，记载于文字的右上方，但是如数学式所示，与记载于文字的正上方是一样的。而且x、w、H、G、K、R、∑等是矩阵，如数学式所示用粗体字记载，但是为了输入的方便，用普通的文字记载。
2.系统估计的硬件和程序本系统估计方法或系统估计装置·系统，可以利用使计算机执行其各步骤用的系统估计程序、记录系统估计程序的计算机可读记录媒体、包含系统估计程序，能够下载于计算机的内部存储器的程序产品、包含该程序的服务器等的计算机等提供。
图1是有关本实施形态的硬件的结构图。
该硬件具有作为中央处理器(CPU)的处理部101、输入部102、输出部103、显示部104以及存储部105。又，处理部101、输入部102、输出部103、显示部104以及存储部105用星形连接或总线连接等适当的连接手段连接。存储部105根据需要存储系统估计的「1.记号的说明」所示的已知的数据。又，未知·已知的数据和计算出的关于超高速H∞滤波器的数据·其他数据，利用处理部101根据需要写入和/或读出。
3.能够最合适决定忘却系数的超高速H∞滤波器定理1考虑下式那样的状态空间xk+1＝Fkxk+Gkwk， wk，xk∈RN(6)yk＝Hkxk+vk， yk，vk∈R(7)zk＝Hkxk，zk∈R，Hk∈R1×N，k＝0，1，...，L (8)对这样的状态空间模型，提出下式所示的H∞评价基准。
满足该H∞评价基准的状态估计值x^k|k(或输出估计值zˇk|k)，利用下面的电平γf的超高速H∞滤波器提供。
x^k|k=Fk-1x^k-1|k-1+Ks,k(yk-HkFk-1x^k-1|k-1)---(11)]]>Ks,k=Σ^k|k-1HkT·(HkΣ^k|k-1HkT+ρ)-1---(12)]]>Σ^k|k=Σ^k|k-1-Σ^k|k-1CkTRe,k-1CkΣ^k|k-1Σ^k+1|k=(FkΣ^k|kFkT)/ρ---(13)]]>其中，
Re,k=Rk+CkΣ^k|k-1CkT,Rk=ρ00-ργf2,Ck=HkHk---(14)]]>0＜ρ＝1-χ(γf)≤1，γf＞1 (15)还有，式(11)表示滤波器方程式，式(12)表示滤波器增益，式(13)表示Riccati微分方程式。
又，驱动矩阵Gk如下所述生成。
GkGkT=χ(γf)ρFkΣ^k|kFkT---(16)]]>又，为了提高上述高速H∞滤波器的跟踪能力，上限值γf尽可能设定得小，以满足下述存在条件。
Σ^i|i-1=Σ^i|i-1-1+1-γf-2ρHiTHi>0,i=0,...,k(17)]]>其中，χ(γf)是满足χ(1)＝1，χ(∞)＝0的γf的单调衰减函数。
定理1的特征在于，状态估计的可靠性和忘却系数ρ的最佳化同时进行。
图2表示H∞滤波器的可靠性与忘却系数ρ的最佳化的流程图。在这里，数据块「EXC＞0」H∞滤波器的存在条件Δγ正实数首先，处理部101从存储部105或输入部102读出或输入上限值γf(S101)。在这一例子中，给予γf＞＞1。处理部101利用式(15)决定忘却系数ρ(S103)。其后，处理部101根据忘却系数ρ执行式(10)～式(13)的超高速H∞滤波(S105)。处理部101对式(17)(或系数式(18))的右边(将其作为EXC)进行计算(S107)，其存在条件如果在所有的时刻得到满足(S109)，则使γf只小Δγ并反复进行相同的处理(S111)。另一方面，在某一γf存在条件不能够得到满足时(S109)，将该γf加Δγ的和作为γf的最佳值γfop输出到输出部103并且/或者存储于存储部105(S113)。还有，在本例子中是加Δγ，但是也可以加除此以外的预先设定的值。该最佳化的过程被称为γ-迭代。还有，处理部101也可以根据需要把H∞滤波计算步骤S105和存在条件的计算步骤S107等各步骤求得的适当的中间值和最终值适当存储于存储部105，也可以从存储部105读出。
在超高速H∞滤波器满足存在条件时，式(9)的不等式通常被满足。因此，在使(9)的分母的干扰能量受到限定的情况下，式(9)的分子的平方估计误差的总和是有限的，某一时刻以后的估计误差为0。这意味着如果能够使γf更小，则估计值x^k|k能够快速跟踪状态xk的变化。
在这里，希望能够注意定理1的超高速H∞滤波器的算法与通常的H∞滤波器的算法的不同。而且，在γf→∞时，ρ＝1，Gk＝0，定理1的H∞滤波器的算法与卡尔门(Kalman)滤波器的算法一致。
图3表示图2中的(超高速)H∞滤波器(S105)的算法的流程图。
超高速H∞滤波器算法可以摘要说明如下。
〔步骤S201〕处理部101从存储部105读出递归式的初始条件，或者，从输入部102输入初始条件，如图所示那样决定。L表示预先决定的最大数据数。
〔步骤S203〕处理部101将时刻k与最大数据数L加以比较。处理部101在时刻k大于最大数据数时结束处理，如果在最大数据数以下则进入下面的步骤。(如果不要，可以取下条件文。又，也可以根据需要再度开始。)〔步骤S205〕处理部101利用式(12)计算滤波器增益Ks，k。
〔步骤S207〕处理部101更新式(11)的超高速H∞滤波器的滤波器方程式。
〔步骤S209〕处理部101用式(13)的Riccati微分方程式计算与误差的协方差矩阵对应的项∑^k|k、∑^k+1|k。
〔步骤S211〕使时刻k前进(k＝k+1)，返回步骤S203，只要有数据就继续。
还有，处理部101根据需要适当在存储器105存储在步骤H∞滤波器计算步骤S205～S209等各步骤求得的合适的中间值和最终值、存在条件等，也可以从存储部105将其读出。
定标器存在条件而式(17)的存在条件的判定需要O(N2)的计算量。但是如果采用下述条件，则根据计算量O(N)能够验证定理1的H∞滤波器的存在、即式(9)。
系1.定标器存在条件如果采用下述存在条件，则能够用计算量O(N)判定超高速H∞滤波器的存在。
其中，
其中，Ks，i是用式(12)求出的滤波器增益。
证明以下对系1的证明进行说明。
如果求解2×2的矩阵Re，k的特征方程式，即|λI-Re,k|=λ-(ρ+HkΣ^k|k-1HkT)-HkΣ^k|k-1HkT-HkΣ^k|k-1HkTλ-(-ργf2+HkΣ^k|k-1HkT)]]> 则能够得到如下的Re，k的固有值λi。
其中， 如果 则矩阵Re，k的两个固有值中的一个为正，另一个为负，矩阵Rk与Re，k具有相同的贯量。因此，如果采用HkΣ^k|k-1HkT=HkK~k1-1-γf-2ρHkK~k,HkK~k=HkρKs,k1-γf-2HkKs,k]]>则能够得到式(18)的存在条件。其中，HkKs，k的计算量为O(N)。
4.数值上稳定的状态估计算法上述超高速H∞滤波器为了更新∑^k|k+1∈RN×N，每单位时间步骤的计算量为O(N2)、即必须进行正比于N2的算术运算。其中，N为状态矢量xk的幂。因此随着xk的幂的增加，本滤波器的执行需要的计算时间急剧增加。而且误差的协方差矩阵∑^k|k+1根据其性质，通常必须是正定的，但是在数值上也有成为负定的情况。特别是在以单精度计算的情况下，这种倾向显著。已经知道，这时滤波器不稳定。因此，为了算法的实用化和低成本化，最好是研究出在单精度(例如32位)也能够工作的状态估计算法。
在这里，着眼于Rk＝R1/2kJ1R1/2k、Re，k＝R1/2e，kJ1R1/2e，k、∑^k|k-1＝∑^1/2k|k-1∑^1/2k|k-1将数值上稳定化的定理1的H∞滤波器(平方根阵列算法)示于定理2。但是，在这里为了简单，以Fk＝I，但是即使Fk≠I的情况下也同样能够求得。下面表示出实现数值上稳定的状态估计算法用的超高速H∞滤波器。
定理2x^k|k=x^k-1|k-1+Ks,k(yk-Hkx^k-1|k-1)---(20)]]>Ke,k=Kk(:,1)/Re,k(1,1),Kk=ρ12(ρ-12KkRe,k-12J1-1)J1Re,k12---(21)]]> 其中，Rk=Rk12J1Rk12,Rk12=ρ1200ρ12γf,J1=100-1,Σ^k|k-1=Σ^k|k-112Σ^k|k-112]]> Θ(k)为J-酉矩阵，即满足Θ(k)JΘ(k)T＝J，J＝(J1_I)、I是单位矩阵。又，Kk(，1)表示矩阵Kk的第1列的列矢量。
还有，在式(21)、(22)中，J1-1和J1可删除。
图4是定理2的平方根阵列算法的说明图。该计算的算法可以在图2所示的定理1的流程图中的H∞滤波器的计算(S105)中使用。
本估计算法用以J-酉变换为依据的更新式求其因数矩阵∑^1/2k|k-1∈RN×N(∑^k|k-1的平方根矩阵)，代替用Riccati型的更新式求∑^k|k-1。从这时产生的1-1数据块(block)矩阵和2-1数据块矩阵如图所示求滤波器增益Ks，k。因此，∑^k|k-1＝∑^1/2k|k-1∑^1/2k|k-1＞0，∑^k|k-1的正定性得到保证，在数值上能够稳定。还有，定理2的H∞滤波器的每单位步骤的计算量仍然保持O(N2)不变。
还有，在图4中，J1-1可删除。
首先，处理部101从存储部105读出或从内部存储器等得到式(22)左边的矩阵的各要素中包含的项，执行J-酉变换(S301)。处理部101根据式(21)从求得的式(22)的右边的矩阵的要素计算系统增益Kk、Ks，k(S303、S305)。处理部101根据式(20)计算状态估计值x^k|k(S307)。
5.状态估计用的数值上稳定的高速算法如上所述，定理2的H∞滤波器的每单位步骤的计算量仍然是O(N2)。因此作为计算量的对策，在Hk＝Hk+1Ψ、Hk＝〔u(k)、...、u(o)、o、...、o〕时，考虑利用xk＝〔xTK、oT〕T的一步骤预测的误差的协方差矩阵∑k+1|k满足下式(24)，即Σ‾k+1|k-ΨΣ‾k|k-1ΨT=-L‾kRr,k-1L‾kT,L‾k=L~k0---(24)]]>更新幂小的Lk(即L~k)，而不更新∑k+1|k。在这里，如果注意到Rr,k=Rr,k12SRr,k12,]]>则能够得到如下所述的定理3。
定理3x^k|k=x^k-1|k-1+K‾s,k(yk-Hkx^k-1|k-1)---(61)]]>K‾s,k=K‾k(:,1)/Re,k(1,1),K‾k=ρ12(K‾kRe,k-12)Re,k12---(62)]]> 在这里，Θ(k)为任意J-酉矩阵，Cˇk＝Cˇk+1ψ成立。其中Rk=Rk12J1Rk12,Rk12=ρ1200ρ12γf,J1=100-1,Σ^k|k-1=Σ^k|k-112Σ^k|k-112]]> 还有，定理3的证明将在下面叙述。
上式也可以将K-，k(＝P-1/2Kk)代之以Kk进行整理。
而且如果使用下面的J-酉矩阵
Θ(k)=(J1Re,k12⊕-Rr,k12)Σ(k)(Re,k+1-12J1-1⊕-Rr,k+1-12)]]>，则能够得到定理4的高速化的状态估计算法。其中ψ表示移位矩阵。
定理4x^k|k=x^k-1|k-1+Ks,k(yk-Hkx^k-1|k-1)---(25)]]>Ks,k=ρ12K‾k(:,1)/Re,k(1,1)---(26)]]> 其中， diag{·}表示对角矩阵，Re，k+1(1、1)表示矩阵Re，k+1的1-1分量。又，上式也可以将K-k代之以Kk进行整理。
本高速算法利用下面的因式分解，即利用Σ‾k+1|k-ΨΣ‾k|k-1ΨT=-L‾kRτ,k-1L‾kT---(32)]]>中的L~k∈R(N+1)×2的更新求滤波器增益Ks，k，因此每单位步骤的计算量为O(N+1)即可。在这里请求、注意下式。
图5表示定理3的数值上稳定的高速算法的流程图。该高速算法加入图2的H∞滤波器的计算步骤(S105)中，利用γ迭代实现最佳化。因此在存在条件得到满足时，γf慢慢减小，而在没有得到满足的时刻如图所示γf增加。
H∞滤波器算法可以摘要叙述如下。
〔步骤S401〕处理部101如图所示那样决定递归式的初始条件。又，L表示最大数据数。
〔步骤S403〕处理部101将时刻k与最大数据数L加以比较。处理部101在时刻k大于最大数据数时结束处理，如果在最大数据数之下则进入下面的步骤。(如果不要，可以取下条件文。又，也可以再度开始。)〔步骤S405〕处理部101利用式(27)、(31)递归地计算与滤波器增益对应的项Kk+1。
〔步骤S406〕处理部101利用式(29)递归地计算Re，k+1。
〔步骤S407〕处理部101还利用式(26)、(31)计算Ks，k。
〔步骤S409〕处理部101如果在这里判定存在条件EXC＞0，使存在条件满足，就进入步骤S411。
〔步骤S413〕另一方面，如果在步骤S409不满足存在条件，处理部101就增加γf，返回步骤S401。
(步骤S411〕处理部101更新式(25)的H∞滤波器的滤波器方程式。
〔步骤S415〕处理部101利用式(30)递归地计算Rr，k+1。又，处理部101利用式(28)、(31)递归地计算L~k+1。
〔步骤S419〕处理部101使时刻k前进(k＝k+1)。返回步骤S403，只要有数据就继续。
还有，处理部101根据需要适当在存储部105存储H∞滤波器计算步骤S405～S415以及存在条件的计算步骤S409等各步骤求得的适当的中间值以及最终值，也可以将其从存储部105读出。
6.回波消除器下面建立回波消除问题的数学模型。
首先，如果考虑接收信号{uk}为向回波通道的输入信号的情况，借助于回波通道的(时变)脉冲响应{hi〔k〕}，回波{dk}的观测值用下式表示。
yk=dk+vk=Σi=0N-1hi[k]uk-i+vk,k=0,1,2,···---(33)]]>
在这里，uk、yk分别表示时刻tk(＝KT；T是抽样周期)的接收信号与回波，vk表示时刻tk的平均值0的线路噪声，hi〔k〕、i＝0、...、N-1是时变脉冲响应，其抽头数N是已知的。这时，脉冲响应的估计值{h^i(k〕}如果能够时时刻刻得到，则能够如下所述使用其生成拟似回波。
d^k=Σi=0N-1h^k[k]uk-i,k=0,1,2···---(34)]]>如果将其从回波中减去(yk-d^k_0)，则可以消除回波。但是在k-i＜0时，使uk-1＝0。
根据上面所述，问题归结为从可直接观测的接收信号{uk}与回波{yk}逐次估计回波通道的脉冲响应{hi〔k〕}的问题。
通常，为了把H∞滤波器使用于回波消除器，首先必须用状态方程式与观测方程式构成的状态空间模型表现式(32)。因此，由于问题是估计脉冲响应{hi〔k〕}，所以以{hi〔k〕}为状态变量xk，如果允许wk左右的变动，则能够对回波通道使下述状态空间模型成立。
xk+1＝xk+Gkwk， xk，wk∈RN(35)yk＝Hkxk+vk， yk，vk∈R(36)zk＝Hkxk，zk∈R，Hk∈R1×N(37)其中，xk＝[h0[k]，…，hN-1[k]]T，wk＝[wk(1)，…，wk(N)]THk＝[uk，…，uk-N+1]对这样的状态空间模型的超高速和高速H∞滤波器算法如上所述。又，估计脉冲响应时，通常检测发送信号的发生并且在其间中止估计。
7.对脉冲响应的评价动作的确认回波通道的脉冲响应在时间上不变(hi〔k〕＝hi)，而且其抽头数N为48的情况下，用仿真对本高速算法进行确认。
yk=Σi=047hiuk-i+vk---(38)]]>另外，图6表示脉冲响应{hi}的值。
在这里，脉冲响应{hi}i＝023采用图示的值，其他{hi}i＝2447采用0。又，vk是平均值0，方差σv2＝1.0×10-6是恒定的高斯白噪声，为了方便，抽样周期T取1.0。
又，接收信号{uk}如下所示用2次的AR模型近似。
uk＝α1uk-1+α2uk-2+wk’其中α1＝0.7，α2＝0.1，wk’是平均值0，方差σw’2＝0.04的恒定的高斯白噪声。
脉冲响应的估计结果图7表示采用定理3的数值上稳定的高速算法得到的脉冲响应的估计结果。在这里，图7(b)的纵轴表示{Σi=047(hi-x^k(i+1))2}]]>由此可知，利用本高速算法能够进行良好的估计。但是，采用ρ＝1-χ(γf)、χ(γf)＝γf-2、x^0|0＝0、∑^1|0＝20I，计算结果以加倍的精度进行。在确认存在条件的同时，设定γf＝5.5。
8.定理2的证明以下关系式成立时，Rk12CkΣ^k|k-1120ρ-12Σ^k|k-112JRk120Σ^k|k-112CkTρ-12Σ^k|k-112]]>=Re,k120ρ-12KkRe,k-12J1-1Σ^k+1|k12JRe,k12ρ-12J1-1Re,k-12KkT0Σ^k+1|k12---(40)]]>将两边的2×2数据块矩阵的各项加以比较，得到下式。
Re,k=Rk+CkΣ^k|k-1CkT---(41)]]>Kk=Σ^k|k-1CkT---(42)]]>Σ^k+1|k+ρ-1KkRe,k-1KkT=ρ-1Σ^k|k-1---(43)]]>这与定理1的Fk＝I时的式(13)的Riccati微分方程式一致。其中，J=(J1⊕I),J1100-1,Ck=HkHk---(44)]]>另一方面，AJAT＝BJBT成立时，B可以利用J-酉矩阵表达为B＝AΘ(k)。因此，根据式(40)，定理1的Riccati微分方程式与下式等价。
Re,k120ρ-12KkRe,k-12J1-1Σ^k+1|k12=Rk12CkΣ^k|k-1120ρ-12Σ^k|k-112Θ(k)···(45)]]>还有，在式(40)、(45)中，J1-1可删除。
8-2.定理3的证明如下所述，假定数据块三角化的J-酉矩阵Θ(k)存在。
这时，如果将上式的两边J＝(J1_-S)-范数加以比较，可以如下所述决定左边的X、Y、Z。S采用对角分量采取1或-1的对角矩阵。
(1，1)-数据块矩阵XJ1XT]]> =Re,k+(Re,k+1-Rk+1)-(Re,k-Rk)=Re,k+1]]>因此，从Re，k+1＝Re1/2，k+1J1Re1/2，k+1，Rk+1＝Rk，得到X＝Re1/2，k+1。在这里请注意以下关系成立，即J1-1=J1(J12=I),S-1=S,Re,k+1T=Re,k+1,Rr,kT=Rr,k,]]> (2，1)-数据块矩阵
YJ1XT]]> =0K‾k+K‾k+10-0K‾k]]>=K‾k+10]]>据此得到Y=K‾k+10Re,k+1-12J1.]]>其中 (2，2)-数据块矩阵-ZSZT+YJ1YT]]>=0K‾kRe,k-10K‾kT-ρ-1L‾kRr,k-1L~kT]]> 因此， 得到Z=L~k+1Rr,k+1-12]]>8-3.定理4的证明观测矩阵Hk具有移位特性，而且J＝(J1_-S)时，利用与定理2相同的方法能够得到下述关系。

Θ(k)=(J1Re,k12⊕-Rr,k12)Σ(k)(Re,k+1-12J1-1⊕-Rr,k+1-12)]]>其中 ，决定Rr，k+1，使下式成立，即∑(k)T(Re，k_-Rr，k)∑(k)＝(Re，k+1_-Rr，k′+1)接着，如果在式(46)的第3行新追加Rr，k+1的更新式，则最终得到下式。
从该两边的3×2数据块矩阵的各项的对应得到下面的增益矩阵K-k的更新式。
工业应用性通常的民间的通信设备等，往往从成本和速度面上以单精度进行计算。因此，本发明作为实用性的状态估计算法应该能够给各产业领域带来其效果。又，本发明能够使用于通信系统和音响系统的回波消除器、音场重放或噪声控制等。
权利要求
1.一种系统估计方法，是对下式表示的状态空间模型，xk+1＝Fkxk+Gkwk，yk＝Hkxk+vk，zk＝Hkxk，其中，xk为状态矢量或只是状态，wk是系统噪声，vk是观测噪声，yk是观测信号，zk是输出信号，Fk是系统的动力学，Gk是驱动矩阵，作为评价基准，从用忘却系数ρ加权的干扰决定，而且将向滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γf的项小的估计算法中，同时进行状态估计的可靠性和忘却系数ρ的最佳化用的系统估计方法，其特征在于，包含下述步骤，即处理部从存储部或输入部输入包含上限值γf、作为滤波器的输入的观测信号yk、观测矩阵Hk的值的步骤、处理部根据上述上限值γf决定关于状态空间模型的忘却系数ρ的步骤、处理部从存储部读取包含初始值或某一时刻的观测矩阵Hk的值，用所述忘却系数ρ执行下式表示的超高速H∞滤波的步骤、x^k|k＝Fk-1x^k-1|k-1+Ks，k(yk-HkFk-1x^k-1|k-1)其中，x^k|k用观测信号yo～yk的时刻k的状态xk的估计值，Ks，k滤波器增益，处理部存储关于超高速H∞滤波器的求得的值的存储步骤、处理部根据求得的观测矩阵Hi或观测矩阵Hi与滤波器增益Ks，i，计算以所述上限值γf及所述忘却系数ρ为依据的存在条件的步骤、以及处理部使γf小下去，反复执行所述超高速H∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部的步骤。
2.如权利要求1所述的系统估计方法，其特征在于，处理部根据下述公式计算所述存在条件，Σ^i|i-1=Σ^i|i-1-1+1-γf-2ρHiTHi>0,i=0,...,k---(17)]]>
3.如权利要求1所述的系统估计方法，其特征在于，处理部根据下述公式计算所述存在条件， 其中，
4.如权利要求1所述的系统估计方法，其特征在于，所述忘却系数ρ与所述上限值γf存在下式所述的关系，即0＜ρ＝1-χ(γf)≤1，其中，χ(γf)为满足χ(1)＝1，χ(∞)＝0的γf的单调衰减函数。
5.如权利要求1所述的系统估计方法，其特征在于，执行所述超高速H∞滤波的步骤，处理部利用下式求得所述滤波器增益Ks，k x^k|k=Fk-1x^k-1|k-1+Ks,k(yk-HkFk-1x^k-1|k-1)---(11)]]>Ks,k=Σ^k|k-1HkT·(HkΣ^k|k-1HkT+ρ)-1---(12)]]>Σ^k|k=Σ^k|k-1-Σ^k|k-1CkTRe,k-1CkΣ^k|k-1Σ^k+1|k=(FkΣ^k|kFkT)/ρ---(13)]]>其中， Re,k=Rk+CkΣ^k|k-1CkT,Rk=ρ00-ργf2,Ck=HkHk---(14)]]>0＜ρ＝1-χ(γf)≤1，γf＞1 (15)还有，式(16)的右边也可以更加一般化，其中，xk为状态矢量或只是状态，yk是观测信号，zk是输出信号，Fk是系统的动力学，Hk为观测矩阵，x^k|k为用观测信号yo～yk的时刻k的状态xk的估计值，∑^k|k对应于x^k|k的误差的协方差矩阵，Ks，k滤波器增益，ef，i滤波器误差，Re，k辅助变量。
6.如权利要求5所述的系统估计方法，其特征在于，执行所述超高速H∞滤波的步骤包含处理部根据初始条件用所述式(12)计算滤波器增益Ks，k的步骤、处理部更新所述式(11)的H∞滤波器的滤波器方程式的步骤、处理部用所述式(13)计算∑^k|k、∑^k+1|k的步骤、以及处理部使时刻k前进，反复执行所述各步骤的步骤。
7.如权利要求1所述的系统估计方法，其特征在于，执行所述超高速H∞滤波的步骤，其所述处理部利用增益矩阵Kk，使用下式求所述滤波器增益Ks，k，x^k|k=x^k-1|k-1+Ks,k(yk-Hkx^k-1|k-1)---(20)]]>Ks,k=Kk(:,1)/Re,k(1,1),Kk=ρ12(ρ-12KkRs,k-12J1-1)J1Re,k12---(21)]]> 其中，Rk=Rk12J1Rk12,Rk12=ρ1200ρ12γf,J1=100-1,Σ^k|k-1=Σ^k|k-112Σ^k|k-112]]> Θ(k)为J-酉矩阵，即满足Θ(k)JΘ(k)T＝J，J＝(J1_I)、I是单位矩阵。又，Kk(，1)表示矩阵Kk的第1列的列矢量，还有，在这里，在式(21)、(22)中，J1-1和J1可删除，其中，x^k|k为用观测信号yo～yk的时刻k的状态xk的估计值，yk是观测信号，Fk是系统的动力学，Ks，k是滤波器增益，Hk为观测矩阵，∑^k|k对应于x^k|k的误差的协方差矩阵，Θ(k)为J-酉矩阵，Re，k是辅助变量。
8.如权利要求7所述的系统估计方法，其特征在于，执行所述超高速H∞滤波的步骤包含处理部用所述式(22)计算Kk、∑^k+1|k1/2的步骤、处理部根据初始条件用所述式(21)计算滤波器增益Ks，k的步骤、处理部更新所述式(20)的H∞滤波器的滤波器方程式的步骤、以及处理部使时刻k前进，反复执行所述各步骤的步骤。
9.如权利要求1所述的系统估计方法，其特征在于，执行所述超高速H∞滤波的步骤，其所述处理部利用增益矩阵Kk，使用下式求所述滤波器增益Ks，k，x^k|k=x^k-1|k-1+K‾s,k(yk-Hkx^k-1|k-1)---(61)]]>K‾s,k=K‾k(:,1)/Re,k(1,1),K‾k=ρ12(K‾kRe,k-12)Re,k12---(62)]]> 其中，Θ(k)为任意的J-酉矩阵，C^k＝C^k+1Ψ成立。其中，Rk=Rk12J1Rk12,Rk12=ρ1200ρ12γf,J1=100-1,Σ^k|k-1=Σ^k|k-112Σ^k|k-112]]> 其中x^k|k为用观测信号yo～yk的时刻k的状态xk的估计值，yk是观测信号，Ks，k是滤波器增益，Hk为观测矩阵，Θ(k)为J-酉矩阵，Re，k是辅助变量。
10.如权利要求9所述的系统估计方法，其特征在于，执行所述超高速H∞滤波的步骤包含处理部用所述式(63)计算Kk-的步骤、处理部根据初始条件用所述式(62)计算滤波器增益K-s，k的步骤、处理部更新所述式(61)的H∞滤波器的滤波器方程式的步骤、以及处理部使时刻k前进，反复执行所述各步骤的步骤。
11.如权利要求1所述的系统估计方法，其特征在于，执行所述超高速H∞滤波的步骤，其所述处理部利用增益矩阵Kk-，使用下式求所述滤波器增益Ks，k，x^k|k=x^k-1|k-1+Ks,k(yk-Hkx^k-1|k-1)---(25)]]>Ks,k=ρ12K‾k(:,1)/Re,k(1,1)---(26)]]> 其中， 还有，上式也可以用Kk代替Kk-进行整理，K，，其中，yk是观测信号，Fk是系统的动力学，Hk为观测矩阵，x^k|k为用观测信号yo～yk的时刻k的状态xk的估计值，Ks，k为滤波器增益，从增益矩阵Kk-求得，Re，k、L~k辅助变量。
2.如权利要求11所述的系统估计方法，其特征在于，执行所述超高速H∞滤波的步骤包含处理部根据预先决定的初始条件用所述式(27)递归地计算K-k+1的步骤、处理部用所述式(26)计算系统增益Ks，k的步骤、处理部计算存在条件的步骤、如果使所述存在条件满足，处理部就更新所述式(25)的H∞滤波器的滤波器方程式，使时刻k前进，反复执行所述各步骤的步骤、以及如果没有使所述存在条件满足，就增加上限值γf的步骤。
13.如权利要求1所述的系统估计方法，其特征在于，而且利用下式从时刻k的状态估计值x^k|k求输出信号的估计值zˇk|k，zˇk|k＝Hkx^k|k。
14.如权利要求1所述的系统估计方法，其特征在于，使用所述H∞滤波器方程式，求状态估计值x^k|k，对拟似回波如下所述进行估计，用求得的拟似回波抵消实际回波，以此实现回波消除器，d^k=Σi=0N-1h^i[k]uk-i,k=0,1,2,···---(34)]]>
15.一种系统估计程序，是对下式表示的状态空间模型，xk+1＝Fkxk+Gkwk，yk＝Hkxk+vk，zk＝Hkxk，其中，xk为状态矢量或只是状态，wk是系统噪声，vk是观测噪声，yk是观测信号，zk是输出信号，Fk是系统的动力学，Gk是驱动矩阵，作为评价基准，从用忘却系数ρ加权的干扰决定，而且将向滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γf的项小的估计算法中，同时使计算机执行状态估计的可靠性和忘却系数ρ的最佳化用的系统估计程序，其特征在于，使计算机执行下述步骤，即处理部从存储部或输入部输入包含上限值γf、作为滤波器的输入的观测信号yk、观测矩阵Hk的值的步骤、处理部根据上述上限值γf决定关于状态空间模型的忘却系数ρ的步骤、处理部从存储部读取包含初始值或某一时刻的观测矩阵Hk的值，用所述忘却系数ρ执行下式表示的的超高速H∞滤波的步骤、x^k|k＝Fk-1x^k-1|k-1+Ks，k(yk-HkFk-1x^k-1|k-1)其中，x^k|k用观测信号yo～yk的时刻k的状态xk的估计值，Fk系统的动力学，Ks，k滤波器增益，处理部存储关于超高速H∞滤波器的求得的值的存储步骤、处理部根据求得的观测矩阵Hi或观测矩阵Hi与滤波器增益Ks，i计算以所述上限值γf及所述忘却系数ρ为依据的存在条件的步骤、以及处理部使γf小下去，反复执行所述超高速H∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部的步骤。
16.一种记录使计算机执行用的系统估计程序的计算机可读取的记录媒体，是对下式表示的状态空间模型，xk+1＝Fkxk+Gkwk，yk＝Hkxk+vk，zk＝Hkxk，其中，xk为状态矢量或只是状态，wk是系统噪声，vk是观测噪声，yk是观测信号，zk是输出信号，Fk是系统的动力学，Gk是驱动矩阵，作为评价基准，从用忘却系数ρ加权的干扰决定，而且将滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γf的项小的估计算法中，同时使计算机执行状态估计的可靠性和忘却系数ρ的最佳化同时使计算机执行的系统估计程序的记录用的计算机可读取记录媒体，其特征在于，记录使计算机执行下述步骤用的系统估计程序的计算机可读取的记录媒体，所述步骤即处理部从存储部或输入部输入包含上限值γ1、作为滤波器的输入的观测信号yk、观测矩阵Hk的值的步骤、处理部根据上述上限值γf决定关于状态空间模型的忘却系数ρ的步骤、处理部从存储部读取包含初始值或某一时刻的观测矩阵Hk的值，用所述忘却系数ρ执行下式表示的的超高速H∞滤波的步骤、x^k|k=x^k-1|k-1+Ks,k(yk-Hkx^k-1|k-1)---(25)]]>Ks,k=ρ12K‾k(:,1)/Re,k(1,1)---(26)]]> 其中，x^k|k用观测信号yo～yk的时刻k的状态xk的估计值，Fk系统的动力学，Ks，k滤波器增益，处理部存储关于超高速H∞滤波器的求得的值的存储步骤、，，处理部根据求得的观测矩阵Hi或观测矩阵Hi与滤波器增益Ks，j，计算以所述上限值γf及所述忘却系数ρ为依据的存在条件的步骤、以及处理部使γf小下去，反复执行所述超高速H∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部的步骤。
17.一种系统估计装置，是对下式表示的状态空间模型，xk+1＝Fkxk+Gkwk，yk＝Hkxk+vk，zk＝Hkxk，其中，xk为状态矢量或只是状态，wk是系统噪声，vk是观测噪声，yk是观测信号，zk是输出信号，Fk是系统的动力学，Gk是驱动矩阵，作为评价基准，从用忘却系数ρ加权的干扰决定，而且将滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γf的项小的估计算法中，同时进行状态估计的可靠性和忘却系数ρ的最佳化用的系统估计装置，其特征在于，具备执行估计算法的处理部、以及利用所述处理部进行读出和/或写入，存储与状态空间模型相关的各观测值、设定值、估计值的存储部，所述处理部从存储部或输入部输入包含上限值γ1、作为滤波器的输入的观测信号yk、观测矩阵Hk的值、所述处理部根据上述上限值γf决定关于状态空间模型的忘却系数ρ、所述处理部从存储部读取包含初始值或某一时刻的观测矩阵Hk的值，用所述忘却系数ρ执行下式表示的的超高速H∞滤波、x^k|k＝Fk-1x^k-1|k-1+Ks，k(yk-HkFk-1x^k-1|k-1)其中，x^k|k用观测信号yo～yk的时刻k的状态xk的估计值，Fk-1系统的动力学，Ks，k滤波器增益，所述处理部存储关于超高速H∞滤波器的求得的值、所述处理部根据求得的观测矩阵Hi或观测矩阵Hi与滤波器增益Ks，i，计算以所述上限值γf及以所述忘却系数ρ为依据的存在条件、而且所述处理部使γf小下去，反复执行所述超高速H∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部。
全文摘要
本发明确立能最佳地在理论上决定忘却系数的估计方法，同时研究其数值上稳定的估计算法和高速算法。首先，处理部从存储部或输入部读出或输入上限值γ
文档编号G10K11/178GK1836372SQ20048002299
公开日2006年9月20日 申请日期2004年8月5日 优先权日2003年8月11日
发明者西山清 申请人:独立行政法人科学技术振兴机构