系统辨识方法及程式、存储媒体和系统辨识装置

摘要

对大型音响系统或通信系统数值上稳定地进行辨识。当输入信号以M(≤N)阶AR模型表示时，可以计算量3N+O(M)执行高速H∞滤波。处理部决定递归式的初始状态(S201)，根据输入uk设定CUk(S205)，递归地决定变量(S207)，更新矩阵GkN，计算辅助增益矩阵KUkN(S209)，分割辅助增益矩阵(S211)，计算变量DkM和反向预测误差ηM，k(S213)，计算增益矩阵Kk(S215)以及更新高速H∞滤波器的滤波器方程式(S217)。为减少计算的复杂程度，把Kk(：，1)/(1+γf－2HkKk(：，1)直接用作滤波增益Ks，k。

说明书

技术领域

本发明涉及系统辨识方法及程序、存储媒体和系统辨识装置，特别是数值上稳定的高速辨识方法。此外，本发明涉及多种系统辨识方法如利用声音特征为输入信号的大型音响系统或通信系统的高速辨识。

背景技术

系统辨识是指根据输入/输出数据估计系统的输入/输出关系的数学模型(传递函数、脉冲响应等)。代表性的应用例包括国际通信的回波消除器、数据通信的自动均衡器、音响系统的回声消除器和声场重放以及汽车等的主动噪声控制等。迄今为止，作为系统辨识的适应算法，LMS(最小均方)、RLS(递归最小平方)或卡尔门滤波器得到广泛使用。通常，系统输出的观测值表示如下：

[数学式1]

$y_k=\underset{i=0}{\overset{N-1}{\Sigma }}{h}_i{u}_{k-i}+{\upsilon }_k---\left(1\right)$

其中，u_k是输入、h_i是系统的脉冲响应、v_k假定是白噪声。

        详细情况在非专利文献1等中进行说明。

1.LMS

在LMS中，从输入u_k和输出y_k估计系统的脉冲响应x_k＝[h₀、...、h_N-1]^T如下：

[数学式2]

${\hat{x}}_k={\hat{x}}_{k-1}+\mu H_k^T(y_k-H_k{\hat{x}}_{k-1,})---\left(2\right)$

其中，H_k＝[u_k、...、u_k-N+1]^T、μ>0。

2.RLS

在RLS中，从输入u_k和输出y_k估计系统的脉冲响应x_k＝[h₀、...、h_N-1]^T如下：

[数学式3]

${\hat{x}}_k={\hat{x}}_{k-1}+K_k(y_k-H_k{\hat{x}}_{k-1})---\left(3\right)$

$K_k=\frac{P_{k-1}{H_k}^T}{\rho +H_k{P}_{k-1}{H_k}^T}---\left(4\right)$

P_k＝(P_k-1-K_kH_kP_k-1)/ρ           (5)

其中，x＾₀＝0、P₀＝ε₀I、ε₀>0、0是零矢量、I是单位矩阵、K_h是滤波器增益、ρ是忘却系数。(还有，＂＾＂、＂v＂是估计值的意思，应该如数学式所示置于字符正上方，但是为了方便输入，置于字符的右上方。下同。)

3.卡尔门滤波器

用状态空间模型表示的线性系统

[数学式4]

$x_{k+1}={\rho }^{-\frac{1}{2}}{x}_k,y_k=H_k{x}_k+{\upsilon }_k---\left(6\right)$

的状态x_k的最小方差估计值x＾_k|k由以下卡尔门滤波器得到。

[数学式5]

${\hat{x}}_{k|k}={\hat{x}}_{k|k-1}+K_k(y_k-H_k{\hat{x}}_{k|k-1})$

${\hat{x}}_{k+1|k}={\rho }^{-\frac{1}{2}}{\hat{x}}_{k|k}---\left(7\right)$

$K_k={\hat{\Sigma }}_{k|k-1}{H_k}^T{(\rho +H_k{\hat{\Sigma }}_{k|k-1}{H_k}^T)}^{-1}---\left(8\right)$

${\hat{\Sigma }}_{k|k}={\hat{\Sigma }}_{k|k-1}-K_k{H}_k{\hat{\Sigma }}_{k|k-1}$

${\hat{\Sigma }}_{k+1|k}={\hat{\Sigma }}_{k|k}/\rho ---\left(9\right)$

其中，

[数学式6]

${\hat{x}}_{1|0}=0,{\hat{\Sigma }}_{1|0}={ϵ}_{0}I,{ϵ}_0>0---\left(10\right)$

x_k：状态矢量或只是状态；未知，这是估计的对象。

y_k：观测信号；是滤波器的输入，已知。

H_k：观测矩阵；已知。

v_k：观测噪声；未知。

ρ：忘却系数；通常由反复试验决定。

K_k：滤波器增益；从矩阵∑＾_k|k-1得出。

∑＾_k|k：对应于误差为x＾_k|k的协方差矩阵；利用Riccati微分方程式得出。

∑＾_k+1|k：对应于误差为x＾_k+1|k的协方差矩阵；利用Riccati微分方程式得出。

∑＾_1|0：对应于初始状态的协方差矩阵；本来是未知的，但是为了方便，使用ε₀^I。

        另外，迄今为止，专利文献1-2和非专利文献2-5中有对于技术的说明。

专利文献1：WO 02/035727，JP-A-2002-135171

专利文献2：WO 2005/015737

非专利文献1：S.Haykin，Adaptive Filter Theory，3rd Edition，Prentice-Hall，1996

非专利文献2：K.Nishiyama，Derivation of a Fast Algorithm of Modified H∞ Filters，Proceedings of IEEE International Conference on Industrial Electronics，Control andInstrumentation，RBC-II，pp.462-467，2000

非专利文献3：K.Nishiyama，An H∞ Optimization and Its Algorithm for Time-VariantSystem Identification，IEEE Transactions on Signal Processing，52，5，pp.1335-1342，2004

非专利文献4：B.Hassibi，A.H.Sayed，and T.Kailath，Indefinite-Quadratic Estimation andControl，1st Edition，SIAM，1999

非专利文献5：G.Glentis，K.Berberidis，and S.Theodoridis，Efficient least squares adaptivealgorithms for FIR transversal filtering，IEEE Signal Processing Magazine，16，4，pp.13-41，1999

本发明的公开

本发明要解决的问题

目前，系统辨识中使用最广泛的适应算法是LMS。LMS的问题是虽然其计算量小，但收敛速度非常低。另一方面，在RLS或卡尔门滤波器中，支配跟踪性能的忘却系数ρ的值必须由反复试验决定，这通常是件非常困难的任务。而且，忘却系数的确定值是否果真是最佳值也没有判定手段。

此外，虽然∑＾_k|k-1或P_k本来是正定矩阵，但是在以单精度(如32比特)进行计算的情况下，其经常成为负定，这是使卡尔门滤波器或RLS数值上不稳定的因素之一。此外，因为计算量为O(N²)，所以当状态矢量x_k的维数(抽头数)大的情况下，每一步骤的算术运算次数急剧增大，不适用于实时处理。

鉴于上述存在问题，本发明的目的在于提供高速且数值上稳定地辨识大型音响系统或通信系统的辨识方法。此外，本发明的目的在于利用声音特征为输入信号导出一种大大减少之前提出的高速H_∞滤波器的计算量的算法。而且，本发明的目的还在于使用反向预测误差提供在数值上稳定高速H_∞滤波器的方法。

解决问题的手段

根据本发明的第一个解决手段，本发明提供了进行时不变或时变系统的高速实时辨识的系统辨识装置，包括：

一个对干扰具有鲁棒性的滤波器，作为评价基准，规定将由干扰引起的滤波器误差的最大能量增益限制在小于预先给出的上限值γ_f²。

其中，

对于下式(11)-(13)所示的状态空间模型，滤波器符合下式(14)所示的H_∞评价基准，

当用M(≤N)阶自回归模型(AR模型)表示输入信号时，滤波器由下式(38)-(44)给出，以及

滤波器满足下式(45)和(46)的纯量存在(scalar existence)条件：

[数学式7]

x_k+1＝x_k+G_kw_k，w_k，x_k∈R^N               (11)

y_k＝H_kx_k+v_k，y_k，v_k∈R                  (12)

z_k＝H_kx_k，z_k∈R，H_k∈R^1×N               (13)

[数学式8]

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(38\right)$

$K_{s,k}=\frac{K_k(:,1)}{1+{\gamma }_f^{-2}{H}_k{K}_k(:,1)}\in R^{N\times 1}$

$K_k=m_k-D_k{\mu }_k,D_k=\left(\begin{array}{c}{0}_{N-M}\\ D_k^M\end{array}\right)---\left(39\right)$

$D_k^M=\frac{D_{k-1}^M-m_k(N-M+1:N,:)W_{{\eta }_{M,k}}}{1-{\mu }_k{W}_{{\eta }_{M,k}}}$

${\eta }_{M,k}=c_{k-N}+{\underset{‾}{C}}_k^M{D}_{k-1}^M---\left(40\right)$

$S_{M,k}=\rho S_{M,k-1}+e_{M,k}^{T}W{\tilde{e}}_{M,k},e_{M,k}=c_k+C_{k-1}^M{A}_k^M$

$A_k^M=A_{k-1}^M-K_{k-1}(1:M,:)W{\tilde{e}}_{M,k},{\tilde{e}}_{M,k}=c_k+C_{k-1}^M{A}_{k-1}^M---\left(43\right)$

这里，G_k^N更新如下：

$\left(\begin{array}{c}{G}_k^N\\ 0_{1\times 2}\end{array}\right)=\left(\begin{array}{c}{0}_{1\times 2}\\ G_{k-1}^N\end{array}\right)+\left(\begin{array}{c}{e}_{M,k}^T/S_{M,k}\\ A_k^M{e}_{M,k}^T/S_{M,k}\\ 0_{(N-M+1)\times 2}\end{array}\right)$

$-\left(\begin{array}{c}{0}_{(N-M+1)\times 2}\\ e_{M,k-N+M-1}^T/S_{M,k-N+M-1}\\ A_{k-N+M-1}^M{e}_{M,k-N+M-1}^T/S_{M,k-N+M-1}\end{array}\right)---\left(44\right)$

其中，

${\underset{‾}{C}}_k^M=C_k(:,N-M+1:N),$$C_{k-1}^M=C_{k-1}(:,1:M),$$C_k=\left(\begin{array}{c}{H}_k\\ H_k\end{array}\right),$$W=\left(\begin{array}{cc}1& 0\\ 0& -{\gamma }_f^{-2}\end{array}\right)$

$e_{f,i}={z^V}_{i|i}-H_i{x}_i,$${z^v}_{i|i}=H_i{x^{^}}_{i|i}$

(c_k∈R^2×1是C_k＝[c_k、...、c_k-N+1]的第一列矢量，假定c_k-1＝0_2×1、k-i<0，初始值设定为K₀＝0_N×2、G₀^N＝0_(N+1)×2、A₀^M＝0_M×1、S_M，0＝1/ε₀、D₀^M＝0_M×1、x＾_0|0＝x^v0＝0_N×1。这里，0_m×n是m×n零矩阵。)

[数学式9]

这里，分别由下式定义：

${\hat{\Xi }}_i=\frac{{\mathrm{\rho H}}_i{K}_{s,i}}{1-H_i{K}_{s,i}}=\frac{{\mathrm{\rho H}}_i{K}_i(:,1)}{1-(1-{\gamma }_f^{-2})H_i{K}_i(:,1)}---\left(46\right)$

其中

x_k：状态矢量或只是状态；未知，这是估计的对象，

x₀：初始状态；未知，

w_k：系统噪声；未知，

v_k：观测噪声；未知，

y_k：观测信号；是滤波器的输入，已知，

z_k：输出信号，未知，

G_k：驱动矩阵；执行时成为已知，

H_k：观测矩阵；已知，

x＾_k|k：用观测信号y₀-y_k在时刻k时对状态x_k的估计值；由滤波器方程式提供，

K_s，k：滤波器增益；从增益矩阵K_k得到，

ρ：忘却系数；在定理1-3的情况下，当γ_f已决定，则根据ρ＝1-χ(γ_f)可自动决定ρ，

∑₀^-1：表示状态不确定性的加权矩阵的逆矩阵；∑₀为已知，

N：状态矢量的维数(抽头数)；之前给出，

M：AR模型的阶次；之前给出，

μ_k：K^U_k^N的第N+1行矢量；从K^U_k^N得到，

m_k：包含K^U_k^N的第1～N行的N×2矩阵；从K^U_k^N得到，

γ_f：衰减水平；设计时给出，

D_k^M：反向预测系数矢量；从m_k、η_M，k和μ_k得到，

W：加权矩阵；由γ_f决定，

A_k^M：前向预测系数矢量；从K_k-1和e~_M，k得到，

C_k^M：包括C_k的第1～M列矢量的2×M矩阵；从观测矩阵H_k决定，以及

C_k^M：包括C_k的第N-M+1～N列矢量的2×M矩阵；从观测矩阵H_k决定。

     根据本发明的第二个解决手段，本发明提供了进行时不变或时变系统的高速实时辨识的系统辨识方法、促使计算机进行此辨识的系统辨识程序以及记录该程序的计算机可读记录媒体，使用：

一个对干扰具有鲁棒性的滤波器，作为评价基准，规定将由干扰引起的滤波器误差的最大能量增益限制在小于预先给出的上限值γ_f²。

其中，

对于下式(11)-(13)所示的状态空间模型，滤波器符合下式(14)所示的H_∞评价基准，

当用M(≤N)阶自回归模型(AR模型)表示输入信号时，滤波器由下式(38)-(44)给出，以及

滤波器满足下式(45)和(46)的纯量存在条件。

[数学式10]

x_k+1＝x_k+G_kw_k，w_k，x_k∈R^N             (11)

y_k＝H_kx_k+v_k，y_k，v_k∈R                (12)

z_k＝H_kx_k，z_k∈R，H_k∈R^1×N             (13)

[数学式11]

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1,|k-1})$

                                 (38)

$K_{s,k}=\frac{K_k(:,1)}{1+{\gamma }_f^{-2}{H}_k{K}_k(:,1)}\in R^{N\times 1}$

$K_k=m_k-D_k{\mu }_k,D_k=\left(\begin{array}{c}{0}_{N-M}\\ D_k^M\end{array}\right)$

                                   (39)

$D_k^M=\frac{D_{k-1}^M-m_k(N-M+1:N,:)W_{{\eta }_{M,k}}}{1-{\mu }_k{W}_{{\eta }_{M,k}}}$

${\eta }_{M,k}=c_{k-N}+{\underset{‾}{C}}_k^M{D}_{k-1}^M---\left(40\right)$

$S_{M,k}=\rho S_{M,k-1}+e_{M,k}^{T}W{\tilde{e}}_{M,k},e_{M,k}=c_k+C_{k-1}^M{A}_k^M$

$A_k^M=A_{k-1}^M-K_{k-1}(1:M,:)W{\tilde{e}}_{M,k},{\tilde{e}}_{M,k}=c_k+C_{k-1}^M{A}_{k-1}^M---\left(43\right)$

这里，G_k^N更新如下：

$\left(\begin{array}{c}{G}_k^N\\ 0_{1\times 2}\end{array}\right)=\left(\begin{array}{c}{0}_{1\times 2}\\ G_{k-1}^N\end{array}\right)+\left(\begin{array}{c}{e}_{M,k}^T/S_{M,k}\\ A_k^M{e}_{M,k}^T>/S_{M,k}>\\ 0_{(N-M+1)\times 2}\end{array}\right)$

$-\left(\begin{array}{c}{0}_{(N-M+1)\times 2}\\ e_{M,k-N+M-1}^T/S_{M,k-N+M-1}\\ A_{k-N+M-1}^M{e}_{M,k-N+M-1}^T/S_{M,k-N+M-1}\end{array}\right)---\left(44\right)$

其中，

${\underset{‾}{C}}_k^M=C_k(:,N-M+1:N),$$C_{k-1}^M=C_{k-1}(:,1:M),$$C_k=\left(\begin{array}{c}{H}_k\\ H_k\end{array}\right),$$W=\left(\begin{array}{cc}1& 0\\ 0& -{\gamma }_f^{-2}\end{array}\right)$

e_f，i＝z^v_i|i-H_ix_i、z^v_i|i＝H_ix＾_i|i

(c_k∈R^2×1是C_k＝[c_k、...、c_k-N+1]的第一列矢量，假定c_k-1＝0_2×1、k-i<0，初始值设定为K₀＝0_N×2、G₀^N＝0_(N+1)×2、A₀^M＝0_M×1、S_M，0＝1/ε₀、D₀^M＝0_M×1、x＾_0|0＝x^v₀＝0_N×1。这里，0_m×n是m×n零矩阵。)

[数学式12]

这里，分别由下式定义：

${\hat{\Xi }}_i=\frac{{\mathrm{\rho H}}_i{K}_{s,i}}{1-H_i{K}_{s,i}}=\frac{{\mathrm{\rho H}}_i{K}_i(:,1)}{1-(1-{\gamma }_f^{-2})H_i{K}_i(:,1)}---\left(46\right)$

其中

x_k：状态矢量或只是状态；未知，这是估计的对象，

x₀：初始状态；未知，

w_k：系统噪声；未知，

v_k：观测噪声；未知，

y_k：观测信号；是滤波器的输入，已知，

z_k：输出信号，未知，

G_k：驱动矩阵；执行时成为已知，

H_k：观测矩阵；已知，

x＾_k|k：用观测信号y₀-y_k在时刻k对状态x_k的估计值；由滤波器方程式提供，

K_s，k：滤波器增益；从增益矩阵K_k得到，

ρ：忘却系数；在定理1-3的情况下，当γ_f已决定，则根据ρ＝1-χ(γ_f)可自动决定ρ，

∑₀^-1：表示状态不确定性的加权矩阵的逆矩阵；∑₀为已知，

N：状态矢量的维数(抽头数)；之前给出，

M：AR模型的阶次；之前给出，

μ_k：K^U_k^N的第N+1行矢量；从K^U_k^N得到，

m_k：包含K^U_k^N的第1～N行的N×2矩阵；从K^U_k^N得到，

γ_f：衰减水平；设计时给出，

D_k^M：反向预测系数矢量；从m_k、η_M，k和μ_k得到，

W：加权矩阵；由γ_f决定，

A_k^M：前向预测系数矢量；从K_k-1和e~_M，k得到，

C_k^M：包括C_k的第1～M列矢量的2×M矩阵；从观测矩阵H_k决定，以及

C_k^M：包括C_k的第N-M+1～N列矢量的2×M矩阵；从观测矩阵H_k决定。

根据本发明，可以提供一种高速且数值上稳定地辨识大型音响系统或通信系统的辨识方法。此外，根据本发明，可以通过利用声音特征为输入信号而导出一种其计算量比以往提出的高速H_∞滤波器的计算量大大减少的算法。而且，根据本发明，可以通过使用反向预测误差提供在数值上稳定高速H_∞滤波器的方法。

附图说明

图1为显示阶递归高速H_∞滤波器的增益矩阵的更新的示意图。

图2为阶递归高速H_∞滤波器的处理流程图。

图3为纯量存在条件的检查流程图。

图4为在时刻t_k-N+M-1预测变量的值的计算处理流程图。

图5为数值上稳定的高速H_∞滤波器处理的流程图。

图6为通信系统与回波的示意图。

图7为回波消除器的原理图。

图8为显示通过阶递归高速H_∞滤波器的脉冲响应的估计结果(性能评价)的视图。

图9为阶递归高速H_∞滤波器关于声音输入的收敛特性的视图。

图10为存在条件的比较视图。

图11为实施例的硬件结构图。

图12为高速H_∞滤波器和阶递归高速H_∞滤波器的计算量的对比结果图。

图13为显示回波通道的脉冲响应的视图。

最佳实施方式

1.记号说明

首先，说明本实施例中使用的主要记号。

x_k：状态矢量或只是状态；未知，这是估计的对象，

x₀：初始状态；未知，

w_k：系统噪声；未知，

v_k：观测噪声；未知，

y_k：观测信号；是滤波器的输入，已知，

z_k：输出信号，未知，

G_k：驱动矩阵；执行时成为已知，

H_k：观测矩阵；已知，

x＾_k|k：用观测信号y₀-y_k在时刻k对状态x_k的估计值；由滤波器方程式提供，

K_s，k：滤波器增益；从增益矩阵K_k得到，

ρ：忘却系数；在定理1-3的情况下，当γ_f已决定，则根据ρ＝1-χ(γ_f)可自动决定ρ，

∑₀^-1：表示状态不确定性的加权矩阵的逆矩阵；∑₀为已知，

N：状态矢量的维数(抽头数)；之前给出，

M：AR模型的阶次；之前给出，

μ_k：K^U_k^N的第N+1行矢量；从K^U_k^N得到，

m_k：包含K^U_k^N的第1～N行的N×2矩阵；从K^U_k^N得到，

γ_f：衰减水平；设计时给出，

D_k^M：反向预测系数矢量；从m_k、η_M，k和μ_k得到，

W：加权矩阵；由γ_f决定，

A_k^M：前向预测系数矢量；从K_k-1和e~_M，k得到，

C_k^M：包括C_k的第1～M列矢量的2×M矩阵；从观测矩阵H_k决定，以及

C_k^M：包括C_k的第N-M+1～N列矢量的2×M矩阵；从观测矩阵H_k决定。

还有，置于记号上的＂＾＂、＂v＂是估计值的意思。此外，＂～＂、＂-＂、＂U＂等是为了方便而添加的记号。这些记号，为了输入方便，置于文字的右上方，但是如数学式所示，与置于文字正上方是一样的。此外，x、w、H、G、K、R、∑等为矩阵，应该如数学式中所示用粗体字表示，但是为了输入的方便，用普通的字母表示。

2.系统估计的硬件和程序

系统估计方法或系统估计装置/系统可以由使计算机执行其各步骤的系统估计程序、记录系统估计程序的计算机可读记录媒体、包含系统估计程序，能够下载于计算机的内部存储器的程序产品、包含该程序的服务器等的计算机等等提供。

图11为该实施例的硬件结构图。

该硬件包括作为中央处理器(CPU)的处理部101、输入部102、输出部103、显示部104以及存储部105。此外，处理部101、输入部102、输出部103、显示部104以及存储部105用如星形连接或总线连接等适当的连接手段连接。可把经过系统估计的“1.记号说明”中所述的已知数据根据需要存储在存储部105。此外，未知和已知数据、计算出的关于超级H_∞滤波器的数据以及其它数据，利用处理部101根据需要写入和/或读出。

3.和本实施例有关的辨识方法的手段

为了和通常的H_∞滤波器相区别，特别地把在H_∞的含义中可以最合适(准最合适)地决定忘却系数的H_∞滤波器称为超级H_∞滤波器(hyper H_∞filter)(非专利文献2、非专利文献3)。该H_∞滤波器的特点是即使状态的动态未知时也可使用。

定理1(超级H_∞滤波器)

关于状态空间模型

[数学式13]

x_k+1＝x_k+G_kw_k，w_k，x_k∈R^N            (11)

y_k＝H_kx_k+v_k，y_k，v_k∈R               (12)

z_k＝H_kx_k，z_k∈R，H_k∈R^1×N            (13)

满足

的状态估计值x＾_k|k(或输出估计值z^v_k|k)由以下γ_f水平的超级H_∞滤波器给出。

[数学式14]

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(16\right)$

$K_{s,k}={\hat{\Sigma }}_{k|k-1}{H}_k^T{(H_k{\hat{\Sigma }}_{k|k-1}{H}_k^T+\rho )}^{-1}---\left(17\right)$

${\hat{\Sigma }}_{k|k}={\hat{\Sigma }}_{k|k-1}{-{\hat{\Sigma }}_{k|k-1}{C}_k^T{R}_{e,k}^{-1}{C}_k\hat{\Sigma }}_{k|k-1}$

${\hat{\Sigma }}_{k+1|k}={\hat{\Sigma }}_{k|k}/\rho ---\left(18\right)$

其中，

${\hat{\Sigma }}_{1|0}={\Sigma }_0$

$R_{e,k}=R+C_k{\hat{\Sigma }}_{k|k-1}{C}_k^T$

$R=\left(\begin{array}{cc}\rho & 0\\ 0& -\rho {\gamma }_f^2\end{array}\right),C_k=\left(\begin{array}{c}{H}_k\\ H_k\end{array}\right)---\left(19\right)$

0<ρ＝1-χ(γ_f)≤1，γ_f>1            (20)

χ(γ_f)是满足χ(1)＝1及χ(∞)＝0的单调下降函数。此外，驱动矩阵G_k生成如下：

[数学式15]

$G_k{G}_k^T=\chi \left({\gamma }_f\right){\hat{\Sigma }}_{k+1|k}=\frac{\chi \left({\gamma }_f\right)}{\rho }{\hat{\Sigma }}_{k|k}---\left(21\right)$

其中，必须满足以下存在条件。

[数学式16]

${\hat{\Sigma }}_{i|i}^{-1}={\hat{\Sigma }}_{i|i-1}^{-1}+\frac{1-{\gamma }_f^{-2}}{\rho }{H}_i^T{H}_i>0,i=0,...,k---\left(22\right)$

超级H_∞滤波器的特点是状态估计中的鲁棒性的生成和忘却系数ρ的最佳化同时进行。

当超级H_∞滤波器满足该存在条件时，式(14)的不等式总是被满足。因此在式(14)中分母的干扰能量受到限定的情况下，式(14)中分子的平方估计误差的总和是有限的，某一时刻以后的估计误差必须成为0。这意味着如果能够使γ_f更小，则估计值x＾_k|k能够更快地跟随状态x_k的变化。

在这里，需注意定理1的超级H_∞滤波器的算法与通常的H_∞滤波器的算法的不同(非专利文献4)。此外，当γ_f→∞时，则ρ＝1，G_k＝0，定理1的超级H_∞滤波器的算法正式与卡尔门滤波器的算法一致。

超级H_∞滤波器的计算量是O(N²)，这不适合进行实时处理。然而，当观测矩阵H_k具有移位特征H_k+1＝[u_k+1、H_k(1)、H_k(2)、...、H_k(N-1)]时，计算量O(N)的高速算法得到开发(非专利文献2、非专利文献3)。

定理2(高速H_∞滤波器)

当观测矩阵H_k满足移位特征时，可采用下式用计算量O(N)执行具备∑₀＝ε₀I>0的超级H_∞滤波。

[数学式17]

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(23\right)$

$K_{s,k}=\frac{K_k(:,1)}{1+{\gamma }_f^{-2}{H}_k{K}_k(:,1)}\in R^{N\times 1}---\left(24\right)$

K_k＝m_k-D_kμ_k∈R^N×2       (25)

$D_k=[D_{k-1}-m_k{W}_{{\eta }_k}]{[1-{\mu }_k{W}_{{\eta }_k}]}^{-1}$

η_k＝c_k-N+C_kD_k-1          (26)

$S_k=\rho S_{k-1}+e_k^{T}W{\tilde{e}}_k,e_k=c_k+C_{k-1}{A}_k$

$A_k=A_{k-1}-K_{k-1}W{\tilde{e}}_k,{\tilde{e}}_k=c_k+C_{k-1}{A}_{k-1}---\left(29\right)$

其中，$C_k=\left(\begin{array}{c}{H}_k\\ H_k\end{array}\right),$$W=\left(\begin{array}{cc}1& 0\\ 0& -{\gamma }_f^{-2}\end{array}\right)$

c_k∈R^2×1是C_k＝[c_k、...、c_k-N+1]的第一列矢量，假定c_k-i＝0_2×1及k-i<0，初始值设定为k₀＝0_N×2、A₀＝0_N×1、S₀＝1/ε₀、D₀＝0_N×1、X＾_0|0＝0_N×1。其中，0_m×n是m×n零矩阵。

这时，必须满足以下纯量存在条件。

[数学式18]

其中，分别由下式定义。

${\hat{\Xi }}_i=\frac{{\mathrm{\rho H}}_i{K}_{s,i}}{1-H_i{K}_{s,i}}=\frac{{\mathrm{\rho H}}_i{K}_i(:,1)}{1-(1-{\gamma }_f^{-2})H_i{K}_i(:,1)}---\left(31\right)$

4.本发明实施例辨识方法的手段

4.1准备

假定输入信号用M(≤N)阶AR模型(自回归模型)表示，而且已导出了从M×2增益矩阵最佳地生成(N+1)×2辅助增益矩阵K^U_k^N＝Q^U_k^-1C^U_k^T的方法。

其中，Q^U_k用

[数学式19]

表示，且C^U_k＝[C_k c_k-N]＝[c_k C_k-1]成立(非专利文献2)。

[辅助定理1]

当用M(≤N)阶AR模型表示输入信号时，辅助增益矩阵K^U_k^N由下式给出。

[数学式20]

其中，e_M，k＝c_K+C_k-1^MA_k^M、C_k-1^M＝C_k-1(：，1：M)、K^U_k-N+M-1^M-1＝K_k-N+M-1(1：M，：)。这里，C_K-1(：，1：M)是包括C_k-1的第1～M列的矩阵。

(证明)参见“7.辅助定理证明”(之后说明)。

[辅助定理2]

当输入信号用M(≤N)阶AR模型表示时，辅助增益矩阵K^U_k^N等于下式。

[数学式21]

其中，

$\left(\begin{array}{c}{G}_k^N\\ 0_{1\times 2}\end{array}\right)$

$=\left(\begin{array}{c}{0}_{1\times 2}\\ G_{k-1}^N\end{array}\right)+\frac{1}{S_{M,k}}\left(\begin{array}{c}{e}_{M,k}^T\\ A_k^M{e}_{M,k}^T\\ 0_{(N-M+1)\times 2}\end{array}\right)$

$-\frac{1}{S_{M,k-N+M-1}}\left(\begin{array}{c}{0}_{(N-M+1)\times 2}\\ e_{M,k-N+M-1}^T\\ A_{k-N+M-1}^M{e}_{M,k-N+M-1}^T\end{array}\right)---\left(35\right)$

G₀^N＝0_(N+1)×2、K^U_k-N+M-1^M＝0_(M+1)×2、k-N+M-1≤0。

(证明)参见“7.辅助定理证明”(之后说明)。

[辅助定理3]

当输入信号用M(≤N)阶AR模型表示时，增益矩阵K_k用下式表示。

[数学式22]

K_k＝m_k-D_kμ_k          (36)

其中，

$D_k=\left(\begin{array}{c}{0}_{N-M}\\ D_k^M\end{array}\right),$

$D_k^M=\frac{D_{k-1}^M-m_k(N-M+1:N,:)W_{{\eta }_{M,k}}}{1-{\mu }_k{W}_{{\eta }_{M,k}}}$

${\eta }_{M,k}=c_{k-N}+{\underset{‾}{C}}_k^M{D}_{k-1}^M$

m_k∈R^N×2，μ_k∈R^1×2

${\underset{‾}{C}}_k^M=C_k(:,N-M+1:N)---\left(37\right)$

这里，m_k(M：N，：)＝[m_k(M，：)^T、...、m_k(N，：)^T]^T，m_k(i，：)是m_k的第i行矢量。

(证明)参见“7.辅助定理证明”(之后说明)。

图1为显示阶递归高速H_∞滤波器的增益矩阵的更新的说明图。

此图显示了以上概要。以此可知K_k-N+M-1(1：M，：)＝K_k-1(N-M+1：N，：)在理论上成立。

这一点将在下面进行说明。

包含N×2增益矩阵K_k-N+M-1第1～M行的K_k-N+M-1^M是利用A_k-i^M和S_M，k-i外推的，而生成(N+1)×2辅助增益矩阵K^U_k^N后，利用D_k^M得到下一时刻的增益矩阵K_k。以后对此进行重复。

4.2阶递归高速H_∞滤波器的导出

接下来说明带有最佳扩大滤波器增益的高速H_∞滤波器。

定理3(阶递归高速H_∞滤波器)

当输入信号用M(≤N)阶AR模型表示时，由下式用计算量3N+O(M)执行高速H_∞滤波。

[数学式23]

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(38\right)$

$K_{s,k}=\frac{K_k(:,1)}{1+{\gamma }_f^{-2}{H}_k{K}_k(:,1)}\in R^{N\times 1}$

$K_k=m_k-D_k{\mu }_k,D_k=\left(\begin{array}{c}{0}_{N-M}\\ D_k^M\end{array}\right)---\left(39\right)$

$D_k^M=\frac{D_{k-1}^M-m_k(N-M+1:N,:)W_{{\eta }_{M,k}}}{1-{\mu }_k{W}_{{\eta }_{M,k}}}$

${\eta }_{M,k}=c_{k-N}+{\underset{‾}{C}}_k^M{D}_{k-1}^M---\left(40\right)$

$S_{M,k}=\rho S_{M,k-1}+e_{M,k}^{T}W{\tilde{e}}_{M,k},e_{M,k}=c_k+C_{k-1}^M{A}_k^M$

$A_k^M=A_{k-1}^M-K_{k-1}(1:M,:)W{\tilde{e}}_{M,k},{\tilde{e}}_{M,k}=c_k+C_{k-1}^M{A}_{k-1}^M---\left(43\right)$

其中，G^N_k更新如下：

$\left(\begin{array}{c}{G}_k^N\\ 0_{1\times 2}\end{array}\right)=\left(\begin{array}{c}{0}_{1\times 2}\\ G_{k-1}^N\end{array}\right)+\left(\begin{array}{c}{e}_{M,k}^T/S_{M,k}\\ A_k^M{e}_{M,k}^T/S_{M,k}\\ 0_{(N-M+1)\times 2}\end{array}\right)$

$-\left(\begin{array}{c}{0}_{(N-M+1)\times 2}\\ e_{M,k-N+M-1}^T/S_{M,k-N+M-1}\\ A_{k-N+M-1}^M{e}_{M,k-N+M-1}^T/S_{M,k-N+M-1}\end{array}\right)---\left(44\right)$

其中，

${\underset{‾}{C}}_k^M=C_k(:,N-M+1:N),$$C_{k-1}^M=C_{k-1}(:,1:M),$$C_k=\left(\begin{array}{c}{H}_k\\ H_k\end{array}\right),$$W=\left(\begin{array}{cc}1& 0\\ 0& -{\gamma }_f^{-2}\end{array}\right)$

e_f，i＝z^v_i|i-H_ix_i、z^v_i|i＝H_ix＾_i|i，c_k∈R^2×1是C_k＝[c_k、...、c_k-N+1]的第一列矢量，c_k-i＝0_2×1，k-i<0，初始值设定为K₀＝0_N×2、G₀^N＝0_(N+1)×2、A₀^M＝0_M×1、S_M，0＝1/ε₀、D₀^M＝0_M×1、x＾_0|0＝x^v0＝0_N×1。

此外，必须满足以下纯量存在条件。

[数学式24]

这里，分别由下式定义。

${\hat{\Xi }}_i=\frac{{\mathrm{\rho H}}_i{K}_{s,i}}{1-H_i{K}_{s,i}}=\frac{{\mathrm{\rho H}}_i{K}_i(:,1)}{1-(1-{\gamma }_f^{-2})H_i{K}_i(:,1)}---\left(46\right)$

(证明)当把辅助定理1、2和3应用于定理2时得到。

4.3阶递归高速H∞滤波器的装置

图2显示阶递归高速H∞滤波器的流程图。

为减少计算量，不使用K_s，k而直接使用K_k(：，1)/(1+γ_f^-2H_kK_k(：，1))。

下文将参照流程图说明阶递归高速H_∞滤波器的处理。

[步骤S201，初始化]处理部101从存储部105读出递归式的初始状态，或者，从输入部102输入初始状态，并如图所示那样确定。

[步骤S203]处理部101将时刻k和最大数据数L进行比较。还有，L是预先决定的最大数据数。在时刻k大于最大数据数时，处理部101结束处理；如果不大于最大数据数则进入下一步骤。(如果非必要，可以删除条件语句。又，也可以根据需要再度开始。)

[步骤S205，输入]处理部101从输入部102输入一输入u_k，并如图所示设定C^U_k。还有，该输入u_k可从存储部105读出。

[步骤S207，正向预测]处理部101用式(43)递归地决定变量e_M，k、A_k^M、S_M，k、K_k(1：M，：)等。

[步骤S209，外推]处理部101用式(44)更新矩阵G_k^N，用式(42)计算辅助增益矩阵K^U_k^N。

[步骤S211，分割]处理部101用式(41)分割辅助增益矩阵K^U_k^N。

[步骤S213，反向预测]处理部101用式(40)计算变量D_k^M和反向预测误差η_M，k。

[步骤S215，增益矩阵]处理部101用式(39)计算增益矩阵K_k。

[步骤S217，滤波]处理部101以式(38)更新高速H_∞滤波器的滤波器方程式。这里为减少计算量，K_k(：，1)/(1+γ_f^-2H_kK_k(：，1)直接用作滤波器增益K_s，k。

[步骤S219]处理部101使时刻k前进(k＝k+1)，返回至步骤S203，继续至某一指定次数，或当还有数据存在时一直继续。

还有，处理部101根据需要可在存储部105存储在H_∞滤波器计算步骤S205至S217的各步骤求得的合适的中间值和最终值、存在条件的值等，也可以从存储部105将其读出。

此外，应注意当M＝N成立时，阶递归高速H_∞滤波器与高速H_∞滤波器完全一致。

图3为纯量存在条件检查的流程图。

通常，图3的存在条件的检查功能被添加到图2的流程图。可在图2或图4(之后说明)流程图的适当的一个或多个步骤之前或之后或使时刻k前进的步骤(S219)之前或之后的适当时机执行这一步骤。EXC(k)相当于前面提到的式(45)的左边。

4.4预测变量的计算

得到图2的外推(步骤S209)中的变量e_M，k-N+M-1、A_k-N+M-1^M、S_M，k-N+M-1和K^U_k-N+M-1^M-1(＝K_k-N+M-1^M＝K_k-N+M-1(1：M，：)的方法，有以下两个方法。

1)在N-M步骤中把各个步骤得到的变量的值保存在存储器上。

2)在各个步骤也再次计算在时刻t_k-N+M-1的变量的值(参见之后图4所述)。

通常，方法1)适合具备大存储容量的DSP(数字信号处理器)，方法2往往适合普通的DSP，可使用两种方法中的任一方法。

图4是在时刻t_k-N+M-1的预测变量的值的计算方法流程图。

下文将参照流程图说明在时刻t_k-N+M-1的预测变量的值的计算处理。处理部101根据流程图计算各个变量，当计算图2的流程图的步骤S209等时，处理部使用各个变量。

[步骤S401，初始化]处理部101从存储部105读出递归式的起始状态或者从输入部102输入起始状态，如图所示那样确定。

[步骤S403]处理部101将时刻k和最大数据数L进行比较。还有，L是预先决定的最大数据数。在时刻k大于最大数据数时，处理部101结束处理；如果不大于最大数据数则进入下一步骤。(如果非必要，可以删除条件语句。又，也可以根据需要再度开始。)

[步骤S405，输入]从输入部102输入一输入u_k，如图所示设定C^U_k。还有，可从存储部105读出输入u_k。

[步骤S407，正向预测]处理部101如图所示递归地决定变量e_M，k、A_k^M、S_M，k、K_k-1^M等(参见式(43))。

[步骤S409，外推]处理部101如图所示计算辅助增益矩阵K^U_k^M(参见式(42))。

[步骤S411，分割]处理部101如图所示分割辅助增益矩阵K^U_k^M(参见式(41))。

[步骤S413，反向预测]处理部101如图所示计算变量D_k^M和反向预测误差η_M，k(参见式(40))。

[步骤S415，增益矩阵]处理部101如图所示计算增益矩阵K_k^M(参见式(39))。

[步骤S417]处理部101使时刻k前进(k＝k+1)，返回步骤S403，继续至某一指定次数或当还有数据存在时一直继续。

还有，处理部101根据需要可在存储部105存储在H_∞滤波器计算步骤S405至S415等各步骤求得的合适的中间值和最终值、存在条件的值等，也可以从存储部105将其读出。

4.5反向预测误差

为了数值上的稳定，可弃用反向预测误差η_k，改为采用

[数学式25]

${\tilde{\eta }}_k={\eta }_k+\beta ({\eta }_k-{\rho }^{-N}{S}_k{\mu }_k^T)---\left(47\right)$

图5为数值上稳定的高速H_∞滤波器(M＝N时为阶递归高速H_∞滤波器)的流程图。

通常，图3的存在条件的检查功能被添加到此流程图。此外，β表示控制参数，且通常使β为β＝1.0。这种稳定方法也可应用于M<N的情况。

下文将参照流程图说明数值上稳定的高速H_∞滤波器的处理。

[步骤S501，初始化]处理部101从存储部105读出递归式的初始状态或者从输入部102输入初始状态，如图所示那样确定。

[步骤S503]处理部101将时刻k和最大数据数L进行比较。还有，L是预先决定的最大数据数。在时刻k大于最大数据数时，处理部101结束处理；如果不大于最大数据数则进入下一步骤。(如果非必要，可以删除条件语句。又，也可以根据需要再度开始。)

[步骤S505，输入]从输入部102输入一输入u_k，并如图所示设定C^U_k。还有，可从存储部105读出输入u_k。

[步骤S507，正向预测]处理部101如图所示递归地决定变量e_M，k、A_k、S_k、K_k-1等(参见式(43))。

[步骤S509，外推]处理部101如图所示计算辅助增益矩阵K^U_k^N(参见式(42))。

[步骤S511，分割]处理部101如图所示分割辅助增益矩阵K^U_k^N(参见式(41))。

[步骤S513，稳定的反向预测]处理部101如图所示计算变量D_k^M和反向预测误差η~_k(参见式(40))。这里，为了数值上的稳定，未采用反向预测误差η_k，而是采用了η~_k(式(47))。

[步骤S515，增益矩阵]处理部101如图所示计算增益矩阵K_k(参见式(39))。

[步骤S517，滤波]处理部101用式(38)更新高速H_∞滤波器的滤波器方程式。这里，为了减少计算量，把K_k(：，1)/(1+γ_f^-2H_kK_k(：，1)直接用作滤波器增益K_s，k。

[步骤S519]处理部101使时刻k前进(k＝k+1)，返回步骤S503，继续至某一指定次数，或当还有数据存在时一直继续。

还有，处理部101根据需要可在存储部105存储在H_∞滤波器计算步骤S505至S517等各步骤求得的合适的中间值和最终值、存在条件的值等，也可以从存储部105将其读出。

5.比较结果

图12为高速H_∞滤波器(FHF)的计算量和阶递归高速H_∞滤波器(OR-FHF)的计算量之间的对比结果图。

不过，这个结果是在根据以上数式计算的情况下得到的，通过措施可进一步减少计算量。此外，通过取纯量的倒数，矢量除以纯量可计为乘法。

依此可见，当M<<N成立时，阶递归高速H_∞滤波器要求的乘法次数变为约3N，这与高速H_∞滤波器相比，其计算量可大大减少。

6.回波消除器

6.1准备

下文采用回波消除器的实例为实施例，确认本辨识算法的运作。在此之前，将简单地说明回波消除器。

图6为通信系统与回波的说明图。

在国际电话等长途电话线路中，由于信号放大等理由，采用4线式线路。另一方面，用户线由于距离比较短，采用2线式线路。在2线式线路和4线式线路的连接部分，如图所示引入混合变压器，进行了阻抗匹配。该阻抗匹配如果是完全的，则来自说话人B的信号(声音)仅到达说话人A。但是通常很难实现完全匹配，会有接收信号的一部分泄漏到4线式线路，发生在放大之后再度返回接收者(说话人A)的现象。这就是回波现象。回波现象随着传输距离的加长(延迟时间的加长)而影响变大，导致通话质量明显变差(在脉冲传输中，即使是近距离，由于回波而使通话质量变差都可能是一个问题)。

图7为回波消除器的原理图。

如图所示，导入回波消除器，用能够直接观测的接收信号与回波逐次估计回波通道的脉冲响应，用其得到一拟似回波，并将其从实际回波减去以消除回波并将其去除。

进行回波通道的脉冲响应的估计以使残留回波e_k的均方误差为最小。这时，妨碍回波通道估计的要素是线路噪声和来自说话人A的信号(声音)。通常两个说话人同时开始说话(双向通话)时，会暂停脉冲响应的估计。此外，混合变压器的脉冲响应长度为例如50[ms]左右，因此采样周期采用125[μs]时，回波通道的脉冲响应的阶级实际上为400左右。

接下来建立此回波消除问题的数学模型。第一，如果考虑接收信号{u_k}成为回波通道的输入信号，通过回波通道的脉冲反应{h_i[k]}，回波{d_k}的观测值{y_k}用下式表示：

[数学式26]

$y_k=d_k+{\upsilon }_k=\underset{i=0}{\overset{N-1}{\Sigma }}{h}_i{\left[k\right]u}_{k-i}+{\upsilon }_k,k=\mathrm{0,1,2},-···--\left(48\right)$

这里，u_k和y_k分别是该回波在时刻t_k(＝kT；T为采样周期)的接收信号与回波的观测值。v_k是平均值为0的线路噪声在时刻t_k的值。假定抽头数N已知。这时，当脉冲响应的估计值{h＾_i[k]}能够在每个时刻得到，则能够使用如下所述的方法生成拟似回波。

[数学式27]

${\hat{d}}_k=\underset{i=0}{\overset{N-1}{\Sigma }}{\hat{h}}_i{\left[k\right]u}_{k-i},k=\mathrm{0,1,2},-···--\left(49\right)$

当从回波减去d＾_k(y_k-d＾_k≈0)，则可消除回波。其中，假定如果k-i<0，则u_k-i＝0。

根据上面所述，消除问题可简化为从可直接观测的接收信号{u_k}与回波的观测值{y_k}逐次估计回波通道的脉冲反应{h_i[k]}的问题。

通常为了把H_∞滤波器应用于回波消除器，首先，式(48)必须用包含状态方程式和观测方程式的状态空间模型表示。然后，如果使{h_i[k]}成为状态变量x_k并允许出现约为w_k的方差，则可为回波通道建立以下的状态空间模型。

[数学式28]

x_k+1＝x_k+G_kw_k，x_k，w_k∈R^N          (50)

y_k＝H_kx_k+v_k，y_k，v_k∈R             (51)

z_k＝H_kx_k，z_k∈R，H_k∈R^1×N         (52)

其中，

x_k＝[h₀[k]、...、h_N-1[k]]^T、w_k＝[w_k(1)、...、w_k(N)]^T

H_k＝[u_k、...、u_k-N+1]

6.2操作的确认

图13为显示回波通道的脉冲响应的图表。

回波通道的脉冲响应在时间上不变(h_i[k]＝h_i)，且例如抽头数N为48的情况下，用模拟对本高速算法进行确认。在这时，回波的观测值y_k由下式表示。

[数学式29]

$y_k=\underset{i=0}{\overset{63}{\Sigma }}{h}_i{\mu }_{k-i}+{\upsilon }_k---\left(53\right)$

其中，脉冲响应{h_i}_i＝0²³使用图13的值，{h_i}_i＝24⁶³则使为0。此外，假定v_k是平均值为0，方差σ_v＝1.0×10^-4的恒定高斯白噪声，为了方便，采样周期T取1.0。

此外，接收信号{u_k}如下所示用2次的AR模型模拟：

[数学式30]

${\mu }_k={\alpha }_1{\mu }_{k-1}+{\alpha }_2{\mu }_{k-2}+{\omega }_k^{\prime }---\left(54\right)$

其中，α₁＝0.7、α₂＝0.1、而w_k′是平均值为0，方差σ_w′²＝0.04的恒定高斯白噪声。

图8为显示根据M＝2时的阶递归高速H_∞滤波器的脉冲响应的估计结果(性能评价)的视图。

这里，图8A显示了k＝256时脉冲响应的估计值(X＾_256|256)(虚线表示真实值)，图8B显示其收敛特征。其中，图8B的纵轴表示‖x_k-x＾_k|k‖²＝∑_i＝0⁶³(h_i-x＾_k(i+1))²。由此可知，阶递归高速H_∞滤波器可以很好地估计系统的脉冲响应。其中，假定ρ＝1-χ(γ_f)、x(γ_f)＝γ_f^-2、x＾_0|0＝0、γ_f＝10.0、∑＾_1|0＝ε₀I及ε₀＝10.0，计算以MATLAB(双倍精度)进行。

图9为关于声音输入的阶递归高速H_∞滤波器的收敛特征视图(在此例中，N＝512、M＝20、γ_f＝54.0)。

这里，假定未知系统的脉冲响应的阶次为512，声音为第20阶AR模型。此外，假定γ_f＝54.0、ε₀＝1.0，计算以MATLAB(双倍精度)进行。由此可知，即使输入是声音，阶递归高速H_∞滤波器也具有良好性能。

图10为存在条件的对比视图。该图显示关于高速H_∞滤波器(M＝N时的阶递归高速H_∞滤波器)(图10A)以及数值上稳定的高速H_∞滤波器(图10B)(在此例中，N＝256、γ_f＝42.0)的纯量存在条件的左边(EXC(k))值随时间的变化。

这时，可见高速H_∞滤波器在45000步左右时EXC(k)的值变为负数，滤波器停止。另一方面，在数值上稳定的高速H_∞滤波器内，EXC(k)在此部分总是保持正值且满足该存在条件。不过，该计算以单精度进行，各个输入使用相同的声音信号。

7.辅助定理证明

7.1辅助定理1的证明

当状态矢量的维数为M时，从式(28)中得到下式。

[数学式31]

当把此应用于维数为M+1的情况时，得到辅助增益矩阵K^U_k^M+1。

[数学式32]

另一方面，当可用M阶AR模型表示输入信号时，对于m≥M，e_M，k＝min{e_m，k}和e~_M，k＝min{e~_m，k}成立。这时，

[数学式33]

$e_{m,k}=c_k+C_{k-1}^m{A}_k^m=c_k+C_{k-1}^M{A}_k^M=e_{M,k}$

${\tilde{e}}_{m,k}=c_k+C_{k-1}^m{A}_{k-1}^m={\tilde{e}}_{M,k}$

$A_k^m{e}_{m,k}^T=\left(\begin{array}{c}{A}_k^M{e}_{M,k}^T\\ 0_{(m-M)\times 1}\end{array}\right)$

$S_{m,k}=S_{M,k}=\rho S_{M,k-1}+e_{M,k}^{T}W{\tilde{e}}_{M,k}---\left(57\right)$

这时候，

成立。因此得到下式。

[数学式34]

同样地，当重复此至N，最终通过使用A_k-1^M，S_M，k-1满足

[数学式35]

即

，辅助增益矩阵K^U_k^N可表示如下：

[数学式36]

7.2辅助定理2的证明

辅助变量G_k^N定义为

[数学式37]

$G_k^N=\underset{i=0}{\overset{N-M}{\Sigma }}\frac{1}{S_{M,k-1}}\left(\begin{array}{c}{0}_{i\times 2}\\ e_{M,k-i}^T\\ A_{k-i}^M{e}_{M,k-i}^T\\ 0_{(N-M-i)\times 2}\end{array}\right)$

这时，注意G_k^N(i-1，：)＝G_k-1^N(i，：)，其成为

[数学式38]

$\left(\begin{array}{c}{G}_k^N\\ 0_{1\times 2}\end{array}\right)-\left(\begin{array}{c}{0}_{1\times 2}\\ G_{k-1}^N\end{array}\right)$

$=\frac{1}{S_{M,k}}\left(\begin{array}{c}{e}_{M,k}^T\\ A_k^M{e}_{M,k}^T\\ 0_{(N-M+1)\times 2}\end{array}\right)-\frac{1}{S_{M,k-N+M-1}}\left(\begin{array}{c}{0}_{(N-M+1)\times 2}\\ e_{M,k-N+M-1}^T\\ A_{k-N+M-1}^M{e}_{M,k-N+M-1}^T\end{array}\right)$

这时，式(35)成立。因此，辅助增益矩阵K^U_k^N等于式(34)(式(42))。

7.3辅助定理3的证明

当输入信号用M(≤N)阶AR模型表示时，其成为

[数学式39]

${\eta }_k=c_{k-N}+C_k{D}_{k-1}=c_{k-N}+{\underset{‾}{C}}_k^M{D}_{k-1}^M$

当注意到此时，式(36)和式(37)由与非专利文献2类似的方法得到。

工业应用性

本发明的高速辨识方法可应用于大型音响系统或通信系统的辨识，对于扬声器或电视会议系统的回波消除器、主动噪声控制、声场重放等有效果。此外，随着数值上稳定方法的发展，即使以单精度也可使计算更加稳定，可用低成本实现较高性能。

一般而言，通常的民间通信设备等往往因成本和速度考虑而以单精度进行计算。因此，本发明作为实用的高速辨识方法应该能够在各产业领域带来其效果。