信号失真消除装置、方法、程序及记录该程序的记录介质

摘要

本发明的信号失真消除装置包括：反向滤波器应用单元，在满足了规定的反复终止条件时，将反向滤波器应用于观测信号的结果作为复原信号输出，在不满足反复终止条件时，将反向滤波器应用于观测信号的结果作为特定信号输出；预测误差滤波器计算单元，将特定信号进行分帧，输出对各帧的特定信号进行线性预测分析所得的每个帧的预测误差滤波器；反向滤波器计算单元，求对于各帧的特定信号拟合了应用与该帧对应的预测误差滤波器所得的各信号(修正估计值)的全帧中的修正估计值在其全样本间独立的反向滤波器并输出；以及控制单元，使反向滤波器应用单元、预测误差滤波器计算单元、反向滤波器计算单元反复执行，直到满足反复终止条件。

说明书

技术领域

该发明涉及消除信号的失真(Distortion)的技术。

背景技术

若信号是在有反射和混响的环境下被观测，则将作为在原来的信号中重叠了反射和混响等的信号而被观测。以下，原来的信号称为“原信号”，被观测的信号称为“观测信号”。此外，由反射和混响等所代表的重叠在原信号中的失真称为“传递特性”。因此，从观测信号中抽出原信号固有的特征将变得困难。为消除这种不便，很早开始就尝试有各种信号失真消除处理技术。信号失真消除处理是从观测信号中去除重叠在原信号的传递特性的处理。

使用图15说明非专利文献1展示的作为以往的信号失真消除处理的一个例子的信号失真消除方法。预测误差滤波器计算部分(901)将观测信号进行分帧处理，对各帧所包含的观测信号进行线性预测分析，并计算预测误差滤波器。在本说明书中，滤波器为数字滤波器，在求对信号的样本值作用的所谓滤波器系数的意义下，有时仅简称为计算滤波器等。预测误差滤波器应用部分(902)将对每个帧进行所述计算的预测误差滤波器应用于该帧的观测信号。反向滤波器计算部分(903)计算反向滤波器，该反向滤波器能够使对于应用了预测误差滤波器之后的信号应用反向滤波器所得的信号的归一化峰值最大。反向滤波器应用部分(904)通过将经过所述计算的反向滤波器应用于观测信号而得到信号失真消除后的信号(复原信号)。

非专利文献1：B.W.Gillespie，，H.S.Malvar，，and D.A.F.Florencio，，“Speechdereverberation via maximum-kurtosis subband adaptive filtering，”IEEEInternational Conference on Acoustics，Speech，and Signal Processing，pp.3701-3704，2001.

发明内容

发明要解决的课题

上述以往的信号失真消除方法假设，观测信号的各帧内的样本间的相关是原信号固有特性的作用较大，而横跨帧的样本间的相关是传递特性的作用较大。上述以往方法基于该假设，将预测误差滤波器应用于经过分帧处理的观测信号，从而降低观测信号中的原信号固有特性的作用。

但是，该假设是粗略的近似，因此反向滤波器的精度是不充分的。即，观测信号求得的预测误差滤波器受到传递特性的影响，因此，无法正确地仅去除原信号固有的特性。因此，由预测误差滤波器应用后的信号所求得的反向滤波器的精度下降。结果，对观测信号应用反向滤波器所得的信号在与原本的原信号的比较下，其精度并不高。

因此，本发明的目的在于，通过求精度高的反向滤波器，从观测信号中消除由来于传递特性的失真而得到精度高的复原信号。

解决课题的方案

为了解决所述课题，本发明的信号失真消除装置包括：反向滤波器应用单元，在满足了规定的反复终止条件时，将适用于观测信号的滤波器(以下，称为反向滤波器)应用于观测信号，并将该结果作为复原信号输出，在不满足反复终止条件时，将反向滤波器应用于观测信号，并将该结果作为特定信号输出；预测误差滤波器计算单元，将特定信号进行分帧，输出对各帧的特定信号进行线性预测分析所得的每个帧的预测误差滤波器；反向滤波器计算单元，求对于各帧的特定信号拟合了应用与该帧对应的预测误差滤波器所得的各信号(以下，称为修正估计值)的全帧中的修正估计值(以下，称为修正估计值系列)在其全样本间独立的反向滤波器，并输出该反向滤波器；以及控制单元，使反向滤波器应用单元、预测误差滤波器计算单元、反向滤波器计算单元反复执行，直到满足反复终止条件。

在该发明中，对于将用于消除传递特性的反向滤波器应用于观测信号所得特定信号，求应用根据该特定信号求出的预测误差滤波器所得的信号(修正估计值系列)在其全样本间独立的反向滤波器。然后，通过将满足规定的反复终止条件时的反向滤波器应用于观测信号而得到复原信号。

在所述的信号失真消除装置中，预测误差滤波器计算单元也可以作为求全帧中各修正估计值的方差的总和最小时的预测误差滤波器，或全帧中各修正估计值的方差的对数值的总和最小时的预测误差滤波器，对各帧的特定信号进行线性预测分析，并输出每个帧的预测误差滤波器；反向滤波器计算单元也可以将所述修正估计值系列在其全样本间独立的所述反向滤波器作为全帧中各修正估计值的归一化峰值的总和最大时的反向滤波器而求，并输出该反向滤波器。

该结构规定相互信息量作为修正系列的独立性的尺度，并通过交替变量法求使该相互信息量最小化的预测误差滤波器和反向滤波器。具体内容将在后面叙述。

或者，在所述的信号失真消除装置中，预测误差滤波器计算单元也可以作为求全帧中各修正估计值的方差的总和最小时的预测误差滤波器，或全帧中各修正估计值的方差的对数值的总和最小时的预测误差滤波器，对各帧的特定信号进行线性预测分析，并输出每个帧的预测误差滤波器；反向滤波器计算单元也可以将所述修正估计值系列在其全样本间独立的所述反向滤波器作为全帧中各修正估计值的方差的总和最小时的反向滤波器或者全帧中各修正估计值的方差的对数值的总和最小时的反向滤波器而求，并输出该反向滤波器。

该结构规定相互信息量作为修正系列的独立性的尺度，并通过交替变量法求使该相互信息量最小化的预测误差滤波器和反向滤波器，但能够不使用信号的高阶统计量，而通过交替变量法求预测误差滤波器和反向滤波器。

所述的信号失真消除装置中可以是如下的装置结构，使预白化处理前置，对由预白化处理所得的白化信号进行所述同样的处理。具体来说，可以是包括以下单元的信号失真消除装置：白化滤波器计算单元，输出对观测信号进行线性预测分析所得的白化滤波器；白化滤波器应用单元，将白化滤波器应用于观测信号而输出白化信号；反向滤波器应用单元，在满足了规定的反复终止条件时，将适用于白化信号的滤波器(以下，称为反向滤波器)应用于白化信号，并将该结果作为复原信号输出，在不满足反复终止条件时，将反向滤波器应用于白化信号，并将该结果作为特定信号输出；预测误差滤波器计算单元，将特定信号进行分帧，输出对各帧的特定信号进行线性预测分析所得的每个帧的预测误差滤波器；反向滤波器计算单元，求对各帧的特定信号拟合了应用与该帧对应的预测误差滤波器所得的各信号(以下，称为修正估计值)的全帧中的修正估计值(以下，称为修正估计值系列)在其全样本间独立的反向滤波器，并输出该反向滤波器；以及控制单元，使反向滤波器应用单元、预测误差滤波器计算单元、反向滤波器计算单元反复执行，直到满足反复终止条件。

为了解决所述课题，本发明的方法是具有以下步骤的信号失真消除方法：反向滤波器应用步骤，在满足了规定的反复终止条件时，反向滤波器应用单元将适用于观测信号的滤波器(以下，称为反向滤波器)应用于观测信号，并将该结果作为复原信号输出，在不满足反复终止条件时，反向滤波器应用单元将反向滤波器应用于观测信号，并将该结果作为特定信号输出；预测误差滤波器计算步骤，预测误差滤波器计算单元将特定信号进行分帧，输出对各帧的特定信号进行线性预测分析所得的每个帧的预测误差滤波器；反向滤波器计算步骤，反向滤波器计算单元，求对于各帧的特定信号拟合了应用与该帧对应的预测误差滤波器所得的各信号(以下，称为修正估计值)的全帧中的修正估计值(以下，称为修正估计值系列)在其全样本间独立的反向滤波器，并输出该反向滤波器；以及控制步骤，控制单元使反向滤波器应用步骤、预测误差滤波器计算步骤、反向滤波器计算步骤反复执行，直到满足所述反复终止条件。

此外，所述的信号失真消除方法中可以是如下的方法，使预白化处理前置，对由预白化处理所得的白化信号进行所述同样的处理。具体来说，是具有以下步骤的信号失真消除方法：白化滤波器计算步骤，白化滤波器计算单元输出对观测信号进行线性预测分析所得的白化滤波器；白化滤波器应用步骤，白化滤波器应用单元将白化滤波器应用于观测信号而输出白化信号；反向滤波器应用步骤，在满足了规定的反复终止条件时，反向滤波器应用单元将适用于白化信号的滤波器(以下，称为反向滤波器)应用于白化信号，并将该结果作为复原信号输出，在不满足所述反复终止条件时，反向滤波器应用单元将反向滤波器应用于白化信号，并将该结果作为特定信号输出；预测误差滤波器计算步骤，预测误差滤波器计算单元将特定信号进行分帧，输出对各帧的特定信号进行线性预测分析所得的每个帧的预测误差滤波器；反向滤波器计算步骤，反向滤波器计算单元，求对于各帧的特定信号拟合了应用与该帧对应的预测误差滤波器所得的各信号(以下，称为修正估计值)的全帧中的修正估计值(以下，称为修正估计值系列)在其全样本间独立的反向滤波器，并输出该反向滤波器；以及控制步骤，控制单元使反向滤波器应用步骤、预测误差滤波器计算步骤、反向滤波器计算步骤反复执行，直到满足反复终止条件。

根据使计算机具有作为本发明的信号失真消除装置功能的信号失真消除程序，能够使计算机作为信号失真消除装置进行动作处理。并且，根据记录了该信号失真消除程序的计算机可读取的程序记录介质，能够使其它计算机具有作为信号失真消除装置的功能，以及使信号失真消除程序流通等。

发明效果

本发明中，不是使用从观测信号求得的预测误差滤波器来降低观测信号中的原信号固有特性的作用，而是使用对观测信号(假设的)应用反向滤波器所得的特定信号(假设的复原信号)所求的预测误差滤波器来降低观测信号中的原信号固有特性。从特定信号求得的预测误差滤波器不易受到传递特性的影响，因此能够更加正确地去除原信号固有的特性。以将这样的预测误差滤波器应用于特定信号所得的信号(修正估计值系列)在全样本间独立所求出的反向滤波器能够高精度地消除传递特性，因此通过将这样的反向滤波器应用于观测信号，能够得到消除了来源于传递特性的失真的精度高的复原信号。

附图说明

图1是表示了用于说明本发明的原理的模型机构的框图。

图2是表示实施方式1的信号失真消除装置(1)的硬件结构例的图。

图3是表示实施方式1的信号失真消除装置(1)的功能结构例的功能方框图。

图4是表示信号失真消除装置(1)的反向滤波器计算部分(13)的功能结构例的功能方框图。

图5是表示实施方式1的信号失真消除处理的流程的处理流程图。

图6是表示实施方式2的信号失真消除装置(1)的功能结构例的功能方框图。

图7是表示实施方式2的信号失真消除处理的流程的处理流程图。

图8是表示将观测信号长度N变为5秒、10秒、20秒、1分、3分时，反复次数R₁和D₅₀值的关系的图。

图9(A)是不包含混响的声音的声谱图，图9(B)是包含混响的声音的声谱图，图9(C)是消除混响后的声音的声谱图。

图10(A)是用于说明消除混响声音的LPC频谱失真的时间变动的图，图10(B)是对应区间的原声音信号的摘录。

图11是表示实施方式3的信号失真消除装置(1)的反向滤波器计算部分(13)的功能结构例的功能方框图。

图12是表示实施方式3的信号失真消除处理的流程的处理流程图。

图13是显示了观测信号长度N设为3秒、4秒、5秒、10秒时的RASTI的值的图。

图14是表示消除混响前后的能量衰减曲线的例子的图。

图15是用于说明以往技术的功能方框图。

具体实施方式

§1本发明的理论

以下，在说明实施方式之前，先说明本发明的理论。

以下的说明中，只要没有特别的声明，设信号源为1了。

1.1信号

作为本发明的对象的信号广泛包含人的声音、音乐、生物信号、用传感器观测测量对象的物理量的电信号等信号。更加优选是，能够作为自回归(Autoregressive：AR)过程来表现的信号，或者表现理想的信号即可。例如语音信号是通常作为阶段平稳的自回归过程所表现的信号，即看作是使独立同一分布信号(i.i.d.：Independent and Identically Distributed)中表示音韵性的AR系统作用的信号(参照参考文献1)。

以下，例举语音信号作为信号的代表例子，说明本发明的理论。

(参考文献1)L.R.Rabiner，R.W.Schafer，“Digital Processing of SpeechSignals”，Bell Laboratories，Incorporated，1978.

1.2语音信号的模型化

首先，将原信号、即语音信号s(t)作为满足以下3个条件的信号进行模型化。

[条件1]语音信号s(t)是由阶段平稳的AR过程生成的。

由该[条件1]，设AR过程的阶数为P、设认为稳定的区间长度为W样本而将语音信号s(t)进行分帧，则第i帧的语音信号s_i(n)表示为如式(1)所示。式(2)表示第i帧的语音信号s_i(n)的样本和分帧前的语音信号s(t)的样本之间的对应。即，第i帧的n号样本在语音信号s(t)中相当于(i-1)W+n号的样本。在式(1)及式(2)中，b_i(k)表示线性预测系数，e_i(n)表示修正。其中，1≤n≤W、1≤t≤N，N是全部样本数。以下，若没有特别的声明，则参数n表示1帧的样本号，参数t表示全部的样本号。此外，设全帧数为F。

$s_i\left(n\right)=\underset{k=1}{\overset{P}{\Sigma }}{b}_i\left(k\right)s_i(n-k)+e_i\left(n\right)---\left(1\right)$

s_i(n)＝s((i-1)W+n)           (2)

再者，对于第i帧中的n号修正e_i(n)，也能够表示第i帧的n号修正e_i(n)和对于分帧前的语音信号s(t)的修正e(t)之间的对应。这种情况下，第i帧的n号修正e_i(n)在修正e(t)中相当于(i-1)W+n号的修正，e_i(n)＝e((i-1)W+n)成立。

对式(1)进行z变换。设左边的z变换为S_i(Z)，右边第二项的z变换为E_i(Z)，B_i(z)＝∑_k＝1^pb_i(k)_z^-k，则右边第一项变成B_i(z)S_i(Z)。因此，式(1)的z变换为，(1-B_i(z))S_i(Z)＝E_i(Z)。再者，z^-1在时域中相当于1抽头(tap)延迟因子。以后，将时域信号(抽头加权系数)用小字符，z域信号(传递函数)用大字符分别表示。1-B_i(z)必须满足最小相位性，要求“1-B_i(z)在复平面的单位圆内有着所有的零点”。

[条件2]属于第i帧的修正系列e_i(1)，...，e_i(W)独立且属于同一分布。修正系列e_i(1)，...，e_i(W)的概率分布的平均及失真度(3阶累计量)为0，峰值(4阶累计量)为正。并且，属于不同帧i、j[i≠j]的修正e_i(n)和e_j(n′)之间也同样是独立的。但是，并不限于它们必须属于同一分布。

[条件3]预测误差滤波器1-B_i(z)在不同的帧之间没有共同的零点。

根据式(1)及式(2)，语音信号s(t)表示为如式(3)所示。[·]是高斯符号。

$s\left(t\right)=\underset{k=1}{\overset{P}{\Sigma }}{b}_i\left(k\right)s(t-k)+e\left(t\right),i=[\frac{t-1}{W}+1]---\left(3\right)$

这时，[条件2]可表达为，“修正过程e(t)是时间上独立的信号。且，它的统计性质(或是统计量)在帧内是恒定的”。此外，[条件3]可表达为，“线性预测系数{b_i(k)}_k＝1^p没有时不变的极点”。

1.3观测信号的模型化

接着，用M个话筒观测语音信号而得到观测信号，并将此时的观测信号进行模型化。其中，M是M≥1的整数。

由m号(1≤m≤M)话筒观测的包含混响的观测信号x_m(t)，利用从声源至m号话筒的路径的传递函数H_m(z)的抽头加权系数{h_m(k)；0≤k≤K；设K为冲击响应的持续时间。}，如式(4)所示进行模型化。这里，例举混响作为语音信号的情况的代表例子，用混响替换传递特性来说明。但是，并不是将传递特性限定为混响。

$x_m\left(t\right)=\underset{k=0}{\overset{K}{\Sigma }}{h}_m\left(k\right)s(t-k)---\left(4\right)$

将M个观测信号统一来表达，则可表示成如式(5)所示。但是，式(5)中，x(t)＝[x₁(t)，...，x_M(t)]^T，h(k)＝[h₁(k)，...，h_M(k)]^T。

$x\left(t\right)=\underset{k=0}{\overset{K}{\Sigma }}h\left(k\right)s(t-k)---\left(5\right)$

1.4消除信号失真的原理

消除信号失真后的复原信号y(t)是利用多通道反向滤波器{G_m(z)；1≤m≤M}的抽头加权系数{g_m(k)；1≤m≤M，0≤k≤L；L是反向滤波器的阶数}，根据式(6)来计算。在本发明中，仅由观测信号x₁(t)，...，x_M(t)估计反向滤波器系数g_m(k)。

$y\left(t\right)=\underset{m=1}{\overset{M}{\Sigma }}\underset{k=0}{\overset{L}{\Sigma }}{g}_m\left(k\right)x_m(t-k)---\left(6\right)$

1.5本发明的基本原理

本发明的基本原理以并列估计传递函数{H_m(z)；1≤m≤M}的反向滤波器{G_m(z)；1≤m≤M}和作为AR系统的滤波器{1/(1-B_i(z))；1≤i≤F}的反向滤波器的预测误差滤波器{1-A_i(z)；1≤i≤F}为主要特征。

为说明该基本原理，图1中表示容纳了所述的模型结构的整个系统的结构框图。根据所述模型化，原信号s(t)可看作对每个帧的修正系列e_i(1)，...，e_i(W)应用AR系统的滤波器1/(1-B_i(z))所得的信号s_i(n)的帧结合，可以说观测信号x(t)是传递函数H(z)对原信号s(t)作用的信号。并且，信号失真消除处理成为对观测信号x(t)作用反向滤波器G(z)而得到复原信号y(t)的处理。这时，将由信号失真消除处理所得的复原信号y(t)进行帧分割，对于各自的帧应用根据各自的信号求得的预测误差滤波器1-A_i(z)所得到的修正估计值d_i(1)，...，d_i(W)与修正系列e_i(1)，...，e_i(W)一致是所期望的。如果预测误差滤波器1-A_i(z)的输出信号d_i(n)满足d_i(n)＝e_i(n)[1≤i≤F，1≤n≤W]，则在[条件3]的条件下将表示成为∑_m＝1^MH_m(z)G_m(z)＝1(关于数学证明，请参照参考文献A)。即，可以说s(t)＝y(t)。这时，1-A_i(z)与1-B_i(z)等价。

(参考文献A)Takuya Yoshioka，Takafumi Hikichi，Masato Miyoshi，Hiroshi G.Okuno：Robust Decomposition of Inverse Filter of Channel and Prediction ErrorFilter of Speech Signal for Dereverberation，Proceedings of the 14th EuropeanSignal Processing Conference (EUSIPCO 2006)，CD-ROM Proceedings，Florence，2006.

但是，实际上，不能将修正e_i(n)[1≤i≤F，1≤n≤W]作为信号失真消除装置的输入信号来利用。图1所示的系统中，从各修正系列e_i(n)得到观测信号x(t)的一系列过程是模型过程，实际上无法知道或很难知道各修正系列e_i(n)、滤波器1/(1-B_i(z))和传递函数H(z)，能够利用的信息只有观测信号x(t)。因此，根据所述[条件2]，估计反向滤波器G_m(z)和预测误差滤波器1-A_i(z)，以将第i帧的修正估计值d_i(1)，...，d_i(W)结合所得的全帧中的修正估计值系列在它的全样本间独立，即修正估计值系列d₁(1)，...，d₁(W)，...，d_i(1)，...，d_i(W)，...d_F(1)，...，d_F(W)独立。

以对比以往的方法的形式来叙述该方法。以往方法是将反向滤波器作为以下问题的解而得到的：“将根据观测信号而求得的预测误差滤波器应用于观测信号，求对于预测误差滤波器应用后的信号应用反向滤波器所得的信号的归一化峰值最大时的反向滤波器”。对此，本发明是将反向滤波器作为以下问题的解而得到的：“对于将反向滤波器应用于观测信号所得的信号，求使应用根据该信号求出的预测误差滤波器所得信号在全样本间为独立时的反向滤波器”。该问题中必须注意的是，预测误差滤波器根据将反向滤波器应用于观测信号所得的信号而求，因此，不仅是反向滤波器，还要一同求预测误差滤波器。

该问题可根据与ICA(Independent Component Analysis)相同的思考方式来公式化。这里是从将相互信息量最小化的观点出发来进行说明，但也可以在例如以使对数似然最大化的方针进行公式化。不管怎样，也只不过是对问题的近似方法的不同。

作为独立性的尺度使用相互信息量(Kullback-Leibler信息量)，则需要解决的问题如式(7)所示被公式化。其中，设g＝[g₁^T，...，g_M^T]^T、g_m＝[g_m(0)，...，g_m(L)]^T、a＝[a₁^T，...，a_F^T]^T、a_i＝[a_i(1)，...，a_i(P)]^T，a_i(k)表示预测误差滤波器的系数。I(U₁，...，U_n)表示概率变元U_i之间的相互信息量。此外，g和a加上符号^的是应当得到的最优解。T表示转置。

$\{\hat{g},\hat{a}\}=\underset{g,a}{\arg \min }I(d_1\left(1\right),...,d_1\left(W\right),...,d_F\left(1\right),...,d_F\left(W\right))---\left(7\right)$

约束条件

[1]‖g‖＝1(其中，‖·‖表示模)

[2]1-A_i(z)在复平面的单位圆内具有所有的零点[1≤i≤F]。

即使修正估计值系列d₁(1)，...，d₁(W)，...，d_i(1)，...，d_i(W)，...d_F(1)，...，d_F(W)的振幅变为常数倍，相互信息量I也不会变化。式(7)的约束条件[1]是用于排除该振幅的不定性的条件。式(7)的约束条件[2]是用于与所述[条件1]对应，将预测误差滤波器限制在最小相位系统的条件。以下，将I看作是以修正估计值系列作为输入，输出它们之间的相互信息量的函数，称为损失函数。

1.6损失函数的推导

为执行式(7)的最优化，必须由有限长度的信号系列{d_i(n)；1≤i≤F，1≤n≤W}估计损失函数I(d₁(1)，...，d_F(W))。若将(多变量)概率变元U的微分熵记为D(U)，则I(d₁(1)，...，d_F(W))由式(8)定义。其中，d＝[d_F^T，...，d₁^T]^T、d_i＝[d_i(W)，...，d_i(1)]^T。

$I(d_1\left(1\right),...,d_F\left(W\right))=\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}D\left(d_i\left(n\right)\right)-D\left(d\right)---\left(8\right)$

设y＝[y_F^T，...，y₁^T]^T、y_i＝[y_i(W)，...，y_i(1)]^T，d使用y，表示为d＝Ay。其中，矩阵A由式(9)和式(10)表示。

[数8]

因此，D(d)表示为式(11)。

D(d)＝D(y)+log det A    (11)

多变量概率变元U的协方差矩阵记为∑(U)，则对于式(11)的右边第二项，∑(d)＝E{dd^T}＝AE{yy^T}A^T＝A∑(y)A^T成立，所以式(12)成立。

$\log \det A=\frac{1}{2}(\log \det \Sigma \left(d\right)-\log \det \Sigma \left(y\right))---\left(12\right)$

将式(11)、式(12)带入式(8)，则得到式(13)。其中，σ(U)²表示概率变元U的方差。

$I(d_1\left(1\right),...,d_1\left(W\right),...,d_F\left(1\right),...,d_F\left(W\right))$

$=\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}D\left(d_i\left(n\right)\right)-\frac{1}{2}\log \det \Sigma \left(d\right)+(\frac{1}{2}\log \det \Sigma \left(y\right)-D\left(y\right))$

$=-\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}(\frac{1}{2}\log \sigma {\left(d_i\left(n\right)\right)}^2-D\left(d_i\left(n\right)\right))$

$+\frac{1}{2}(\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}\log \sigma {\left(d_i\left(n\right)\right)}^2-\log \det \Sigma \left(d\right))$

$+(\frac{1}{2}\log \det \Sigma \left(y\right)-D\left(y\right))$

$=-\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}J\left(d_i\left(n\right)\right)+C(d_1\left(1\right),...,d_F\left(W\right))+J\left(y\right)---\left(13\right)$

式(13)中J(U)是(多变量)概率变元U的负熵(negentropy)。负熵取表示U的非高斯性程度的非负值，只有在U为高斯分布的情况下取0。C(U₁，...，U_n)由式(14)定义。C(U₁，...，U_n)取表示概率变元U_i之间的相关程度的非负值，只有在这些概率变元之间没有相关的情况下取0。

$C(U_1,...,U_n)=\frac{1}{2}(\underset{i=1}{\overset{n}{\Sigma }}\log \sigma {\left(U_i\right)}^2-\log \det \Sigma \left({[U_1,...,U_n]}^T\right))---\left(14\right)$

但是，若设s＝[s_F^T，...，s₁^T]^T、s_i＝[s_i(W)，...，s_i(1)]^T，则J(y)＝J(s)＝constant(证明略)，因此，式(13)可以进一步简化如式(15)所示。

$I(d_1\left(1\right),...,d_F\left(W\right))$

$=-\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}J\left(d_i\left(n\right)\right)+C(d_1\left(1\right),...,d_F\left(W\right))+\mathrm{const}---\left(15\right)$

只要根据以上求式(16)的最优化问题即可。

$\{\hat{g},\hat{a}\}=\underset{g,a}{\arg m\mathrm{in}}(-\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}J\left(d_i\left(n\right)\right)+C(d_1\left(1\right),...,d_F\left(W\right)))---\left(16\right)$

[约束条件]

[1]‖g‖＝1(其中，‖·‖表示模)

[2]1-A_i(z)在复平面的单位圆内具有所有的零点[1≤i≤F]。

1.7根据交替变量法的最优化

对于式(16)，根据交替变量法将g和a进行最优化。即，若将第r次反复中的g和a的估计值分别表示为g^^(r)、a^^(r)，则根据式(17)和式(18)的交互的最优化得到被更新的估计值g^^(r+1)、a^^(r+1)。再者，g^和a^表示符号^加在了g、a的各自的上面。例如，若设反复次数的上限为R₁，则第R₁次得到的g^^(R1+1)、a^^(R1+1)为式(16)的最优解。上标字符R1为R₁。

${\hat{a}}^{(r+1)}=\underset{a}{\arg \min }(-\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}J\left(d_i\left(n\right)\right)+C(d_1\left(1\right),...,d_F\left(W\right)))---\left(17\right)$

[约束条件]

[1]g＝g^^(r)

[2]1-A_i(z)在复平面的单位圆内有着所有的零点[1≤i≤F]。

${\hat{g}}^{(r+1)}=\underset{g}{\arg \min }(-\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}J\left(d_i\left(n\right)\right)+C(d_1\left(1\right),...,d_F\left(W\right)))---\left(18\right)$

[约束条件]

[1]a＝a^^(r+1)

[2]‖g‖＝1

式(17)的目的在于，根据用于消除传递特性的反向滤波器的当前的估计值来估计用于消除原信号的固有特性的预测误差滤波器。同样，式(18)的目的在于，根据预测误差滤波器的当前的估计值来估计反向滤波器。通过反复这两类的最优化，以使修正估计值系列d₁(1)，...，d₁(W)，...，d_i(1)，...，d_i(W)，...d_F(1)，...，d_F(W)相互独立，从而能够将反向滤波器和预测误差滤波器并行估计。因此，这里的反复对于反向滤波器的高精度的估计很重要。但是，由图8明显看出，随着处理的观测信号长度变长，即使反复次数为1次也能得到一定程度的效果。因此，本发明中，反复次数可以是1次。

1.8a的最优化

本发明中如下进行式(17)的最优化。

首先需要注意的是，相对于C(d₁(1)，...，d_F(W))是与d_i(n)的二阶统计量关联的值，J(d_i(n))是与d_i(n)的高阶统计量关联的值。二阶统计量只提供信号的振幅信息，但高阶统计量还提供相位信息。因此，一般，包含高阶统计量的最优化有推导非最小相位系统的可能性。所以，由1-A_i(z)为最小相位系统的约束条件，在a的最优化中解决式(19)的最优化问题。

${\hat{a}}^{(r+1)}=\underset{a}{\arg \min }C(d_1\left(1\right),...,d_F\left(W\right))---\left(19\right)$

约束条件

[1]g＝g^^(r)

[2]1-A_i(z)在复平面的单位圆内具有全部的零点[1≤i≤F]。

C(d₁(1)，...，d_F(W))由式(20)给出。

$C(d_1\left(1\right),...,d_F\left(W\right))=\frac{1}{2}(\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}\log \sigma {\left(d_i\left(n\right)\right)}^2-\log \det \Sigma \left(d\right))---\left(20\right)$

这里，矩阵A如式(9)和式(10)所示是上三角矩阵，其对角分量全部为1，因此，log detA＝0。通过将该式代入式(12)，得到式(21)的关系。

log det∑(d)＝log det∑(y)＝constant    (21)

因此，式(19)与式(22)的最优化问题是等价的。再者，必须注意到式(22)是反映了所述[条件2]的表现。要说明式(22)的话，可以说式(22)是求将第i帧中的修正估计值d_i(1)，...，d_i(W)的方差的对数值在全帧中相加的值最小时的a。

${\hat{a}}^{(r+1)}=\underset{a}{\arg \min }\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}\log \sigma {\left(d_i\left(n\right)\right)}^2---\left(22\right)$

约束条件

[1]g＝g^^(r)

[2]1-A_i(z)在复平面的单位圆内有着全部的零点[1≤i≤F]。

解决由式(22)表示的最优化问题与对于将g^^(r)给出的反向滤波器应用于观测信号所得的特定信号在各帧中进行线性预测分析是等价的，一定能够得到最小相位预测误差滤波器。关于线性预测分析，请参照所述参考文献1。

再者，在式(22)中，将第i帧的修正估计值d_i(1)，...，d_i(W)的方差的对数值在全帧中的总和最小时的a作为a^^(r+1)来求，但不是限定于此。所述各式中没有清楚写明对数函数的底数(base)，一般来说将底数设为10或者纳披尔数是惯例，不管怎样，底数都要比1大。这时，对数函数是单调递增函数，所以能够将第i帧的修正估计值d_i(1)，...，d_i(W)的方差在全帧中的总和最小时的a作为a^^(r+1)来求。

1.9g的最优化

本发明中，如下进行式(18)的最优化。

如上所述，C(d₁(1)，...，d_F(W))是有关{d_i(n)；1≤i≤F，1≤n≤W}的相关程度的指标，但在第(r+1)次的a的最优化中被最小化，因此与∑_i＝1^F∑_n＝1^WJ(d_i(n))相比，可以忽略。所以在g的最优化中，解决式(23)的最优化问题。

${\hat{g}}^{(r+1)}=\underset{g}{\arg \min }(-\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}J\left(d_i\left(n\right)\right))---\left(23\right)$

约束条件

[1]a＝a^^(r+1)

[2]‖g‖＝1

J(d_i(n))可根据[条件2]，由式(24)进行近似。该详细内容请参照参考文献2。但是，对于概率变元U，κ₄(U)表示U的峰值(4阶累计量)。式(24)的右边称为第i帧中的归一化峰值。

(参考文献2)A.Hyvarinen，J.Karhunen，E.Oja，″INDEPENDENTCOMPONENT ANALYSIS″，John Wiley & Sons，Inc.2001.

$J\left(d_i\left(n\right)\right)\approx \frac{{\kappa }_4{\left(d_i\left(n\right)\right)}^2}{\sigma {\left(d_i\left(n\right)\right)}^8}---\left(24\right)$

根据[条件2]，语音信号的修正的峰值为正，所以κ₄(d_i(n))/σ(d_i(n))⁴为正。因此，式(23)的最优化问题归结为式(25)的最优化问题。σ(d_i(n)，κ₄(d_i(n))是根据[条件1]中说明的语音信号的局部恒定性，由各帧内的样本进行计算。式(26)中，添加了1/W，不过这只是为了后续计算的方便，并不会对由式(25)求g的最优解产生影响。由式(25)和式(26)，归一化峰值在全帧中的总和最大时的g成为g^^(r+1)。再者，必须注意到式(25)和式(26)是反映了[条件2]的表现。若说明式(25)和式(26)，则可以说它们求将第i帧中的归一化峰值在全帧中相加的值最大时的g。

${\hat{g}}^{(r+1)}=\underset{g}{\arg \min }Q---\left(25\right)$

$Q=\frac{1}{W}\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}\frac{{\kappa }_4\left(d_i\left(n\right)\right)}{\sigma {\left(d_i\left(n\right)\right)}^4}---\left(26\right)$

约束条件

[1]a＝a^^(r+1)

[2]‖g‖＝1

根据式(25)求g的最优解，只要求用g将Q进行微分的式子等于零时的值即可。一般来讲，该解根据式(27)表示的更新规则而求。g′除以g′的模是为了施加所述约束条件[2]。η(u)表示学习率。u表示g的最优化中的更新次数。

$g^{\prime }\leftarrow g^{<u>}+\eta \left(u\right)·▿Q_{g=g^{<u>}}$

$g^{<u+1>}\leftarrow \frac{g^{\prime }}{\left|\right|g^{\prime }\left|\right|}---\left(27\right)$

式(27)中，由式(28)和式(29)给出。

$▿Q_g=[\frac{\partial Q}{\partial g_1\left(0\right)},...,\frac{\partial Q}{\partial g_1\left(L\right)},...,\frac{\partial Q}{\partial g_M\left(0\right)},...,\frac{\partial Q}{\partial g_M\left(L\right)}]---\left(28\right)$

$\frac{\partial Q}{\partial g_m\left(k\right)}=\underset{i=1}{\overset{F}{\Sigma }}\frac{4R}{E{\left\{d_i{\left(n\right)}^2\right\}}^4}$

R＝E{d_i(n)³v_mi(n-k)}E{d_i(n)²}²                 (29)

-E{d_i(n)⁴}E{d_i(n)²}E{d_i(n)v_mi(n-k)}

式(29)中，d_i(n)由式(30)给出，v_mi(n)是由式(31)和式(32)给出。x_mi(n)是由m号的话筒观测的观测信号的i号帧的信号。

$d_i\left(n\right)=y_i\left(n\right)-\underset{k=1}{\overset{P}{\Sigma }}{a}_i\left(k\right)y_i(n-k)---\left(30\right)$

$v_{\mathrm{mi}}\left(n\right)=x_{\mathrm{mi}}\left(n\right)-\underset{k=1}{\overset{P}{\Sigma }}{a}_i\left(k\right)x_{\mathrm{mi}}(n-k)---\left(31\right)$

x_mi(n)＝x_m((i-1)W+n)                            (32)

§2基于二阶统计量的信号失真消除

所述的以往方法的信号失真消除方法需要较长时间的观测信号(例如是20秒左右。)。这是因为一般计算像归一化峰值这样的高阶统计量需要大量的观测信号的样本。但是，实际上，很多情况下是无法利用这样长时间的观测信号的。因此，以往方法的信号失真消除方法，其应用领域是极为有限的。

再者，高阶统计量的计算是较为复杂的，因此，用以往方法的信号失真消除方法，装置的结构容易变得复杂。

因此，说明在观测信号为较短时间(例如是3秒到5秒程度。)的情况下也有效，且计算比以往容易的信号失真消除原理。该原理只利用信号的二阶统计量，由§1中说明的本发明的基本原理所派生。

2.1基于二阶统计量的信号失真消除原理

基于二阶统计量的信号失真消除中，在所述3个条件中，设定以下两个条件。

[条件4]M≥2。即，使用多个话筒。

[条件5]H_m＝{h_m(k)}_k＝0^k在不同的话筒m之间没有共同的零点。

在所述式(16)的最优化问题中求出使包含有关高阶统计量的值的负熵J和表示概率变元之间的相关程度的指标C的值最小化的a和g。

表示概率变元之间的相关程度的指标C由二阶统计量所规定。因此，将需要解决的最优化问题由式(33)公式化。

$\{\hat{g},\hat{a}\}=\arg \min C(d_1\left(1\right),...,d_F\left(W\right))$

$=\underset{g,a}{\overset{g,a}{\arg \min }}(\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}\log \sigma {\left(d_i\left(n\right)\right)}^2-\log \det \Sigma \left(d\right))---\left(33\right)$

参考式(21)，则式(33)的最优化问题将转化为式(34)的最优化问题。再者，需要注意式(34)是反映了所述[条件2]的表现。若说明式(34)，则可以说式(34)是求将第i帧中的修正估计值d_i(1)，...，d_i(W)的方差的对数值在全帧中相加的值最小时的g和a。

$\{\hat{g},\hat{a}\}=\underset{g,a}{\arg \min }\left(\underset{i=1}{\overset{F}{\Sigma }}\underset{n=1}{\overset{W}{\Sigma }}\log \sigma {\left(d_i\left(n\right)\right)}^2\right)---\left(34\right)$

可是，所述的[条件4]和[条件5]成立时，多通道的观测信号可以作为根据来自声源的原信号而被驱动的AR系统来捕捉(参照参考文献3)。这意味着能够将反向滤波器G的最前头抽头如式(35)进行固定。但是，与m＝1对应的话筒是最接近声源的话筒。

(参考文献3)K.Aded-Meraim，E.Moulines，and P.Loubaton.Prediction errormethod for second-order blind identification.IEEE Trans.Signal Processing，Vol.45，No.3，pp.694-705，1997.

$g_m\left(0\right)=\left(\begin{array}{cc}1& m=1\\ 0& m=2,...,M\end{array}\right)---\left(35\right)$

将由式(34)和式(35)规定的g作为系数的反向滤波器G按照式(6)应用于观测信号x(t)，从而得到传递特性被消除的复原信号y(t)。

2.2a的最优化

对于式(34)，根据交替变量的方法将g和a进行最优化。

在固定了反向滤波器的系数g_m(k)的状态下，关于预测误差滤波器的系数a_i(k)，将式(34)的损失函数最小化。

这时，注意以下两点。第一点在于g＝[g₁^T，...，g_M^T]^T是被固定的，因此作为反向滤波器G的输出的复原信号y(t)在预测误差滤波器的最优化中是不变的。第二点在于，i号帧中的预测误差滤波器的系数a_i(1)，...，a_i(P)是仅对d_i(1)，...，d_i(W)有用。

因此，在各个帧中，只要估计预测误差滤波器的系数a_i(1)，...，a_i(P)以使∑_n＝1^Wlogσ(d_i(n))²最小化即可。由[条件2]，修正估计值d_i(1)，...，d_i(W)的方差在帧内是恒定的，因此∑_n＝1^Wlogσ(d_i(n))²的最小化与W*σ(d_i(n))²的最小化等价。*符号表示乘法。方差σ(d_i(n))²可作为<d_i(n)²>_n＝1^W计算。但是，<d_i(n)²>_n＝1^W表示使用相当于1帧的修正估计值d_i(1)，...，d_i(W)所计算的d_i(n)的2次方平均。结果，系数a_i(k)作为W*<d_i(n)²>_n＝1^W、即d_i(n)的2次方总和最小时的a_i(k)而被估计。这样的系数a_i(k)根据线性预测分析的方法而计算。

再者，这里是对将第i帧中的修正估计值d_i(1)，...，d_i(W)的方差的对数值在全帧中的总和最小时的a作为a^^(r+1)来求进行了说明，但并不是限定于此。如上所述，可将第i帧中的修正估计值d_i(1)，...，d_i(W)的方差在全帧中的总和最小时的a作为a^^(r+1)来求。

2.3g的最优化

在固定了预测误差滤波器的系数a_i(k)的状态下，关于反向滤波器的系数g_m(k)，将式(34)的损失函数最小化。

与反向滤波器的系数g_m(k)有关的损失函数的最小化中使用梯度法。若使用[条件2]，则式(34)的最优化问题被转化为式(36)的最优化问题。

$\hat{g}=\underset{g}{\arg \min }\left(\underset{i=1}{\overset{F}{\Sigma }}\log \left({<d_i{\left(n\right)}^2>}_{n=1}^W\right)\right)---\left(36\right)$

就根据式(36)求g的最优解来说，求将∑_i＝1^Flog<d_i(n)²>_n＝1^W用g微分的式子等于0时的解即可。该解一般可根据式(37)表示的更新规则而求。δ表示学习率。1≤m≤M、1≤k≤L。再者，在式(37)中不得不注意到，由式(35)的条件，‖g‖＝1的约束条件没有被施加。同样，由式(35)的条件，k的取值范围是1≤k≤L。

$g_m{\left(k\right)}^{\prime }=g_m\left(k\right)+\delta \underset{i=1}{\overset{F}{\Sigma }}\frac{{<d_i\left(n\right)v_{\mathrm{mi}}(n-k)>}_{n=1}^W}{{<d_i{\left(n\right)}^2>}_{n=1}^W}---\left(37\right)$

$v_{\mathrm{mi}}\left(n\right)=x_{\mathrm{mi}}\left(n\right)-\underset{k=1}{\overset{P}{\Sigma }}{a}_i\left(k\right)x_{\mathrm{mi}}(n-k)---\left(38\right)$

与所述式(29)或者所述非专利文献1所记载的式(3)相比，很明显，式(37)的右边第二项表示为二阶统计量，该计算中不需要高阶统计量。因此，计算高阶统计量在短时间的观测信号的情况下也同样有效，计算本身也容易。

再者，式(36)中，将第i帧的修正估计值d_i(1)，...，d_i(W)的方差的对数值在全帧中的总和最小时的g作为g^来求，但不是限定于此。所述各式中没有清楚写明对数函数的底数(base)，一般来说将底数设为10或者纳披尔数是惯例，不管怎样，底数都要比1大。这时，对数函数是单调递增函数，所以能够将第i帧的修正估计值d_i(1)，...，d_i(W)的方差在全帧中的总和最小时的g作为g^来求。这种情况下，不能应用式(37)所示的根新规则，重新求将∑_i＝1^F<d_i(n)²>_n＝1^W用g微分的式子等于0时的解即可。由于该结果得到的更新规则可用与ICA同样的思考方式被公式化，所以省略。

§3预白化

本发明的信号失真消除中可以应用预白化。通过将观测信号进行预白化，最优化计算的稳定化，特别是更新规则的高速收敛变得可能。

通过X阶的线性预测分析来计算用于白化各话筒得到的观测信号系列全体{x_m(t)；1≤t≤N}的滤波器(白化滤波器)的系数{f_m(k)；0≤k≤X}。

根据式(39)，将所述白化滤波器应用于各话筒所得的观测信号x_m(t)。w_m(t)表示将m号话筒所得的观测信号x_m(t)白化后的信号。

$w_m\left(t\right)=\underset{k=0}{\overset{X}{\Sigma }}{f}_m\left(k\right)x_m(t-k)---\left(39\right)$

这时，式(31)和式(38)变更为式(40)、式(32)变更为式(41)即可。

$v_{\mathrm{mi}}\left(n\right)=w_{\mathrm{mi}}\left(n\right)-\underset{k=1}{\overset{P}{\Sigma }}{a}_i\left(k\right)w_{\mathrm{mi}}(n-k)---\left(40\right)$

w_mi(n)＝w_m((i-1)W+n)           (41)

§4实施方式

以下，参照附图说明本发明的实施方式。作为本发明的实施方式，并不限定于后述的各实施方式，只要是实现各节中说明的原理的实施方式既可。

<实施方式1>

实施本发明的实施方式1时，按照以下步骤处理由传感器得到的观测信号。这里，以具体说明实施方式的观点出发，将语音信号作为信号例举进行说明。

再者，在说明实施方式1之前，先概括说明观测信号以及分帧处理。

((观测信号))

通过未图示的传感器(例如，话筒)得到的模拟信号(该模拟信号中叠加着来源于传递特性的失真。)，例如以8000Hz的采样频率而被采样，变换成被适当量化的离散信号。以下，将该离散信号设为观测信号。执行从模拟信号到观测信号的A/D变换等所需要的构成元素(方式)，不论是哪一种都根据公知技术的常用方式来完成，因此省略说明及图示。

((分帧处理))

未图示的信号截取单元从离散信号在时间轴方向以一定时间幅度移动起点，并将规定时间长度的离散信号作为帧进行截取。例如将200样本点(8000Hz×25ms)长度的离散信号每隔80样本点(8000Hz×10ms)一边移动起点一边进行截取。各帧的截取，只要将公知的窗函数(例如，哈明窗、高斯窗、矩形窗等)应用于离散信号即可。基于应用窗函数的帧截取根据公知技术的常用方式来完成。

说明用计算机(通用机)执行本发明的实施方式1的信号失真消除装置(1)时的硬件构成例子。

如图2所示，信号失真消除装置(1)包括：可连接键盘、指示器等的输入部分(11)；可连接液晶显示器、CRT(Cathode Ray Tube)显示器等的输出部分(12)；可连接能够与信号失真消除装置(1)外部的通信装置(例如，通信电缆、LAN卡、路由器、调制解调器等)进行通信的通信部分(13)；DSP(Digital Signal Processor)(14)[也可以是CPU(Central Processing Unit)。此外，也可以包括高速缓冲存储器和寄存器(19)等。]；作为存储器的RAM(15)、ROM(16)；作为硬盘、光盘、半导体存储器等的外部存储装置(17)；以及用于能够在这些输入部分(11)、输出部分(12)、通信部分(13)、DSP(14)、RAM(15)、ROM(16)、外部存储装置(17)之间进行数据的交换而连接的总线(18)。此外根据需要，也可以在信号失真消除装置(1)中设置可读取CD-ROM(Compact Disc Read Only Memory)、DVD(Digital VersatileDisc)等存储介质的装置(设备)等。

信号失真消除装置(1)的外部存储装置(17)中存储着用于消除信号失真的程序以及该程序的处理中所需要的数据(观测信号)等[不限定于外部存储装置，例如也可以是将程序存储到读出专用存储装置、即ROM。]。此外，通过这些程序的处理所得到的数据等被适当存储在RAM和外部存储装置等，在提供给其他程序的处理时，从RAM和外部存储装置等读入。

更具体地说，信号失真消除装置(1)的外部存储装置(17)[或者是ROM等]中存储着用于将反向滤波器应用于观测信号的处理的程序；从将反向滤波器应用于观测信号所得的信号求预测误差滤波器的处理的程序；用于从预测误差滤波器求反向滤波器的处理的程序；以及这些程序的处理中所需要的数据(帧单位的观测信号等)。此外还适当保存用于控制基于这些程序的处理的控制程序。

在实施方式1的信号失真消除装置(1)中，存储在外部存储装置(17)[或者是ROM等]的各程序和这些程序的处理所需要的数据根据需要而被读入RAM(15)，在DSP(14)被解释执行/处理。其结果，通过DSP(14)实现规定的功能(反向滤波器应用部分、预测误差滤波器计算部分、反向滤波器计算部分、控制部分)，从而实现信号失真消除。

下面参照图3～图5顺序说明信号失真消除装置(1)中的信号失真消除处理的流程。

大致的处理步骤为：(a)求对观测信号x(t)应用反向滤波器的信号(以下，称为特定信号)；(b)从特定信号求预测误差滤波器；(c)从该预测误差滤波器求反向滤波器；(d)反复所述(a)、(b)、(c)的处理而求最优的反向滤波器；(e)将最优化的反向滤波器对观测信号应用的信号作为复原信号y(t)来得到。

(b)相当于所述的a的最优化，(c)相当于所述的g的最优化，(d)相当于式(17)及式(18)。(d)的处理的反复次数设为预先决定的次数R₁。即，设为1≤r≤R₁。此外，(c)的处理中，最优化g的更新规则的更新次数设为预先决定的次数R₂。即，设为1≤u≤R₂。(d)的处理，即，在每进行一次(a)、(b)、(c)的一连串的处理，就进行根据更新规则的R₂次的更新。实施方式中，次数R₁设为预先决定的次数，但并不限定于此，例如也可以在计算出r次的g时式(26)的Q值和计算出r+1次的g时式(26)的Q值的差的绝对值成为规定的正的微小值ε以下(或者低于ε)时终止反复。同样的，次数R₂设为预先决定的次数，但并不限定于此，例如也可以在计算出u次的g时式(26)的Q值和计算出u+1次的g时式(26)的Q值的差的绝对值成为规定的正的微小值ε以下(或者不到ε)时终止反复。

(步骤S100)

反向滤波器应用部分(14)根据式(42)将反向滤波器应用于观测信号x(t)＝[x₁(t)，...，x_m(t)，...，x_M(t)]^T，从而求特定信号y(t)。特定信号y(t)在计算上与复原信号完全相同，但这里为了明示并非是经过后述的R₁次的处理而求得的复原信号，因而称之为特定信号。这里，t表示所有的样本号，设为1≤t≤N。N是全样本数。实施方式1中，话筒的个数M可以是1以上。

$y\left(t\right)=\underset{m=1}{\overset{M}{\Sigma }}\underset{k=0}{\overset{L}{\Sigma }}{g}_m\left(k\right)x_m(t-k)---\left(42\right)$

作为反向滤波器的系数列{g_m(k)；0≤k≤L}，在反复次数R₁的第一次中，使用预先决定的初始值，而在第二次以后则使用由后述的反向滤波器计算部分(13)求得的反向滤波器g^^(r+1)。

(步骤S101)

预测误差滤波器计算部分(15)由进行分帧处理的分帧处理部分(151)和帧预测误差滤波器计算部分(152)构成。再者，帧预测误差滤波器计算部分(152)由从第i帧的特定信号求预测误差滤波器的第i帧用预测误差滤波器计算部分(152i)。但是，i是满足1≤i≤F的整数。

分帧处理部分(151)将反向滤波器应用部分(14)求得的特定信号{y(t)；1≤t≤N}进行分帧处理。例如式(43)所示，分帧处理通过每W点移动并应用截取W点的窗函数来进行。{y_i(n)；1≤n≤W}表示i号帧所包含的特定信号列。

y_i(n)＝y((i-1)W+n)             (43)

再者，第i帧用预测误差滤波器计算部分(152i)按照式(22)对第i帧的特定信号列{y_i(n)；1≤n≤W}进行P阶的线性预测分析，计算预测误差滤波器的系数列{a_i(k)；1≤k≤P}。该计算方法请参照所述参考文献1。这里得到的a₁(1)，...，a₁(P)，...，a_i(1)，...，a_i(P)，...，a_F(1)，...，a_F(P)供给式(22)的a^^(r+1)。

(步骤S102)

参照图4说明反向滤波器计算部分(13)的功能结构例。反向滤波器计算部分(13)由梯度计算部分(131)、反向滤波器更新部分(132)以及更新用反向滤波器应用部分(133)构成。并且，梯度计算部分(131)构成为包括：具有作为对观测信号的预测误差滤波器应用部分功能的第1预测误差滤波器应用部分(1311)；具有作为对观测信号应用更新用反向滤波器所得的信号(更新用反向滤波器应用后的信号)的预测误差滤波器应用部分功能的第2预测误差滤波器应用部分(1312)；以及梯度矢量计算部分(1313)。这里，更新用反向滤波器相当于式(27)的g^<u>。

第1预测误差滤波器应用部分(1311)将由m号[1≤m≤M]话筒所观测的观测信号x_m(t)进行分帧，对于各帧，将步骤S101的处理所得的i号预测误差滤波器a_i(k)应用于i号帧的信号x_mi(n)而计算预测误差滤波器应用后的信号v_mi(n)(参照式(31))。这里叙述的处理的详细例子将在后述的实施方式3中说明。

第2预测误差滤波器应用部分(1312)将更新用反向滤波器应用后的信号y(t)进行分帧，对于各帧，将步骤S101的处理所得的i号预测误差滤波器a_i(k)应用于i号帧的信号y_i(n)而计算预测误差滤波器应用后的修正估计值d_i(1)，...，d_i(W)(参照式(30))。再者，更新用反向滤波器应用后的信号y(t)的初始值设为由步骤S100的处理所得的信号即可。然后，第2预测误差滤波器应用部分(1312)输入后述的更新用反向滤波器应用部分(133)输出的更新用反向滤波器应用后的信号y(t)。这里叙述的处理的详细例子将在后述的实施方式3中说明。

梯度矢量计算部分(1313)使用信号v_mi(n)和修正估计值d_i(n)，计算当前的更新用反向滤波器g^<u>的梯度矢量(参照式(28)及式(29))。使用有限个样本v_mi(n)和d_i(n)运算式(29)时，由样本求期望值E即可。这里叙述的处理的详细例子将在后述的实施方式3中说明。

反向滤波器更新部分(132)以当前的更新用反向滤波器g^<u>、学习率η(u)、梯度矢量作为输入，按照式(27)求第u+1次的更新用反向滤波器g^<u+1>。式(27)意味着将求得的g^<u+1>看作新的g^<u>来进行更新。

更新用反向滤波器应用部分(133)以通过反向滤波器更新部分(132)所得到的g^<u+1>，即新的g^<u>和观测信号x(t)作为输入，按照式(42)求更新用反向滤波器应用后的信号y(t)。具体来讲，作为式(42)的g_m(k)，使用由第u+1次的更新所得到的g进行计算。该计算所得到的更新用反向滤波器应用后的信号y(t)成为第2预测误差滤波器应用部分(1312)的输入。再者，更新用反向滤波器应用后的信号y(t)计算上与复原信号是完全相同的，但这里为了明示并非是经过后述的R₁次的处理而求得的复原信号，而是为了进行更新规则而计算的信号，因而称之为更新用反向滤波器应用后的信号。

通过控制部分(600)的控制进行R₂次的更新，作为该结果所得到的g^<R2+1>相当于式(25)的g^^(r+1)。上标字符R2是R₂。反向滤波器计算部分(13)输出g^^(r+1)。

根据控制部分(500)的控制，每进行一次所述的一连串的处理，r就加1，直到r等于R₁，即通过反复R₁次所述的一连串的处理(步骤S103)而得到g^^(R1+1)。上标字符R1是R₁。该g^^(R1+1)被设为式(16)的最优解。因此，在得到g^^(R1+1)的阶段，反向滤波器应用部分(14)通过按照式(42)将反向滤波器g^^(R1+1)应用于观测信号x(t)＝[x₁(t)，...，x_M(t)]^T，能够得到复原信号y(t)(步骤S104)。

<实施方式2>

实施方式2相当于实施方式1的变形例。具体来讲，是进行§3中叙述的预白化的方式。因此，对于与实施方式1不同的部分，参照图6和图7附加说明。再者，预白化是对观测信号进行的预处理，所以这里说明的进行预白化的方式也可以适用于后述的实施方式3。

实施方式2中，信号失真消除装置(1)的外部存储装置(17)[或者ROM等]中还存储着用于求白化滤波器的处理的程序，以及用于进行将白化滤波器应用于观测信号的处理的程序。

实施方式2的信号失真消除装置(1)中，存储在外部存储装置(17)[或者ROM等]的各程序和该各程序的处理所需要的数据根据需要而被读入RAM(15)，被DSP(14)解释执行/处理。其结果，通过DSP(14)实现规定的功能(反向滤波器应用部分、预测误差滤波器计算部分、反向滤波器计算部分、白化滤波器计算部分、白化滤波器应用部分)，实现消除信号失真。

(步骤S100a)

白化滤波器计算部分(11)根据X阶的线性预测分析来计算用于将各话筒所得到的观测信号全体{x_m(t)；1≤t≤N}进行白化的滤波器(白化滤波器)的系数{f_m(k)；0≤k≤X}。该计算与线性预测分析相同，请参考所述参考文献1。白化滤波器的系数成为白化滤波器应用部分(12)的输入。

(步骤S100b)

白化滤波器应用部分(12)按照式(39)，将所述白化滤波器应用于各话筒所得到的观测信号而得到白化信号w_m(t)。如上所述，式(31)变更为式(40)即可，所以在实施方式1中，将根据反向滤波器计算部分(13)，特别是第1预测误差滤波器应用部分(1311)的处理改为根据式(40)的计算处理而不是根据式(31)的计算处理即可。此外，实施方式1中，将根据反向滤波器应用部分(14)的处理改为根据式(44)的计算处理而不是根据式(42)的计算处理即可。步骤S100b的处理之后，进行实施方式1的步骤S100～S104的处理，但这些处理中，将实施方式1的各处理中的观测信号替换成由步骤S100b的处理所得到的白化信号而进行与实施方式1相同的处理。在此意义下，在图7中，在表示相当于与实施方式1的步骤S100～S104的各处理的处理的符号上附加了记号′。

$y\left(t\right)=\underset{m=1}{\overset{M}{\Sigma }}\underset{k=0}{\overset{L}{\Sigma }}{g}_m\left(k\right)w_m(t-k)---\left(44\right)$

<实施例1>

发明人进行了实施方式2的验证实验，因此表示该实验结果。作为实验条件，设话筒数量M＝4、白化滤波器的阶数X＝500、反向滤波器的阶数L＝1000、窗函数的截取样本数(1帧的样本数)W＝200、预测误差滤波器的阶数P＝16、反复次数R₁＝10、反向滤波器计算部分的更新次数R₂＝20。学习率η(u)，设其初始值为0.05，如果根据式(27)而式(26)的值减少，则通过将η(u)的值依次减半，使得式(26)的值必然增大。输入到图6所示的反向滤波器应用部分(14)的初始反向滤波器设定为如式(45)所示。

作为消除信号失真的指标，使用D₅₀值(直到对于脉冲响应的全能量的初期50msec的能量比)评价本发明的实施方式2的效果。从连续发声数据库取出男女各一名的发声，通过叠加在混响时间为0.5秒的混响室测量的脉冲响应而合成观测信号。

图8表示对于男声及女声将观测信号长度N变为5秒、10秒、20秒、1分、3分时，反复次数R₁(使具有各信号长度N的观测信号执行图6所示的反向滤波器应用部分(14)、预测误差滤波器计算部分(15)、反向滤波器计算部分(13)的一轮处理而求反向滤波器的次数)和D₅₀值的关系。无论是哪一种情况，增加反复次数则D₅₀值提高，明显看出反复处理的效果。特别是即使在观测信号长度为5～10秒左右的较短的长度时，也能看出根据反复处理而使得D₅₀值大大提高。

此外，通过比较声音频谱来验证了本发明的实施方式2的效果。

图9A表示观测信号长度为1分钟时不包含混响的声音(原声音)的声谱图的摘录，图9B表示观测信号长度为1分钟时包含混响的声音(观测声音)的声谱图的摘录，图9C表示观测信号长度为1分钟时消除混响后的声音(复原信号)的声谱图的摘录。由图9A和图9C的对比以及图9B和图9C的对比可知，观测信号所包含的混响被抑制，作为原声音固有的特征的谐波构造和共振峰构造被恢复。

此外，使用LPC频谱失真验证了本发明的实施方式2的效果。

图10B表示原声音的波形，图10A表示原声音和观测声音的LPC频谱失真的时间序列(图中的虚线)以及原声音和复原声音的LPC频谱失真的时间序列(图中的实线)。图10A和图10B的各横轴表示单位为秒的时间，两图中使时标一致。图10B的纵轴表示振幅值。但是，只要知道原信号的相对大小即可，由该意旨而在图10B的纵轴中没有明示其单位。图10A的纵轴表示LPC频谱失真SD(dB)。

由图10A可知，原声音和复原声音的LPC频谱失真的时间序列(图中的实线)通常比原声音和观测声音的LPC频谱失真的时间序列(图中的虚线)还要小。再者，相对于观测信号中LPC频谱失真的平均为5.39dB、方差为4.20dB，而复原声音中平均为2.38dB、方差为2.00dB。

此外，由图10A与图10B的对比可知，在原声音和复原声音的LPC频谱失真的时间序列(图中的实线)表示较大值的区间(例如参照约1.0秒～约1.2秒的区间)，原声音的波形的振幅值大致为0。实际上，该区间中没有发声，是无声区间。因此，实际上所感觉的失真相当小。即，可得到以下结论：发声区间中的原声音和复原声音的LPC频谱失真的时间序列(图中的实线)比原声音和观测声音的LPC频谱失真的时间序列(图中的虚线)还要相当小，因此，能够以高精度复原原声音的频谱。

<实施方式3>

实施方式3相当于实施方式1的变形例。具体来讲，是进行§2中叙述的基于二阶统计量的信号失真消除处理的方式。因此，对于与实施方式1不同的部分，参照图11和图12附加说明。但是，实施方式3中，设话筒的个数M为2以上。

步骤S100的处理和步骤S101的处理与实施方式1相同。

接着步骤S101的处理，进行步骤S102a的处理。

参照图11说明实施方式3的反向滤波器计算部分(13)的功能构成例。

反向滤波器计算部分(13)的构成为：具有作为对观测信号的预测误差滤波器应用部分功能的第1预测误差滤波器应用部分(1311)；具有作为对观测信号应用更新用反向滤波器所得的信号(更新用反向滤波器应用后信号)的预测误差滤波器应用部分功能的第2预测误差滤波器应用部分(1312)；梯度矢量计算部分(1313)；反向滤波器更新部分(132)；以及更新用反向滤波器应用部分(133)。这里，更新用反向滤波器相当于式(37)的g_m(k)。

第1预测误差滤波器应用部分(1311)将由m号[1≤m≤M]话筒所观测的观测信号x_m(t)进行分帧，对于各帧，将步骤S101的处理所得的i号预测误差滤波器a_i(k)应用于i号帧的信号x_mi(n)而计算预测误差滤波器应用后的信号v_mi(n)(参照式(38))。具体来讲，分帧处理部分(402B)对所输入的观测信号x_m(t)进行分帧处理，输出观测信号x_m(t)的i号帧的信号x_mi(n)。然后，预测误差滤波器应用部分(404i)将信号x_mi(n)作为输入，按照式(38)输出信号v_mi(n)。其中，1≤i≤F。

第2预测误差滤波器应用部分(1312)将更新用反向滤波器应用后的信号y(t)进行分帧，对于各帧，将步骤S101的处理所得的i号预测误差滤波器a_i(k)应用于各帧而计算预测误差滤波器应用后的修正估计值d_i(1)，...，d_i(W)(参照式(30))。再者，更新用反向滤波器应用后的信号y(t)的初始值设为由步骤S100的处理所得的信号即可。具体来讲，除了初始值的情况以外，分帧处理处理部分(402A)对于后述的更新用反向滤波器应用部分(133)所输出的更新用反向滤波器应用后的信号y(t)进行分帧处理，将更新用反向滤波器应用后的信号y(t)进行分帧，输出i号帧的信号y_i(n)。然后，预测误差滤波器应用部分(403i)将信号y_i(n)作为输入，按照式(30)输出修正估计值d_i(1)，...，d_i(W)。其中，1≤i≤F。

梯度矢量计算部分(1313)使用信号v_mi(n)和修正估计值d_i(n)，计算当前的更新用反向滤波器g_m(k)的梯度矢量(参照式(37)的右边第二项)。具体来讲，关于各帧号i(1≤i≤F)，相互相关计算部分(405i)计算信号v_mi(n)和修正估计值d(n)之间的相关<d_i(n)v_mi(n-k)>_n＝1^W。此外，关于各帧号i(1≤i≤F)，方差计算部分(406i)求修正估计值d_i(1)，...，d_i(W)的方差<d_i(n)²>_n＝1^W。关于各帧号i(1≤i≤F)，除法部分(407i)求<d_i(n)v_mi(n-k)>_n＝1^W/<d_i(n)²>_n＝1^W。加法部分(407)求除法部分(4071)～(407F)在全帧中的总和，即式(37)的右边第二项。

反向滤波器更新部分(132)以当前的更新用滤波器g_m(k)、学习率δ、梯度矢量作为输入，按照式(37)求u+1次的更新用反向滤波器g_m(k)′。式(37)意味着将求得的g_m(k)′看作新的g_m(k)来进行更新。

更新用反向滤波器应用部分(133)以通过反向滤波器更新部分(132)所得到的g_m(k)′，即新的g_m(k)和观测信号x(t)作为输入，按照式(42)求更新用反向滤波器应用后的信号y(t)。具体来讲，作为式(42)的g_m(k)，使用由u+1次的更新所得到的g进行计算。该计算所得到的更新用反向滤波器应用后的信号y(t)成为第2预测误差滤波器应用部分(1312)的输入。

接着步骤S102a的处理，进行步骤S103和步骤S104的处理，但由于与实施方式1相同而省略说明。

<实施例2>

发明者们进行了实施方式3的验证实验，因此表示该实验结果。作为实验条件，设M＝4、L＝1000、W＝200、P＝16、R₁＝6、R₂＝50。设定学习率δ，初始值设为0.05，如果∑_i＝1^Flog<d_i(n)²>_n＝1^W的值增加，则通过将学习率δ的值依次减半，使得∑_i＝1^Flog<d_i(n)²>_n＝1^W的值必然减小。反向滤波器的初始估计值设定为g_m(k)＝0、1≤m≤M、1≤k≤L。

以表示声音清晰度的RASTI(参照参考文献5)作为消除混响的指标来评价本发明的实施方式3的效果。从连续发声数据库取出男女各五名的发声，通过叠加在混响时间为0.5秒的混响室测量的脉冲响应而合成观测信号。

(参考文献5)H.kuttruff.Room acoustics.Elsevier Applied Science，thirdedition，P.2371991.

图13是显示了设观测信号长度N为3秒、4秒、5秒、10秒时的RASTI的值的图。如图13所示，可知即使是观测信号为3～5秒的短时间的情况下，也显示着较高的混响消除性能。

图14是消除混响前后的能量衰减曲线的例子。可知在直接声到达50毫秒后的反射声的能量减少了15dB。

工业上的可利用性

本发明是用于提高各种各样的信号处理系统的性能的基础技术，例如可利用于语音识别系统、视频会议系统、助听器、音乐信息处理系统等。