信号歪み除去装置、方法、プログラム及びそのプログラムを記録した記録 媒体
技術分野
[0001] この発明は、信号の歪み(ひずみ: distortion)を除去する技術に関する。
背景技術
[0002] 信号は反射や残響などが存在する環境で観測されると、本来の信号に反射や残響 などが重畳された信号として観測される。以下、本来の信号を「原信号」と云い、観測 された信号を「観測信号」と云うことにする。また、反射や残響などに代表される、原信 号に重畳された歪みを「伝達特性」と云うことにする。このため、観測信号から原信号 固有の特徴を抽出することが困難になる。この不都合を解消すベぐかねてより信号 歪み除去処理技術が各種試みられてきた。信号歪み除去処理は、原信号に重畳し た伝達特性を観測信号力 取り除く処理である。
[0003] 従来の信号歪み除去処理の一例として非特許文献 1に開示されて!、る信号歪み除 去方法を図 15を用いて説明する。予測誤差フィルタ計算部(901)は、観測信号をフ レーム化処理して、各フレームに含まれる観測信号に対して線形予測分析を行 、、 予測誤差フィルタを計算する。この明細書では、フィルタはディジタルフィルタであり、 信号のサンプル値に対して作用するいわゆるフィルタ係数を求める意味で単にフィ ルタを計算するなどということがある。予測誤差フィルタ適用部(902)は、各フレーム ごとに上記計算された予測誤差フィルタを当該フレームの観測信号に適用する。逆 フィルタ計算部(903)は、予測誤差フィルタ適用後の信号に対して逆フィルタを適用 して得られる信号の正規ィ匕尖度が最大となるような逆フィルタを計算する。逆フィルタ 適用部(904)は、上記計算された逆フィルタを観測信号に適用することで信号歪み 除去後の信号 (復元信号)を得る。
特 S午文献 1 : B.W.Giliespie,, H.S.Malvar,,and D.A.F.Florencio,, speech aerever beration via maximum— kurtosis subband adaptive filtering, IEEE International Conf erence on Acoustics, Speech, and Signal Processing, pp.3701— 3704, 2001.
発明の開示
発明が解決しょうとする課題
[0004] 上記の従来的な信号歪み除去方法は、観測信号の各フレーム内ではサンプル間 の相関は原信号固有の特性の寄与が大きぐフレームを跨ぐサンプル間の相関は伝 達特性による寄与が大きいことを仮定している。上記従来方法は、この仮定に基づい て、フレーム化処理された観測信号に予測誤差フィルタを適用して観測信号中の原 信号固有の特性の寄与を低減して 、る。
[0005] しかし、この仮定は粗い近似であるため、逆フィルタの精度は不十分である。つまり 、観測信号力 求まる予測誤差フィルタは伝達特性の影響を受けているので、原信 号固有の特性のみを正しく取り除くことができない。このため、予測誤差フィルタ適用 後の信号力 求める逆フィルタの精度は劣化する。結果として、観測信号に逆フィル タを適用して得る信号は、本来の原信号との比較で精度の良いものではな力つた。 そこで本発明は、精度の良い逆フィルタを求めることで、伝達特性に由来する歪み を観測信号力も除去して精度の良い復元信号を得ることを目的とする。
課題を解決するための手段
[0006] 上記課題を解決するため、本発明の信号歪み除去装置は、所定の繰り返し終了条 件を満たした場合には、観測信号に適用するためのフィルタ(以下、逆フィルタという 。)を、観測信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を 満たさない場合には、観測信号に逆フィルタを適用して、この結果をアドホック信号と して出力する逆フィルタ適用手段と、アドホック信号をフレーム化して、各フレームの アドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力す る予測誤差フィルタ計算手段と、各フレームのアドホック信号に対して当該フレーム に対応する予測誤差フィルタを適用して得る各信号 (以下、イノベーション推定値と いう。)を結合した全フレームでのイノベーション推定値(以下、イノベーション推定値 系列という。)が、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを 出力する逆フィルタ計算手段と、繰り返し終了条件を満たすまで逆フィルタ適用手段 、予測誤差フィルタ計算手段、逆フィルタ計算手段を繰り返し実行させる制御手段と 、を備える。
この発明では、観測信号に伝達特性を除去するための逆フィルタを適用して得るァ ドホック信号に対して、このアドホック信号に基づ 、て求めた予測誤差フィルタを適用 して得る信号 (イノベーション推定値系列)力 その全サンプル間で独立となるような 逆フィルタを求める。そして、所定の繰り返し終了条件を満たしたときの逆フィルタを 観測信号に適用することで復元信号を得る。
[0007] 上記の信号歪み除去装置では、予測誤差フィルタ計算手段は、各イノベーション推 定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、各 イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測 誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない 、各フレームごとの予測誤差フィルタを出力するものであり、逆フィルタ計算手段は、 上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、 各イノベーション推定値の正規ィ匕尖度の全フレームでの総和が最大となるときの逆フ ィルタとして求め、この逆フィルタを出力するものであるとしてもよい。
この構成は、イノベーション系列の独立性の尺度として相互情報量を規定し、これ を最小化する予測誤差フィルタと逆フィルタを交代変数法で求めるものである。この 詳細は後述する。
[0008] あるいは、上記の信号歪み除去装置では、予測誤差フィルタ計算手段は、各イノべ ーシヨン推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、 または、各イノベーション推定値の分散の対数値の全フレームでの総和が最小となる ときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分 析を行ない、各フレームごとの予測誤差フィルタを出力するものであり、逆フィルタ計 算手段は、上記イノベーション推定値系列がその全サンプル間で独立となる上記逆 フィルタを、各イノベーション推定値の分散の全フレームでの総和が最小となるときの 逆フィルタ、または、各イノベーション推定値の分散の対数値の全フレームでの総和 が最小となるときの逆フィルタとして求め、この逆フィルタを出力するものであるとして ちょい。
この構成は、イノベーション系列の独立性の尺度として相互情報量を規定し、これ を最小化する予測誤差フィルタと逆フィルタを交代変数法で求めるものである力 信
号の高次統計量を用いることなく予測誤差フィルタと逆フィルタを交代変数法で求め ることがでさる。
[0009] 上記の信号歪み除去装置では、プリ ·ホワイトユング処理を前置させ、プリ ·ホワイト ニング処理で得られた白色化信号に対して、上記同様の処理を行う装置構成とする ことができる。具体的には、観測信号を線形予測分析して得た白色化フィルタを出力 する白色化フィルタ計算手段と、白色化フィルタを観測信号に適用して白色化信号 を出力する白色化フィルタ適用手段と、所定の繰り返し終了条件を満たした場合に は、白色化信号に適用するためのフィルタ(以下、逆フィルタという。)を、白色化信号 に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない場合 には、白色化信号に逆フィルタを適用して、この結果をアドホック信号として出力する 逆フィルタ適用手段と、アドホック信号をフレーム化して、各フレームのアドホック信号 を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィ ルタ計算手段と、各フレームのアドホック信号に対して当該フレームに対応する予測 誤差フィルタを適用して得る各信号 (以下、イノベーション推定値という。)を結合した 全フレームでのイノベーション推定値(以下、イノベーション推定値系列という。)が、 その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィル タ計算手段と、繰り返し終了条件を満たすまで逆フィルタ適用手段、予測誤差フィル タ計算手段、逆フィルタ計算手段を繰り返し実行させる制御手段と、を備えた信号歪 み除去装置とすることができる。
[0010] 上記課題を解決するため、本発明の方法は、逆フィルタ適用手段が、所定の繰り返 し終了条件を満たした場合には、観測信号に適用するためのフィルタ (以下、逆フィ ルタという。)を、観測信号に適用して、この結果を復元信号として出力し、繰り返し終 了条件を満たさない場合には、観測信号に逆フィルタを適用して、この結果をアドホ ック信号として出力する逆フィルタ適用ステップと、予測誤差フィルタ計算手段が、ァ ドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た 各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、逆フ ィルタ計算手段力 各フレームのアドホック信号に対して当該フレームに対応する予 測誤差フィルタを適用して得る各信号 (以下、イノベーション推定値という。)を結合し
た全フレームでのイノベーション推定値(以下、イノベーション推定値系列という。)が 、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィ ルタ計算ステップと、制御手段が、繰り返し終了条件を満たすまで逆フィルタ適用ス テツプ、予測誤差フィルタ計算ステップ、逆フィルタ計算ステップを繰り返し実行させ る制御ステップと、を有する信号歪み除去方法とする。
[0011] また、上記の信号歪み除去方法では、プリ ·ホワイトユング処理を前置させ、プリ ·ホ ワイトニング処理で得られた白色化信号に対して、上記同様の処理を行う方法とする ことができる。具体的には、白色化フィルタ計算手段が、観測信号を線形予測分析し て得た白色化フィルタを出力する白色化フィルタ計算ステップと、白色化フィルタ適 用手段が、白色化フィルタを観測信号に適用して白色化信号を出力する白色化フィ ルタ適用ステップと、逆フィルタ適用手段が、所定の繰り返し終了条件を満たした場 合には、白色化信号に適用するためのフィルタ(以下、逆フィルタという。)を、白色化 信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない 場合には、白色化信号に逆フィルタを適用して、この結果をアドホック信号として出力 する逆フィルタ適用ステップと、予測誤差フィルタ計算手段が、アドホック信号をフレ ーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの 予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、逆フィルタ計算手段 力 S、各フレームのアドホック信号に対して当該フレームに対応する予測誤差フィルタ を適用して得る各信号 (以下、イノベーション推定値系列という。)を結合した全フレ ームでのイノベーション推定値(以下、イノベーション推定値系列という。)が、その全 サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算 ステップと、制御手段が、繰り返し終了条件を満たすまで逆フィルタ適用ステップ、予 測誤差フィルタ計算ステップ、逆フィルタ計算ステップを繰り返し実行させる制御ステ ップと、を有する信号歪み除去方法とする。
[0012] 本発明の信号歪み除去装置としてコンピュータを機能させる信号歪み除去プロダラ ムによって、コンピュータを信号歪み除去装置として作動処理させることができる。そ して、この信号歪み除去プログラムを記録した、コンピュータに読み取り可能なプログ ラム記録媒体によって、他のコンピュータを信号歪み除去装置として機能させること
や、信号歪み除去プログラムを流通させることなどが可能になる。
発明の効果
[0013] 本発明では、観測信号力 求まる予測誤差フィルタを用いて観測信号中の原信号 固有の特性の寄与を低減するのではなぐ観測信号に (仮の)逆フィルタを適用して 得られるアドホック信号 (仮の復元信号)力 求まる予測誤差フィルタを用いて観測信 号中の原信号固有の特性を低減する。アドホック信号から求まる予測誤差フィルタは 、伝達特性の影響を受けにくいので、原信号固有の特性をより正確に取り除くことが 可能である。このような予測誤差フィルタをアドホック信号に適用して得る信号 (イノべ ーシヨン推定値系列)が全サンプル間で独立となるように求められた逆フィルタは精 度良く伝達特性を除去可能なものであるから、このような逆フィルタを観測信号に適 用することで、伝達特性に由来する歪みが除去された精度の良い復元信号を得るこ とがでさる。
図面の簡単な説明
[0014] [図 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 値の関係を示す図。
1 50
[図 9]Aは残響を含まな 、音声のスペクトログラム、 Bは残響を含む音声のスぺタトログ ラム、 Cは残響除去後の音声のスペクトログラム。
[図 10]Aは残響除去音声の LPCスペクトル歪みの時間変動を説明するためのグラフ
、 Bは対応する区間における原音声信号の抜粋。
[図 11]第 3実施形態に係る信号歪み除去装置(1)の逆フィルタ計算部(13)の機能 構成例を示す機能ブロック図。
[図 12]第 3実施形態における信号歪み除去処理の流れを示す処理フロー図。
[図 13]観測信号長 Nを 3秒、 4秒、 5秒、 10秒としたときの、 RASTIの値を表示した図
[図 14]残響除去前後におけるエネルギー減衰曲線の例を示した図。
[図 15]従来技術を説明するための機能ブロック図。
発明を実施するための最良の形態
[0015] § 1 本発明の理餘
以下、実施形態の説明に先立ち、本発明の理論を説明する。
以下の説明では、特に断りのない限り、信号源は 1つとする。
[0016] 1. 1 信号
本発明の対象となる信号は、人の音声、音楽、生体信号、測定対象物の物理量を センサで観測した電気信号などの信号を広く包含する。より好ましくは、自己回帰 (A utoregressive: AR)過程として表現することができる、ある 、は表現することが好まし い信号であればよい。例えば音声信号は、通常、区分定常な自己回帰過程として表 現される信号、すなわち独立同一分布(U.d.: Independent and Identically Distribut ed)信号に音韻性を表す AR系を作用させた信号として看做される (参考文献 1参照) 以下、信号の代表例として音声信号を挙げて本発明の理論を説明する。 考文献 1) L.R.Rabiner, R.W.bchafer, Digital Processing of Speech Signals , Be 11 Laboratories, Incorporated, 1978.
[0017] 1. 2 音声信号のモデル化
まず、原信号である音声信号 s (t)を、以下の 3つの条件を満足する信号としてモデ ル化する。
[0018] [条件 1]音声信号 s (t)は、区分定常な AR過程で生成される。
この [条件 1]から、 AR過程の次数を P、定常とみなせる区間長を Wサンプルとして
音声信号 s(t)をフレーム化すると、第 iフレームの音声信号 Si (n)は、式(1)のように 表される。式(2)は、第 iフレームの音声信号 s (n)のサンプルと、フレーム化前の音 声信号 s (t)のサンプルとの対応を示している。つまり、第 iフレームの n番目のサンプ ルは、音声信号 s(t)において、(i l)W+n番目のサンプルに相当する。式(1)お よび式(2)において、 b (k)は線形予測係数、 e (n)はイノベーションを表す。但し、 1 ≤n≤W、 l≤t≤N、 Nは全サンプル数である。以下、特に断りの無い限り、パラメ一 タ nは 1フレームのサンプル番号を表し、パラメータ tは全てのサンプル番号を表す。 また、全フレーム数は Fとする。
[数 1]
P
ん =1
[0019] なお、第 iフレームにおける n番目のイノベーション e (n)についても、第 iフレームの n番目のイノベーション e (n)と、フレーム化前の音声信号 s (t)に対するイノべーショ ン e (t)との対応を示すことができる。この場合、第 iフレームの n番目のイノベーション e (n)は、イノベーション e (t)にお!/、て(i— 1) W+n番目のイノベーションに相当し、 e (n)=e( (i 1) W+n)が成り立つ。
[0020] 式(1)を z変換する。左辺の z変換を S (Z)とし、右辺第二項の z変換を E (Z)とし、 B
(z) =∑ Pb (k)z_kとすれば、右辺第一項は、 B (z)S (Z)となる。従って、式(1)の z変換は、(l— B (z))S (Z) =E (Z)である。なお、 z_1は時間領域では 1タップ遅延 素子に相当する。以降、時間領域信号 (タップ重み係数)を小文字で、 z領域信号( 伝達関数)を大文字でそれぞれ表す。 1 Bi(z)は最小位相性を満足しなければなら ず、『1— B (z)は、複素平面上で単位円の内部に全ての零点をもつ』ことが要求され る。
[0021] [条件 2]第 iフレームに属するイノベーション系列 e (1), ···, e (W)は独立且つ同一 分布に属する。イノベーション系列 e (1), ···, e (W)の確率分布の平均及び歪度(3
次キュムラント)は 0、尖度 (4次キュムラント)は正である。さらに、異なるフレーム i、 j〔i ≠j]に属するイノベーション e (n)と e (η' )同士も独立である。ただし、これらは必ず
1 ]
しも同一分布に属するとは限らない。
[条件 3]予測誤差フィルタ 1—Β (ζ)は、相異なるフレーム間で共通する零点をもたな い。
[0022] 式(1)および式(2)から、音声信号 s (t)は、式(3)のように表される。 [·]はガウス記 号である。
[数 2]
P t - l
s(i) =∑bi (k) s(t - k) + e(t) i = + 1 ( 3 )
k=\ W
[0023] このとき、 [条件 2]は、『イノベーション過程 e (t)は時間的に独立な信号である。また 、その統計的性質 (あるいは統計量)はフレーム内では定常である。』と表現できる。 また、 [条件 3]は、『線形予測係数 {b (k) } Pは、時不変な極を持たな!、』と表現でき
i k= l
る。
[0024] 1. 3 観測信号のモデル化
次に、 M個のマイクロホンで音声信号を観測して観測信号を得たときの観測信号を モデル化する。但し Mは、 M≥lの整数である。
m番目(l≤m≤M)のマイクロホンで観測される残響を含む観測信号 X (t)を、音 源力 m番目のマイクロホンに至る経路の伝達関数 H (z)のタップ重み係数 {h (k)
m m
; 0≤ k≤ Κ ; Κはインパルス応答の持続時間とする。 }を用いて式 (4)のようにモデル 化する。ここでは、音声信号の場合の伝達特性の代表例として残響を挙げて、伝達 特性を残響に言い換えて説明する。但し、伝達特性を残響に限定する趣旨ではない
[数 3]
¾ (り (4 )
[0025] M個の観測信号にっ 、てまとめて表現すれば、式(5)のように表すことができる。
但し、式(5)において、 x(t) = [x (k)]
T
である。
画 κ
x(t)= ∑h(k) s(t-k) (5)
k=0
[0026] 1.4 信号歪み除去の原理
信号歪み除去後の復元信号 y(t)は、多チャネル逆フィルタ {G (z); l≤m≤M}の タップ重み係数 {g (k);l≤m≤M, 0≤k≤L;Lは逆フィルタの次数 }を用いて式(6 )により計算される。本発明においては、逆フィルタ係数である g (k)を観測信号 X (t
m l
), ···, X (t)のみ力 推定する。
M
[数 5]
M L
y(t)= ∑ ∑gm(k)xm(t-k) (6)
m=\ k=0
[0027] 1.5 本発明の基本原理
本発明の基本原理は、伝達関数 {H (z);l≤m≤M}の逆フィルタ {G (z);l≤m
m m
≤ M }と AR系のフィルタ {lZ(l B.(z));l≤i≤F}の逆フィルタである予測誤差フ ィルタ {1— A (z) ;l≤i≤F}とを並行して推定することを主な特徴とする。
[0028] この基本原理を説明するため、上述のモデル機構を組み込んだ系全体の構成線 図を図 1に示す。上述のモデルィ匕に拠れば、原信号 s(t)は、フレームごとのイノべ一 シヨン系列 e(l), ···, e (W)に対して AR系のフィルタ 1Z(1— B (z))を適用して得 られる信号 s (n)のフレーム結合と看做すことができ、観測信号 X (t)は、原信号 s (t) に対して伝達関数 H(z)が作用したものと言える。そして、信号歪み除去処理は、観 測信号 X (t)に対して逆フィルタ G (z)を作用させて復元信号 y (t)を得る処理となる。 このとき、信号歪み除去処理で得られた復元信号 y(t)をフレーム分割して、それぞ れに対して、それぞれの信号に基づいて求めた予測誤差フィルタ 1 A (z)を適用し て得られるイノベーション推定値 d(l), ···, d (W)はイノベーション系列 e (l), ···, e (W)に一致することが望ましい。もし、予測誤差フィルタ 1— A (z)の出力信号 d (n)
力 (n) =e (n) [l≤i≤F, l≤n≤W〕を満たすならば、 [条件 3]の条件下で∑ M i i m= 1
H (z) G (z) = lとなることが示せる (数学的証明については、参考文献 Aを参照さ m m
る。
考文献 Aj Takuya Yoshioka, Takafumi Hikicni, Masato iyoshi, Hiroshi G. Ok uno: Robust Decomposition of Inverse Filter of channel and Prediction Error Filter of Speech Signal for D ereverberation , Proceedings of the 14th European Signal Proc essing Conference (EUSIPCO 2006), CD-ROM Proceedings, Florence, 2006.
[0029] しかし、実際には、イノベーション e^n) [l≤i≤F, l≤n≤W〕を信号歪み除去装置 への入力信号として利用できない。図 1に示す系において、各イノベーション系列 ( n)から観測信号 x(t)を得る一連の過程は、モデル過程であって、実際には各イノべ ーシヨン系列 e (n)、フィルタ 1Z (1— B (z) )や伝達関数 H (z)を知ることはできない か知ることが困難であり、利用できる情報は観測信号 x (t)のみである。そこで、上記 [ 条件 2]に基づいて、第 iフレームにおけるイノベーション推定値 d (1) , · ··, d (W)を 結合して得る全フレームでのイノベーション推定値系列力 その全サンプル間で独立 になるように、つまりイノベーション推定値系列 d (1) , · ··, d (W) , · ··, d (l) , · ··, d (W) , · ··, d (1) , · ··, d (W)が独立となるように逆フィルタ G (z)と予測誤差フィルタ
F F m
1 A (z)を推定する。
[0030] このことを、従来手法との対比で述べる。従来手法は、逆フィルタを、「観測信号に 基づいて求めた予測誤差フィルタを観測信号に適用し、予測誤差フィルタ適用後の 信号に対して逆フィルタを適用して得られる信号の正規化尖度が最大となる逆フィル タを求めよ」という問題の解として得ていた。これに対して、本発明は、逆フィルタを、「 観測信号に逆フィルタを適用して得る信号に対して、当該信号に基づいて求めた予 測誤差フィルタを適用して得る信号力 全サンプル間で独立となる逆フィルタを求め よ」という問題の解として得る。この問題で留意しなければならないことは、予測誤差 フィルタが、観測信号に逆フィルタを適用して得る信号に基づ 、て求められるため、 逆フィルタだけでなく予測誤差フィルタも一緒に求めることになるということである。 この問題は、 ICA (Independent Component Analysis)と同様の考え方によって定式
化することができる。ここでは相互情報量を最小化する観点から説明を行うが、例え ば対数尤度を最大化する方針で定式ィ匕することも可能である。 V、ずれにしても問題 に対するアプローチ方法の違 ヽに過ぎな、、。
[0031] 独立性の尺度として相互情報量 (Kullback-Leibler情報量)を用いると、解くべき問 題は式 (7)のように定式化される。ただし、 g=[gT, ···, g T]T, g =[g (0), ···, g
1 M m m m
(L)]T, &=[&ι τ, ···, ap T]T, a = [a(l), ···, ^(Ρ)]τとし、 ^(k)は予測誤差フィルタ 係数を表す。 KU^ ···, Un)は確率変数 間の相互情報量を表す。また gおよび aに 記号'を付したものは、得るべき最適解である。 Tは転置を表す。
[数 6]
[g,a] = arg min I\dx (1),···,^ (W),---,dF (1),■•■,dF{W)) (7) 拘束条件
[l] II g II =i (但し II · IIはノルムを表す。)
[2] 1— (2)は、複素平面上で単位円の内部に全ての零点をもつ〔l≤i≤F〕。
[0032] 相互情報量 Iは、イノベーション推定値系列 d (1), ···, d (W), ···, d(l), ···, d( W), ···, d (1), ···, d (W)の振幅が定数倍されても変化しない。式 (7)の拘束条件
F F
[1]は、この振幅の不定性を排除するための条件である。式 (7)の拘束条件 [2]は、上 記 [条件 1]に対応して、予測誤差フィルタを最小位相系に制限するための条件である 。以下、 Iを、イノベーション推定値系列を入力としそれらの間の相互情報量を出力す る関数と看做して、損失関数と呼称することにする。
[0033] 1.6 損失関数の導出
式 (7)の最適化を実行するためには、損失関数 I(d (1), ···, d (W))を有限長の
1 F
信号系列 {d(n);l≤i≤F, l≤n≤W}から推定しなければならない。(多変量)確率 変数 Uの微分エントロピーを D(U)と表記すると、 I(d (1), ···, d (W))は式 (8)で定
1 F
義される。ただし、 d=[d τ, ···, dTl\d = [d (W), ···, d(l)]Tである。
[数 7]
F W
I{dx (l),-,dFm) =∑∑。( (n))-D{d) (8)
i—\ n=\
[0034] y=[yF,…, γΎ, y = [y.(W), ···, y^l)]1とおくと、 dは yを用いて、 d=Ayと表 される。ただし、行列 Aは、式(9)および式(10)で表される。
[数 8]
[0035] よって、 D(d)は式(11)のように表される。
[数 9]
D(d) = D(y) + \ogdetA (11)
[0036] 多変量確率変数 Uの共分散行列を∑ (U)と表記すると、式(11)右辺第二項につ いて、∑ (d)=E{ddT}=AE{yyT}AT=A∑ (y) Ατが成立するから、式(12)が成り 立つ。
[数 10] log detA = - (log det∑ (め― log det∑ ( ) (12)
ム
[0037] 式(11)、式(12)を式 (8)に代入すると、式(13)を得る。ただし、 σ (U)2は確率変
数 uの分散を表す。
[数 11] - β (
F W
ニー∑∑ Adi (")) + C{d, (l),-,dF(W)) + J(y) (1 3)
i=l n=\ 式( 13)で J (U)は(多変量)確率変数 Uのネゲントロピー(negentropy)である。ネゲ ントロピーは Uの非ガウス性の度合いを表す非負の値をとり、 Uがガウス分布に従う場 合に限り 0をとる。 C(U , ···, U )は式(14)で定義される。 C(U , ···, U )は確率変 数 U間の相関の度合いを表す非負の値をとり、これらが無相関の場合に限り 0をとる
[数 12] _logdet∑ ([ひい…,ひ„] ) (14)
ところで、 s=[s T, ···, sつ"1, s = [s (W) , ···, s (l)]1とおく ^J(y)=J(s)= const antとなるため(証明略)、式( 13)は更に式( 15)のように簡単化できる。
[数 13] l(dx{\),---,dF{W))
F W (1 5)
= - ∑ゾ ( (")) +C(dx{\),---,dF(W))+ const
i-\n=\
[0040] 以上から、式(16)の最適化問題を解けばよいことになる。
[数 14] g,a) - ゾ ( ( ) +C ( (1),…,^( ) (l 6)
[1] II g II =1 (但し II · IIはノルムを表す。)
[2] 1 A (z)は、複素平面上で単位円の内部に全ての零点をもつ〔l≤i≤F〕。
[0041] 1.7 交代変数法による最適化
式(16)について、交代変数の方法により、 gと aを最適化する。すなわち、 r回目の 繰り返しにおける g及び aの推定値をそれぞれ g"W、 ωと表せば、式(17)および式 (18)の交互の最適化により更新された推定値 g"(I+1), a"(I+1)を得る。なお、 g"およ び a "は、記号"が g、 aのそれぞれの上に付されたものを表す。例えば繰り返し回数の 上限を Rとすれば、 R回目で得られる g'(R1+1)、a'(R1+1)が式(16)の最適解である。 上付き文字の R1は、 Rである。
[1] g=g"(r)
[2] 1-A (z)は、複素平面上で単位円の内部に全ての零点をもつ〔l≤i≤F〕< [数 16] g d F, (1 8)
拘束条件
[1] a = a"(r+1)
[2] II g || =i
[0042] 式(17)の意図するところは、伝達特性を打ち消すための逆フィルタの現在の推定
値に基づいて原信号に固有の特性を打ち消すための予測誤差フィルタを推定するこ とである。同様に、式(18)の意図するところは、予測誤差フィルタの現在の推定値に 基づいて逆フィルタを推定することである。イノベーション推定値系列 d (1), ···, d ( W), ···, d (1), ···, d (W), ···, d (1), ···, d (W)が互いにより独立になるようにこ
i i F F
れら 2種類の最適化を繰り返すことで、逆フィルタと予測誤差フィルタを並行して推定 することが可能になっている。したがって、ここでの繰り返しは逆フィルタの高精度な 推定のために重要である。但し、図 8から明らかなように処理する観測信号長が長く なる程、繰り返し回数は 1回でも或る程度の効果が得られることが見て取れる。従って 、この発明では、繰り返し回数は 1回でもよい。
[0043] 1.8 aの最適化
本発明では、式(17)の最適化を以下のように行う。
まず注意すべきことは、 C(d (1), ···, d (W))は d(n)の二次の統計量に関連する
1 F i
のに対して、 J (d (n) )は d (n)の高次の統計量に関連する値である。二次の統計量 は信号の振幅情報のみ提供するが、高次の統計量は位相情報も提供する。したがつ て、一般に、高次統計量を含む最適化は、非最小位相系を導く可能性がある。そこ で、 1— A (z)が最小位相系であるという拘束条件から、 aの最適化においては式(19 )の最適化問題を解く。
[数 17] a (r+1)= argminC ( (1),···,ί/ ( )) (1 9)
a 拘束条件
[l] g=g"(r)
[2] 1 A (z)は、複素平面上で単位円の内部に全ての零点をもつ〔l≤i≤F〕。
[0044] C(d (1), ···, d (W))は式(20)で与えられる。
1 F
[数 18]
[0045] ここで、行列 Aは式(9)および式(10)に示すように上三角行列でその対角成分が すべて 1であるから、 log det A=0である。これを式(12)に代入することで式(21)の 関係を得る。
[数 19] logdet∑ (め = logdet∑0) = constant ( 2 1 )
[0046] よって、式(19)は、式(22)の最適化問題と等価である。なお、式(22)は、上記 [条 件 2]を反映した表現になっていることに留意しなければならない。式 (22)を説明す れば、式(22)は、第 iフレームにおけるイノベーション推定値 d (1) , · · ·, d (W)の分 散の対数値を全フレームで加算した値が最小となる aを求めよ、と云っている。
[数 20]
F W
ά argmin∑∑logcr ( ("))2 ( 2 2 )
a i=\ η-ϊ 拘束条件
[1] g=g" (r)
[2] 1 A (z)は、複素平面上で単位円の内部に全ての零点をもつ〔l≤i≤F〕。
[0047] 式 (22)で表される最適化問題を解くことは、観測信号に g" Wで与えられる逆フィル タを適用して得られるアドホック信号に対して、各フレームにおいて線形予測分析を 行うことと等価であり、必ず最小位相予測誤差フィルタを得ることができる。線形予測 分析に関しては、上記参考文献 1を参照されたい。
[0048] なお、式(22)では、第 iフレームにおけるイノベーション推定値 d (l) , · · ·, d. (W)の 分散の対数値の全フレームでの総和が最小となるときの aを a" として求める力 こ れに限定する趣旨ではな 、。上記各式では対数関数の底 (base)を明記して!/、な!/、 力 一般的には底を 10ないしネィピア数とするのが慣例であり、いずれにしても底は 1よりも大きい。この場合、対数関数は単調増加関数であるから、第 iフレームにおけ るイノベーション推定値 d (l) , · · ·, d (W)の分散の全フレームでの総和が最小となる ときの aを として求めることができる。
[0049] 1. 9 gの最適化
本発明では、式(18)の最適化を以下のように行う。
前述したとおり、 C(d (1), ···, d (W))は {d (n) ;l≤i≤F, l≤n≤W}の相関の
1 F i
度合いに関わる指標であるが、(r + 1)回目の aの最適化にぉ 、て最小化されて 、る ため、∑ F∑ wJ(d (n))に比べて無視できる。そこで gの最適化においては、式 ( 23)の最適化問題を解く。
[数 21] g (23)
[1] a = a"(r+1)
[2] II g || =ι
[0050] J (d (n) )は、 [条件 2]に基づ 、て、式(24)によって近似できる。この詳細は参考文 献 2を参照されたい。ただし、確率変数 Uについて、 κ (U)は Uの尖度 (4次キュムラ
4
ント)を表す。式(24)の右辺を第 iフレームにおける正規ィ匕尖度と!/、う。
(参考文献 2) A.Hyvarinen, J.Karhunen, E.Oja, "INDEPENDENT COMPONENT A
NALYSIS", John Wiley & Sons, Inc.2001.
[数 22] ) ¾^ 4)
σ ( ( )
[0051] [条件 2]力 音声信号のイノベーションの尖度は正であるため、 κ (d (η))/σ (d
4 i i
(n)) 4は正である。従って、式(23)の最適化問題は、式(25)の最適化問題に帰着 する。 σ (d (η)), κ (d(n))は、 [条件 1]で述べた音声信号の局所的な定常性に
i 4 i
基づいて、各フレーム内のサンプル力も計算される。式(26)では、 1ZWを付カ卩して いるが、これは後の計算の便宜に過ぎず、式(25)で gの最適解を求めるにあたり影 響を及ぼすものではない。式(25)および式(26)から、正規化尖度の全フレームで の総和が最大となるときの gが、 g'(I+1)となる。なお、式(25)および式(26)は、上記 [ 条件 2]を反映した表現になっていることに留意しなければならない。式(25)および
式(26)を説明すれば、これらは、第 iフレームにおける正規ィ匕尖度を全フレームで加 算した値が最大となる gを求めよ、と云っている。
[数 23] g (r+1)=argmaxg (25)
[1] a = a"(r+1)
[2] II g || =1
[0052] 式(25)に従って gの最適解を求めるには、 Qを gで微分してこれをゼロとしたときの 解を求めればよい。この解は、一般的には、式(27)で表される更新則に従って求め られる。 g' を g' のノルムで除しているのは上記拘束条件 [2]を課すためである。 r? ( u)は学習率を表す。 uは、 gの最適化における更新回数を表す。
[数 24]
[0053] 式(27)〖こおいて、 VQは式(28)および式(29)で与えられる
g
[数 25] dQ dQ dQ dQ
(28)
¾i(0) dgl(L) dgM(0) dgM(L)j
— (")4 (")2 } ("一 ))
[0054] 式(29)において、 d (n)は式(30)で、 v (n)は式(31)および式(32)で与えられ
i mi
る。 X (n)は、 m番目のマイクロホンで観測された観測信号の i番目のフレームの信 号である。
[数 26]
P
di (n) = yt (n) -∑ at (k) yt (n - k) ( 3 0 )
k=\
P
v - (") = xmi {n)― X α,· (k) xmi (n - k) ( 3 1 )
[0055] § 2 二次統計量に基づく信号歪み除去
上述の従来的手法の信号歪み除去方法は、比較的長時間の観測信号 (例えば 20 秒程度である。)を要する。これは、一般に、正規化尖度のような高次統計量を計算 するためには大量の観測信号のサンプルが必要となるからである。しかし、実際には そうした長時間の観測信号を利用できない場合が多い。このため、従来的手法の信 号歪み除去方法の適用分野は極めて限られて!/ヽた。
また高次統計量の計算は比較的複雑であるため、従来的手法の信号歪み除去方 法では装置の構成が複雑になりやす 、。
そこで、観測信号がより短時間(例えば 3秒から 5秒程度である。)の場合にも有効 であり、かつ計算が従来に比して容易な信号歪み除去の原理を説明する。この原理 は、信号の二次統計量のみを用いるものであり、 § 1で説明した本発明の基本原理よ り派生する。
[0056] 2. 1 二次統計量に基づく信号歪み除去の原理
二次統計量に基づく信号歪み除去では、上述の 3つの条件に、次の 2つの条件を 設定する。
[条件 4] M≥2である。すなわち、複数本のマイクロホンを用いる。
[条件 5] H = {h (k) } κは相異なるマイクロホン mの間で共通の零点を持たない
m m k=0
[0057] 上記の式(16)の最適化問題では、高次の統計量に関する値であるネゲントロピー Jおよび確率変数間の相関の度合いを示す指標 Cを含む値を最小化する gおよび aを 求めた。
確率変数間の相関の度合いを示す指標 Cは、二次の統計量で規定される。そこで 、解くべき最適化問題を式 (33)で定式ィ匕する。
[数 27] g, a} = argmin ^(Ι),· ··,^^))
- ogCT ( ("))2 - logdet (め ( 3 3 )
[0058] 式(21)を参酌すれば、式(33)の最適化問題は、式(34)の最適化問題に転ィ匕さ れる。なお、式(34)は、上記 [条件 2]を反映した表現になっていることに留意しなけ ればならない。式(34)を説明すれば、式(34)は、第 iフレームにおけるイノべーショ ン推定値 d (l) , · · ·, d (W)の分散の対数値を全フレームで加算した値が最小となる gおよび aを求めよ、と云っている。
[数 28]
F W
g,a) = argmin ∑∑logび ( ( 3 4 )
g," =l "=l
[0059] ところで、上記の [条件 4]および [条件 5]が成立する場合、多チャンネルの観測信号 は、音源力 の原信号によって駆動される AR系として捉えることができる (参考文献 3参照)。このことは、逆フィルタ Gの先頭タップを式(35)のように固定できることを意 味する。但し、 m= lに相当するマイクロホンは、最も音源に近いマイクロホンである。
考文献 3) K. Aded-Meraim, E. Moulines, and P. Loubaton. Prediction error me thod for second-order blind identification. IEEE Trans. Signal Processing, Vol. 45, No.3, pp. 694-705, 1997.
[0060] 式(34)および式(35)で規定される gを係数とする逆フィルタ Gを、式(6)に従って 観測信号 x(t)に適用することで伝達特性が除去された復元信号 y(t)を得る。
[0061] 2.2 aの最適化
式(34)について、交代変数の方法により、 gと aを最適化する。
逆フィルタの係数 g (k)を固定した状態で予測誤差フィルタの係数 a (k)に関して 式 (34)の損失関数を最小化する。
このとき、次の 2点に注意する。 1点目は、 g=[gT, ···, g τ]τは固定されているの で、逆フィルタ Gの出力である復元信号 y(t)は予測誤差フィルタの最適化において 不変であるということである。 2点目は、 i番目のフレームにおける予測誤差フィルタの 係数 a(l), ···, a(P)は、 d(l), ···, d(W)にのみ寄与するということである。
[0062] 従って、各々のフレームにお 、て、∑ wlog σ (d (n) ) 2を最小化するように、予測 誤差フィルタの係数 a (1), ···, a (P)を推定すればよい。 [条件 2]からイノベーション 推定値 d (1), ···, d (W)の分散はフレーム内で定常であるから、∑ wloga (d (n) )2の最小化は、 W* σ (d (n)) 2の最小化と等価である。記号 *は乗算を表す。分散 σ (d (n))2は、く d (n)2> wとして計算できる。但し、 <d (n)2> wは、 1フレー ム分のイノベーション推定値 d (1), ···, d (W)を使って計算した d (n)の 2乗平均を 表す。結局、係数 a (k)は、 W* <d (n)2> w、つまり d (n)の 2乗総和が最小にな るときの &i(k)として推定される。このような係数 &i(k)は、線形予測分析の手法によつ て計算される。
[0063] なお、ここでは、第 iフレームにおけるイノベーション推定値 d (1), ···, d (W)の分 散の対数値の全フレームでの総和が最小となるときの aを a" として求めるとして説 明したが、これに限定する趣旨ではない。既述のとおり、第 iフレームにおけるイノべ ーシヨン推定値 d(l), ···, d (W)の分散の全フレームでの総和が最小となるときの a を a' として求めることができる。
[0064] 2.3 gの最適化
予測誤差フィルタの係数 a (k)を固定した状態で逆フィルタの係数 g (k)に関して
式 (34)の損失関数を最小化する。
逆フィルタの係数 g (k)に関する損失関数の最小化には、勾配法を用いる。 [条件 2]を用いると、式(34)の最適化問題は、式(36)の最適化問題に転化される。
[数 30]
( F W ヽ
g = ar min ∑log ( (" ) (36)
g V/=i ヽ ln=l J
[0065] 式(36)に従って gの最適解を求めるには、∑ 1'log<d(n)"> wを gで微分して これをゼロとしたときの解を求めればよい。この解は、一般的には、式(37)で表され る更新則に従って求められる。 δは学習率を表す。 l≤m≤M、 l≤k≤Lである。な お、式 (37)では、式 (35)の条件から II g II =1なる拘束条件が課されないことに留 意しなければならない。同様に式(35)の条件力 kが取る値の範囲は l≤k≤Lであ る。
[数 31]
Sm ( ) = Sm (ん) + δ (37)
P
vmi (") = xmi (n) -∑ ( ) Xmi ("-ん) (38)
k=\
[0066] 上記式(29)ある 、は上記非特許文献 1に記載された式(3)と比べて明らかなように 、式(37)の右辺第二項は、二次の統計量で表されており、この計算に高次の統計量 を必要としない。このため、高次統計量を計算するには短い時間の観測信号の場合 にも有効であり、計算自体も容易である。
[0067] なお、式(36)では、第 iフレームにおけるイノベーション推定値 d(l), ···, d(W)の 分散の対数値の全フレームでの総和が最小となるときの gを g"として求める力 これ に限定する趣旨ではな 、。上記各式では対数関数の底 (base)を明記して!/、な 、が、 一般的には底を 10ないしネィピア数とするのが慣例であり、いずれにしても底は 1より
も大である。この場合、対数関数は単調増加関数であるから、第 iフレームにおけるィ ノベーシヨン推定値 d(l), ···, d(W)の分散の全フレームでの総和が最小となるとき の gを g"として求めることができる。なお、この場合には式(37)で示した更新則は適 用できなくなり、改めて∑ F<d(n)2> wを gで微分してこれをゼロとしたときの解 を求めればよい。この結果得られる更新則は ICAと同様の考え方で定式ィ匕できるか ら略 る ο
[0068] §3 プリ'ホワイトユング
本発明の信号歪み除去には、プリ'ホワイトユングを適用することができる。観測信 号をプリ'ホワイトユングすることで、最適化計算の安定化、とくに更新則の高速な収 束が可能となる。
各マイクロホンで得られた観測信号系列全体 {x (t) ;l≤t≤N}を白色化するフィ ルタ(白色化フィルタ)の係数 {f (k) ;0≤k≤X}を X次の線形予測分析によって計算 する。
式(39)に従って、上記白色化フィルタを各マイクロホンで得られた観測信号 X (t) に適用する。 w (t)は、 m番目のマイクロホンで得た観測信号 X (t)を白色化した信
m m
号を表す。
[数 32]
X
^( =∑fm(k)xm(t-k) (39)
k=0
[0069] このとき、式(31)および式(38)は式(40)に、式(32)は式(41)に、変更すればよ い。
[数 33]
P
vmi(n) = wmi(n)― ^ ( ) wmi^n-k) (40)
k=l
wmi(n) = wm((i-l)W + n) (41)
[0070] §4 実施形態
以下、本発明の実施形態を図面を参照して説明する。本発明の実施形態として、
後述の各実施形態に限定するものではなぐ各セクションで説明した原理を実現する 実施形態であればよい。
[0071] <第 1実施形態 >
本発明の第 1実施形態を実施する場合、以下の手順に従ってセンサで得た観測信 号を処理する。ここでは、実施形態を具体的に説明する観点力 信号として音声信 号を例に挙げて説明する。
なお、第 1実施形態の説明に先立ち、観測信号およびフレーム化処理について概 説する。
[0072] ((観測信号))
図示しないセンサ(例えばマイクロホン)によって得られたアナログ信号 (このアナ口 グ信号には伝達特性に由来する歪みが重畳されている。)は、例えば 8,000Hzのサ ンプリングレートでサンプリングされ、適宜量子化された離散信号に変換される。以下 、この離散信号を観測信号ということにする。アナログ信号から観測信号への AZD 変換などを実行するために必要となる構成要素(手段)は、 V、ずれも公知技術の常套 手段によって達成されるから、説明および図示を略する。
[0073] ((フレーム化処理))
図示しない信号切出手段が、離散信号から、時間軸方向に一定時間幅で始点を 移動させながら、所定時間長の離散信号をフレームとして切り出す。例えば 200サン プル点(8,000Hz X 25ms)長の離散信号を、 80サンプル点(8,000Hz X 10ms)ず つ始点を移動させながら切り出す。各フレームの切り出しは、離散信号に公知の窓 関数 (例えば、ハミング窓、ガウス窓、方形窓など)を適用すればよい。窓関数の適用 によるフレームの切り出しは公知の常套手段によって達成される。
[0074] 本発明の第 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 (1 4)、 RAM ( 15)、 ROM ( 16)、外部記憶装置( 17)間のデータのやり取りが可能なよ うに接続するバス(18)を有している。また必要に応じて、信号歪み除去装置(1)に、 CD-ROM (Compact Disc Read Only Memory)、 DVD (Digital Versatile Disc)など の記憶媒体を読み書きできる装置 (ドライブ)などを設けるとしてもよ!/、。
[0075] 信号歪み除去装置(1)の外部記憶装置(17)には、信号歪み除去のためのプログ ラムおよびこのプログラムの処理にぉ 、て必要となるデータ (観測信号)などが記憶さ れている〔外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置である ROMに記憶させておくなどでもよい。〕。また、これらのプログラムの処理によって得 られるデータなどは、 RAMや外部記憶装置などに適宜に記憶され、他のプログラム の処理に供されるときに、 RAMや外部記憶装置など力 読み込まれる。
[0076] より具体的には、信号歪み除去装置(1)の外部記憶装置(17)〔あるいは ROMなど 〕には、観測信号に逆フィルタを適用する処理のためのプログラム、観測信号に逆フ ィルタを適用して得られた信号力 予測誤差フィルタを求める処理のためのプロダラ ム、予測誤差フィルタから逆フィルタを求める処理のためのプログラム、およびこれら のプログラムの処理にぉ 、て必要となるデータ(フレーム単位の観測信号など)が記 憶されている。その他、これらのプログラムに基づく処理を制御するための制御プログ ラムも適宜に保存しておく。
[0077] 第 1実施形態に係る信号歪み除去装置(1)では、外部記憶装置(17)〔あるいは R OMなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必 要に応じて RAM (15)に読み込まれて、 DSP ( 14)で解釈実行 '処理される。その結 果、 DSP (14)が所定の機能 (逆フィルタ適用部、予測誤差フィルタ計算部、逆フィル タ計算部、制御部)を実現することで、信号歪み除去が実現される。
[0078] そこで次に、図 3〜図 5を参照して、信号歪み除去装置(1)における信号歪み除去 処理の流れを順次説明する。
大ま力な処理の手順は、(a)観測信号 x (t)に対して逆フィルタを適用した信号 (以 下、アドホック信号という。)を求め、(b)アドホック信号から予測誤差フィルタを求め、 (c)この予測誤差フィルタ力 逆フィルタを求め、(d)前記(a)、(b)、(c)の処理を繰り 返して最適な逆フィルタを求め、(e)最適化された逆フィルタを観測信号に対して適 用した信号を復元信号 y (t)として得る。
(b)は上述の aの最適化に相当し、(c)は上述の gの最適化に相当し、(d)は、式(1 7)および式(18)に相当する。(d)の処理の繰り返し回数は予め定めた回数 Rとする 。つまり、 l≤r≤Rとする。また、(c)の処理で gを最適化する更新則の更新回数は 予め定めた回数 Rとする。つまり、 l≤u≤Rとする。(d)の処理、つまり(a)、(b)、 (c
2 2
)の一連の処理を 1回行うたびに、更新則による R回の更新が行なわれる。実施形態
2
では、回数 Rは、予め定めた回数とするが、これに限定されず、例えば、 r回目の gを 算出したときの式 (26)の Qの値と r+ 1回目の gを算出したときの式(26)の Qの値との 差の絶対値が所定の正の微小値 ε以下 (あるいは ε未満)になったときに繰り返しを 中止するようにしてもよい。同様に、回数 Rは、予め定めた回数とする力 これに限定
2
されず、例えば、 u回目の gを算出したときの式(26)の Qの値と u+ 1回目の gを算出 したときの式(26)の Qの値との差の絶対値が所定の正の微小値 ε以下 (あるいは ε 未満)になったときに繰り返しを中止するようにしてもょ 、。
[0079] (ステップ S 100)
逆フィルタ適用部(14)は、式 (42)に従って逆フィルタを観測信号 X (t) = [X (t) , · · ·, χ (t) , · · ·, x (t) ]Tに適用することで、アドホック信号 y (t)を求める。アドホック m M
信号 y (t)は、計算上は復元信号と全く同じであるが、ここでは後述の R回の処理を 経て求められた復元信号ではな 、ことを明示するためアドホック信号と呼称する。ここ で tは、全てのサンプル番号を示し、 l≤t≤Nとする。 Nは全サンプル数である。第 1 実施形態では、マイクロホンの個数 Mは 1以上であればよい。
[数 34]
M L
y(t) = ∑ ∑gm (k) xm (t - k) ( 4 2 )
m=l k=0
[0080] 逆フィルタの係数列 {g (k) ; 0≤k≤L}として、繰り返し回数 Rの初回には予め定
められた初期値を、 2回目以降には後述する逆フィルタ計算部(13)によって求めら れた逆フィルタ g' を用いる。
[0081] (ステップ S 101)
予測誤差フィルタ計算部(15)は、フレーム化処理を行うフレーム化処理部(151)と フレーム予測誤差フィルタ計算部(152)によって構成される。そして、フレーム予測 誤差フィルタ計算部(152)は、第 iフレームのアドホック信号力も予測誤差フィルタを 求める第 iフレーム用予測誤差フィルタ計算部(152i)力もなる。ただし、 iは、 l≤i≤F を満たす整数である。
[0082] フレーム化処理部(151)は、逆フィルタ適用部(14)で求められたアドホック信号 {y
(t) ; l≤t≤N}をフレーム化処理する。フレーム化処理は、例えば式 (43)のように、 W点分を切り出す窓関数を W点ずつシフトさせて適用することにより行う。 {y (n) ; 1 ≤n≤W}は i番目のフレームに含まれるアドホック信号列を表す。
[数 35] yi(n) = y((i ~ W + n) ( 4 3 )
[0083] そして、第 iフレーム用予測誤差フィルタ計算部(152i)は、式 (22)に従って、第 iフ レームのアドホック信号列 {y (n) ; l≤n≤W}に対して P次の線形予測分析を行い、 予測誤差フィルタの係数列 {a (k) ; l≤k≤P}を計算する。この算出方法は、上記参 考文献 1を参照されたい。ここで得られた a (1) , · ··, a (P) , · ··, a (1) , · ··, a (P) , · ··, a (1) , · ··, a (P)は、式(22)の a' 1)を与える。
F F
[0084] (ステップ S 102)
逆フィルタ計算部(13)の機能構成例を、図 4を参照して説明する。逆フィルタ計算 部(13)は勾配計算部(131)、逆フィルタ更新部(132)および更新用逆フィルタ適用 部(133)によって構成される。更に、勾配計算部(131)は、観測信号への予測誤差 フィルタ適用部として機能する第 1の予測誤差フィルタ適用部(1311)と、観測信号 に更新用逆フィルタを適用して得られる信号 (更新用逆フィルタ適用後信号)への予 測誤差フィルタ適用部として機能する第 2の予測誤差フィルタ適用部(1312)と、勾 配ベクトル計算部(1313)とを備えて構成される。ここで更新用逆フィルタは、式(27)
の g<U〉に相当する。
[0085] 第 1の予測誤差フィルタ適用部(1311)は、 m番目〔l≤m≤M〕のマイクロホンで観 測された観測信号 X (t)をフレーム化して、各フレームにっき、 i番目のフレームの信 号 X (n)に対してステップ S 101の処理で得られた i番目の予測誤差フィルタ a (k)を 適用して予測誤差フィルタ適用後の信号 v (n)を計算する (式 (31)を参照)。ここで 述べた処理の詳細の一例は、後述の第 3実施形態の説明に譲る。
[0086] 第 2の予測誤差フィルタ適用部(1312)は、更新用逆フィルタ適用後信号 y(t)をフ レーム化して、各フレームにっき、 i番目のフレームの信号 y (n)に対してステップ SI 01の処理で得られた i番目の予測誤差フィルタ a (k)を適用して予測誤差フィルタ適 用後のイノベーション推定値 d (1) , · ··, d (W)を計算する (式 (30)を参照)。なお、 更新用逆フィルタ適用後信号 y(t)の初期値は、ステップ S 100の処理で得られた信 号とすればよい。爾後、第 2の予測誤差フィルタ適用部(1312)は、後述する更新用 逆フィルタ適用部(133)が出力した更新用逆フィルタ適用後信号 y(t)を入力とする 。ここで述べた処理の詳細の一例は、後述の第 3実施形態の説明に譲る。
[0087] 勾配ベクトル計算部(1313)は、信号 V (n)とイノベーション推定値 d (n)とを用い て現在の更新用逆フィルタ g<u〉の勾配ベクトル VQを計算する(式(28)および式(2
9)を参照)。有限個のサンプル V (n)および d (n)を用いて式 (29)を演算するときは
、期待値 Eをサンプルから求めればよい。ここで述べた処理の詳細の一例は、後述の 第 3実施形態の説明に譲る。
[0088] 逆フィルタ更新部(132)は、現在の更新用逆フィルタ g<u〉、学習率 7? (u)、勾配べク トル VQを入力として、式(27)に従って、 u+ 1回目の更新用逆フィルタ 8 <11+1〉を求め る。式 (27)は、求められた g<u+1〉を新たな g<u〉と見立てて更新を行なうことを意味する
[0089] 更新用逆フィルタ適用部(133)は、逆フィルタ更新部(132)によって得られた g<u+1 、、つまり新たな g<u〉および観測信号 x (t)を入力として、式 (42)に従って、更新用逆フ ィルタ適用後信号 y(t)を求める。具体的には、式 (42)の g (k)として u+ 1回目の更 新で得られた gを用いて計算する。この計算で得られた更新用逆フィルタ適用後信号 y (t)は、第 2の予測誤差フィルタ適用部(1312)の入力となる。なお、更新用逆フィ
ルタ適用後信号 y(t)は、計算上は復元信号と全く同じであるが、ここでは後述の 回の処理を経て求められた復元信号ではなぐ更新則を行なうために算出される信 号であることを明示するため更新用逆フィルタ適用後信号と呼称する。
[0090] 制御部(600)の制御によって R回の更新が行なわれた結果として得られた g<R2+1〉
2
は、式(25)の g' (I+1)に相当する。上付き文字の R2は、 Rである。逆フィルタ計算部(
2
13)は、 g' (I+1)を出力する。
[0091] 制御部(500)の制御によって、上述の一連の処理を 1回行うごとに rに 1を加算して rが Rに等しくなるまで、つまり上述の一連の処理を R回繰り返すことで (ステップ S1 〇3)、 g' (R1 +1)を得る。上付き文字の R1は、 Rである。この g' (R1+1)が、式 6)の最 適解とされる。そこで、 g" (R1+1)を得た段階で、逆フィルタ適用部(14)は、式 (42)に 従って逆フィルタ g" (R1+1)を観測信号 x (t) = [x (t) , · ··, x (t) ]Tに適用することで、
1 M
復元信号 y (t)を得ることができる (ステップ S 104)。
[0092] <第 2実施形態 >
第 2実施形態は、第 1実施形態の変形例に相当する。具体的には、 § 3で述べたプ リ 'ホワイトユングを行なう形態である。そこで、第 1実施形態と異なる部分について図 6および図 7を参照して説明を加える。なお、プリ'ホワイトニングは観測信号に対して 行なうプリ ·プロセスであるから、ここで説明するプリ ·ホワイトユングを行なう形態は、 後述の第 3実施形態にも適用可能である。
第 2実施形態では、信号歪み除去装置(1)の外部記憶装置(17)〔あるいは ROM など〕に、白色化フィルタを求める処理のためのプログラム、白色化フィルタを観測信 号に適用する処理のためのプログラムも記憶されている。
第 2実施形態に係る信号歪み除去装置(1)では、外部記憶装置(17)〔あるいは R OMなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必 要に応じて RAM (15)に読み込まれて、 DSP ( 14)で解釈実行 '処理される。その結 果、 DSP (14)が所定の機能 (逆フィルタ適用部、予測誤差フィルタ計算部、逆フィル タ計算部、白色化フィルタ計算部、白色化フィルタ適用部)を実現することで、信号歪 み除去が実現される。
[0093] (ステップ S 100a)
白色化フィルタ計算部(11)は、各マイクロホンで得られた観測信号全体 {x (t) ; 1 m
≤t≤N}を白色化するフィルタ(白色化フィルタ)の係数 {f (k) ;0≤k≤X}を X次の 線形予測分析によって計算する。この計算は線形予測分析と同じであり上記参考文 献 1を参照されたい。白色化フィルタの係数は、白色化フィルタ適用部(12)の入力と なる。
[0094] (ステップ S 100b)
白色化フィルタ適用部(12)は、式(39)に従って、上記白色化フィルタを各マイクロ ホンで得られた観測信号に適用して、白色化信号 w (t)を得る。既述のとおり、式(3 1)は式 (40)に変更すればよいので、第 1実施形態において、逆フィルタ計算部(13 )、とくに第 1の予測誤差フィルタ適用部(1311)による処理を式 (31)ではなく式 (40 )による計算処理に改めればよい。また、第 1実施形態において、逆フィルタ適用部( 14)による処理を、式 (42)ではなく式 (44)による計算処理に改めればよい。ステップ SlOObの処理の後、第 1実施形態のステップ S100〜S104の処理を行う力 これら の処理では第 1実施形態の各処理における観測信号をステップ SlOObの処理で得 られた白色化信号に読み替えて第 1実施形態と同様の処理を行う。この意味で、図 7 では、第 1実施形態のステップ S100〜S104の各処理に相当する処理を示す符号 に記号' を付している。
[数 36]
M L
y(t) = ∑ ∑gm (k) wm (t - k) ( 4 4 )
m=\ =0
[0095] <実施例 1 >
発明者らは第 2実施形態の実証実験を行ったので、その実験結果を示す。実験条 件として、マイクロホンの数 M=4、白色化フィルタの次数 X= 500、逆フィルタの次 数 L= 1000、窓関数の切出しサンプル数(1フレームのサンプル数) W= 200、予測 誤差フィルタの次数 P= 16、繰り返し回数 R = 10、逆フィルタ計算部の更新回数 R
1 2
= 20とした。学習率 7? (u)は、初期値を 0. 05に設定し、もし式(27)によって式(26) の値が減少するならば、 r? (u)の値を順次半減させることで、式(26)の値が必ず増 大するようにした。図 6に示した逆フィルタ適用部(14)へ入力する初期逆フィルタは、
式 (45)のように設定した。
[数 37] l/ = 200
gm(k) = \ ユ ,、, ^ , , ≤m≤M ( 4 5 )
S m [0 それ以外のとき
[0096] 本発明の第 2実施形態の効果を、信号歪み除去の指標として D 値 (インパルス応
50
答の全エネルギーに対する初期の 50msecまでのエネルギーの比)を用いて評価し た。連続発話データベース力も男女各一名の発話を取り出し、残響時間 0. 5秒の残 響室で測定したインパルス応答を畳み込むことで観測信号を合成した。
図 8は、男声および女声について観測信号長 Nを 5秒、 10秒、 20秒、 1分、 3分に 変化させたときの、繰り返し回数 R (各信号長 Nを持つ観測信号を図 6に示す逆フィ ルタ適用部(14)と、予測誤差フィルタ計算部(15)と、逆フィルタ計算部(13)を一巡 する処理を実行して逆フィルタを求める回数)と D 値の関係を示している。いずれの
50
場合においても、繰り返し回数を増加させると D 値が向上しており、繰り返し処理の
50
効果が顕著に見て取れる。特に観測信号長が 5〜10秒程度の比較的短い長さであ つても、繰り返し処理によって D 値が大きく向上したことが分かる。
50
[0097] また、本発明の第 2実施形態の効果を、音声スぺタトグラムの比較力も検証した。
図 9Aは観測信号長が 1分の場合の残響を含まない音声 (原音声)、図 9Bは観測 信号長が 1分の場合の残響を含む音声 (観測音声)、図 9Cは観測信号長が 1分の場 合の残響除去後の音声 (復元音声)のスペクトログラムの抜粋を示して 、る。図 9Aと 図 9Cとの対比および図 9Bと図 9Cとの対比から、観測信号に含まれる残響が抑制さ れ、原音声固有の特徴である調波構造やフォルマント構造が回復されたことが分力る
[0098] また、本発明の第 2実施形態の効果を、 LPCスペクトル歪みを用いて検証した。
図 10Bは原音声の波形、図 10Aは原音声と観測音声との LPCスペクトル歪みの時 系列(図中の点線)および原音声と復元音声との LPCスペクトル歪みの時系列(図中 の実線)を示している。図 10Aおよび図 10Bの各横軸は秒単位の時間を表し、両図 でタイムスケールを揃えている。図 10Bの縦軸は振幅値を表している。但し、原信号
の相対的な大きさが分かればよいので、この趣旨から図 10Bの縦軸では単位を明示 していない。図 10Aの縦軸は、 LPCスペクトル歪み SD (dB)を表している。
図 10Aから、原音声と復元音声との LPCスペクトル歪みの時系列(図中の実線)は 、原音声と観測音声との LPCスぺクトル歪みの時系列(図中の点線)よりも常に小さ いことがわかる。なお、観測音声では LPCスペクトル歪みの平均が 5. 39dB、分散が 4. 20dBであったのに対して、復元音声では平均が 2. 38dB、分散が 2. OOdBであ つた o
また、図 10Aと図 10Bとの対比から、原音声と復元音声との LPCスペクトル歪みの 時系列(図中の実線)が大きな値を示す区間 (例えば約 1. 0秒〜約 1. 2秒の区間を 参照)は、原音声の波形の振幅値がほぼ 0であることがわかる。実際、この区間では 発声がなく無音区間である。このため、実際に知覚される歪みはかなり小さくなつてい た。つまり、発声区間における原音声と復元音声との LPCスペクトル歪みの時系列( 図中の実線)は、原音声と観測音声との LPCスペクトル歪みの時系列(図中の点線) よりもかなり小さぐこのため原音声のスペクトルを高い精度で復元できたことが結論 付けられる。
[0099] <第 3実施形態 >
第 3実施形態は、第 1実施形態の変形例に相当する。具体的には、 § 2で述べた二 次統計量に基づく信号歪み除去処理を行なう形態である。そこで、第 1実施形態と異 なる部分について図 11および図 12を参照して説明を加える。但し、第 3実施形態で は、マイクロホンの個数 Mは 2以上とする。
[0100] ステップ S100の処理およびステップ S101の処理は、第 1実施形態と同じである。
[0101] ステップ S101の処理に続いて、ステップ S 102aの処理を行う。
第 3実施形態に係る逆フィルタ計算部(13)の機能構成例を、図 11を参照して説明 する。
逆フィルタ計算部(13)は、観測信号への予測誤差フィルタ適用部として機能する 第 1の予測誤差フィルタ適用部(1311)と、観測信号に更新用逆フィルタを適用して 得られる信号 (更新用逆フィルタ適用後信号)への予測誤差フィルタ適用部として機 能する第 2の予測誤差フィルタ適用部(1312)と、勾配べ外ル計算部(1313)と、逆
フィルタ更新部(132)および更新用逆フィルタ適用部(133)によって構成される。こ こで更新用逆フィルタは、式(37)の g (k)に相当する。
[0102] 第 1の予測誤差フィルタ適用部(1311)は、 m番目〔l≤m≤M〕のマイクロホンで観 測された観測信号 X (t)をフレーム化して、各フレームにっき、 i番目のフレームの信 号 X (n)に対してステップ S 101の処理で得られた i番目の予測誤差フィルタ a (k)を 適用して予測誤差フィルタ適用後の信号 v (n)を計算する (式 (38)を参照)。具体 的には、フレーム化処理部 (402B)力 入力された観測信号 X (t)に対してフレーム 化処理を行 ヽ、観測信号 X (t)の i番目のフレームの信号 X (n)を出力する。そして
、予測誤差フィルタ適用部 (404i)が信号 x (n)を入力として、式 (38)に従って信号
V (n)を出力する。但し、 l≤i≤Fである。
mi
[0103] 第 2の予測誤差フィルタ適用部(1312)は、更新用逆フィルタ適用後信号 y(t)をフ レーム化して、各フレームにっき、各フレームに対してステップ S101の処理で得られ た i番目の予測誤差フィルタ a (k)を適用して予測誤差フィルタ適用後のイノべーショ ン推定値 d (1) , · ··, d (W)を計算する (式 (30)を参照)。なお、更新用逆フィルタ適 用後信号 y (t)の初期値は、ステップ S 100の処理で得られた信号とすればよい。具 体的には、フレーム化処理部 (402A)が、初期値の場合を除き、後述する更新用逆 フィルタ適用部(133)が出力した更新用逆フィルタ適用後信号 y(t)に対してフレー ム化処理を行い、更新用逆フィルタ適用後信号 y(t)をフレーム化して、 i番目のフレ ームの信号 y (n)を出力する。そして、予測誤差フィルタ適用部 (403i)が信号 y (n) を入力として、式(30)に従ってイノベーション推定値 d (1) , · ··, d (W)を出力する。 但し、 l≤i≤Fである。
[0104] 勾配ベクトル計算部(1313)は、信号 V (n)とイノベーション推定値 d (n)とを用い て現在の更新用逆フィルタ g (k)の勾配ベクトルを計算する(式(37)の右辺第二項 を参照)。具体的には、各フレーム番号 i (l≤i≤F)に関して、相互相関計算部 (405 i)は信号 V (n)とイノベーション推定値 d (n)との相互相関く d (n)v (n—k)〉 wを
=l 計算する。また、各フレーム番号 i (l≤i≤F)に関して、分散計算部 (406i)は、イノべ ーシヨン推定値 d (1) , · ··, d (W)の分散く (11) 2〉 wを求める。各フレーム番号 i(l ≤i≤F)に関して、除算部 (407i)は、〈d (n)v (n-k) > V<d (n) 2> wを求め
る。加算部 (407)は、除算部 (4071)〜(407F)の全フレームに亘る総和、つまり式( 37)の右辺第二項を求める。
[0105] 逆フィルタ更新部(132)は、現在の更新用逆フィルタ g (k)、学習率 δ、勾配べタト ルを入力として、式(37)に従って、 u+ l回目の更新用逆フィルタ g (k) ' を求める。 式 (37)は、求められた g (k) ' を新たな g (k)と見立てて更新を行なうことを意味す m m
る。
[0106] 更新用逆フィルタ適用部(133)は、逆フィルタ更新部(132)によって得られた g (k m 丫 、つまり新たな g (k)および観測信号 x (t)を入力として、式 (42)に従って、更新 用逆フィルタ適用後信号 y(t)を求める。具体的には、式 (42)の g (k)として u+ l回 目の更新で得られた gを用いて計算する。この計算で得られた更新用逆フィルタ適用 後信号 y (t)は、第 2の予測誤差フィルタ適用部( 1312)の入力となる。
[0107] ステップ S 102aの処理に続いて、ステップ S 103およびステップ S 104の処理を行う 力 第 1実施形態と同じであるから説明を略する。
[0108] <実施例 2>
発明者らは第 3実施形態の実証実験を行ったので、その実験結果を示す。実験条 件として、 M=4, L= 1000, W= 200, P= 16, R =6, R = 50とした。学習率 δ
1 2
は、初期値を 0. 05に設定し、∑ Flog< d (n) 2> wの値が増加するならば、学習 率 δの値を順次半減させることで∑ Flog< d (n) 2> wの値が必ず減少するよう に設定した。逆フィルタの初期推定値は、 g (k) =0, l≤m≤M, l≤k≤Lとして m 設 し 7こ。
[0109] 本発明の第 3実施形態の効果を、音声明瞭度を表す RASTI (参考文献 5を参照) を残響除去の指標として評価した。連続発話データベースから男女各五名の発話を 取り出し、残響時間 0. 5秒の残響室で測定したインパルス応答を畳み込むことで観 測信号を合成した。
考文献 5) H. kuttruff. Room acoustics. Elsevier Applied Science, tnird edition, P.237 1991.
[0110] 図 13は、観測信号長 Nを 3秒、 4秒、 5秒、 10秒としたときの、 RASTIの値を表示し たものである。図 13に示すように、観測信号が 3〜5秒のように短時間の場合でも、高
V、残響除去性能を示して 、ることが分かる。
図 14は、残響除去前後におけるエネルギー減衰曲線の例である。直接音が到達し てから 50ミリ秒後の反射音のエネルギーが 15dB低減されていることが分かる。
産業上の利用可能性
本発明は、様々な信号処理システムの性能向上に寄与する要素技術であるところ 、例えば音声認識システム、テレビ会議システム、補聴器、音楽情報処理システム等 に利用することができる。