(57)【要約】【課題】パワースペクトル微細構造も考慮するようにしてより高い雑音抑圧効果を発揮する雑音抑圧フィルタ算出方法を提供する。【解決手段】最大事後確率パラメータ推定過程とウィーナーフィルタ計算過程と状態事後確率算出過程と重みつき和計算過程を含む。最大事後確率パラメータ推定過程は、観測短時間スペクトルを入力として、混合型包絡事前分布の包絡状態と、取り得る全てのパワースペクトル微細構造が記録されている微細構造符号帳の微細構造符号の各組み合わせについて、観測短時間スペクトルに最も合致する全極パラメータの推定値である全極パラメータ状態依存最大事後確率推定値を計算する。重みつき和計算過程は、状態依存ウィーナーフィルタと結合状態事後確率とを乗じた値を、包絡状態と微細構造符号の全ての組み合わせについて加算して雑音抑圧フィルタを計算する。【選択図】図３
2011164126.tif 000002
【特許請求の範囲】
【請求項１】
観測短時間スペクトルと、外部から入力される雑音のパワースペクトルとを入力として、パワースペクトル包絡を表す全極パラメータの事前分布である混合型包絡事前分布の包絡状態と、取り得る全てのパワースペクトル微細構造と微細構造符号事前確率が記録されている微細構造符号帳の微細構造符号の各組み合わせについて、上記観測短時間スペクトルに最も合致する全極パラメータの推定値である全極パラメータ状態依存最大事後確率推定値を計算する最大事後確率パラメータ推定過程と、
上記雑音のパワースペクトルと、上記全極パラメータ状態依存最大事後確率推定値を入力として、上記包絡状態と上記微細構造符号の各組み合わせについて、状態依存ウィーナーフィルタを計算するウィーナーフィルタ計算過程と、
上記観測短時間スペクトルと、上記雑音のパワースペクトルとを入力として、上記包絡状態と上記微細構造符号の各組み合わせについて、結合状態事後確率を計算する状態事後確率算出過程と、
上記状態依存ウィーナーフィルタと上記結合状態事後確率とを乗じた値を、包絡状態と微細構造符号の全ての組み合わせについて加算して雑音抑圧フィルタを計算する重みつき和計算過程と、
を含む雑音抑圧フィルタ算出方法。
【請求項２】
請求項１に記載した雑音抑圧フィルタ算出方法において、
上記混合型包絡事前分布は、
その要素分布である各状態依存包絡事前分布が相異なる正規−逆ガンマ分布で規定され、
当該正規−逆ガンマ分布が第一信頼度、第二信頼度、線形予測係数平均ベクトル、線形予測係数共分散行列、残差パワースケールからなる全極超パラメータで規定される、
ことを特徴とする雑音抑圧フィルタ算出方法。
【請求項３】
請求項２に記載した雑音抑圧フィルタ算出方法において、
上記第一信頼度と上記第二信頼度は、予め定められた定数であり、
上記線形予測係数平均ベクトルは、予め与えられた包絡符号帳の符号語が表す線形予測係数ベクトルに設定され、
上記線形予測係数共分散行列は、上記包絡符号帳の符号語が表す線形予測係数ベクトルと残差パワーを用いて、逆フーリエ変換とテプリッツ行列化を介して設定され、
上記残差パワースケールは、上記包絡符号帳の符号語が表す残差パワーと第二信頼度の積に設定される、
ことを特徴とする雑音抑圧フィルタ算出方法。
【請求項４】
請求項２又は３に記載した雑音抑圧フィルタ算出方法において、
上記最大事後確率パラメータ推定過程は、
Ｅ−ステップとＭステップの二つのステップで上記包絡状態と上記微細構造符号の各組み合わせについての最大事後確率パラメータを推定する過程であり、
上記Ｅ−ステップは、上記観測短時間スペクトルと上記雑音のパワースペクトルと上記全極パラメータ状態依存最大事後確率推定値を入力として、クリーン音声信号の短時間スペクトルの最小平均二乗誤差推定値と最小平均二乗誤差を計算し、
上記Ｍ−ステップは、上記クリーン音声信号の短時間スペクトルの最小平均二乗誤差推定値と最小平均二乗誤差を入力として、上記状態依存全極パラメータ最大事後確率推定値を計算する、
ステップであることを特徴とする雑音抑圧フィルタ算出方法。
【請求項５】
請求項１乃至４の何れかに記載した雑音抑圧フィルタ算出方法において、
上記状態事後確率算出過程は、
包絡状態と微細構造符号の各組み合わせについて、上記観測短時間スペクトルと上記雑音のパワースペクトルと上記全極パラメータ状態依存最大事後確率推定値を用いて、上記結合状態事後確率を近似計算する過程であることを特徴とする雑音抑圧フィルタ作成方法。
【請求項６】
請求項１乃至５の何れかに記載した雑音抑圧フィルタ算出方法において、
上記微細構造符号帳は、
その符号語が、相異なる調波パラメータを持つ調波混合ガウスモデルが表すパワースペクトルで規定される、
ことを特徴とする雑音抑圧フィルタ算出方法。
【請求項７】
パワースペクトル包絡を表す全極パラメータの事前分布である混合型包絡事前分布の各包絡状態の超パラメータと包絡状態事前確率が記録されている混合型包絡事前分布超パラメータ記録部と、
取り得る全てのパワースペクトル微細構造と微細構造符号事前確率が記録されている微細構造符号帳と、
観測短時間スペクトルと、外部から入力される雑音のパワースペクトルと、上記混合型包絡事前分布の各包絡状態の超パラメータと、上記各パワースペクトル微細構造を入力として、上記包絡状態と微細構造符号の各組み合わせについて、上記観測短時間スペクトルに最も合致する全極パラメータの推定値である全極パラメータ状態依存最大事後確率推定値を計算する最大事後確率パラメータ推定部と、
上記雑音のパワースペクトルと、上記全極パラメータ状態依存最大事後確率推定値を入力として、上記包絡状態と上記微細構造符号の各組み合わせについて、状態依存ウィーナーフィルタを計算するウィーナーフィルタ計算部と、
上記観測短時間スペクトルと、上記雑音のパワースペクトルと、上記混合型包絡事前分布の各包絡状態の超パラメータと、上記各パワースペクトル微細構造を入力として、上記包絡状態と上記微細構造符号の各組み合わせについて、結合状態事後確率を計算する状態事後確率算出部と、
上記状態依存ウィーナーフィルタと、上記結合状態事後確率を入力として、上記状態依存ウィーナーフィルタと上記結合状態事後確率とを乗じた値を、包絡状態と微細構造符号の全ての組み合わせについて加算して雑音抑圧フィルタを計算する重みつき和計算部と、
を具備する雑音抑圧フィルタ算出装置。
【請求項８】
請求項７に記載した雑音抑圧フィルタ算出装置において、
上記混合型包絡事前分布超パラメータ記録部は、
その要素分布である各状態依存包絡事前分布が相異なる正規−逆ガンマ分布で規定され、
当該正規−逆ガンマ分布が第一信頼度、第二信頼度、線形予測係数平均ベクトル、線形予測係数共分散行列、残差パワースケールからなる全極超パラメータで規定される、
ことを特徴とする雑音抑圧フィルタ算出装置。
【請求項９】
請求項８に記載した雑音抑圧フィルタ算出装置において、
上記第一信頼度と上記第二信頼度は、予め定められた定数であり、
上記線形予測係数平均ベクトルは、予め与えられた包絡符号帳の符号語が表す線形予測係数ベクトルに設定され、
上記線形予測係数共分散行列は、上記包絡符号帳の符号語が表す線形予測係数ベクトルと残差パワーを用いて、逆フーリエ変換とテプリッツ行列化を介して設定され、
上記残差パワースケールは、上記包絡符号帳の符号語が表す残差パワーと第二信頼度の積に設定される、
ことを特徴とする雑音抑圧フィルタ算出装置。
【請求項１０】
請求項８又は９に記載した雑音抑圧フィルタ算出装置において、
上記最大事後確率パラメータ推定部は、
Ｅ−ステップ実行手段とＭ−ステップ実行手段を備え、上記包絡状態と上記微細構造符号の各組み合わせについての最大事後確率パラメータを推定するものであり、
上記Ｅ−ステップ実行手段は、上記観測短時間スペクトルと上記雑音のパワースペクトルと上記全極パラメータ状態依存最大事後確率推定値を入力として、クリーン音声信号の短時間スペクトルの最小平均二乗誤差推定値と最小平均二乗誤差を出力し、
上記Ｍ−ステップ実行手段は、上記クリーン音声信号の短時間スペクトルの最小平均二乗誤差推定値と最小平均二乗誤差を入力として、上記状態依存全極パラメータ最大事後確率推定値を出力する、
ものであることを特徴とする雑音抑圧フィルタ算出装置。
【請求項１１】
請求項７乃至１０の何れかに記載した雑音抑圧フィルタ算出装置において、
上記状態事後確率算出部は、
包絡状態と微細構造符号の各組み合わせについて、上記観測短時間スペクトルと上記雑音のパワースペクトルと上記全極パラメータ状態依存最大事後確率推定値を用いて、上記結合状態事後確率を近似計算する、
ものであることを特徴とする雑音抑圧フィルタ作成装置。
【請求項１２】
請求項７乃至１１の何れかに記載した雑音抑圧フィルタ算出装置において、
上記微細構造符号帳は、
その符号語が、相異なる調波パラメータを持つ調波混合ガウスモデルが表すパワースペクトルで規定される、
ことを特徴とする雑音抑圧フィルタ算出装置。
【請求項１３】
請求項７乃至１２のいずれか１項に記載の雑音抑圧フィルタ算出装置の各部の機能を、コンピュータに実行させるための装置プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、音声信号に含まれる雑音を抑圧する雑音抑圧フィルタ算出方法と、その装置と、プログラムに関する。
【背景技術】
【０００２】
図６に一般的な雑音抑圧方法の処理の流れを示す。最初に、短時間フーリエ変換過程において、雑音を含む観測音声信号ｘ（t）から短時間フーリエ係数ｘn,qを計算する（ステップＳ６０）。ｔは時間、nは短時間フレーム、qは周波数ビンのインデックスである。周波数ビンｑの総数をＱと表記する。
【０００３】
次に雑音パワースペクトル推定過程において、観測短時間スペクトルＸn＝{ｘn,q}0≦q≦Q-1に含まれる雑音のパワースペクトルの推定値{λVn（q）}0≦q≦Q-1を計算する（ステップＳ６１）。続いて、雑音抑圧フィルタ算出過程において、雑音抑圧フィルタ{Ｇn,q}0≦q≦Q-1を計算する（ステップＳ６２）。フィルタ適用過程では、雑音抑圧フィルタ{Ｇn,q}0≦q≦Q-1を観測短時間スペクトルに適用して出力音声信号の短時間スペクトル{s^n,q}0≦q≦Q-1を計算する（式（１））（ステップＳ６３）。なお、＾等は式中の表記が正しい。
【０００４】
【数１】
2011164126.tif 000003
【０００５】
そして最後に、逆短時間フーリエ変換過程において、短時間スペクトルs^n,qから出力音声信号ｓ＾（t）を計算する（ステップＳ６４）。
【０００６】
雑音抑圧フィルタ算出過程のより具体的な動作を、非特許文献１に開示された例で説明する。その方法は、音声信号を線形予測分析して得られる線形予測係数のベクトルと残差パワー（これらをまとめて全極パラメータと称する）に関する大きさＭＥの符号帳（これを包絡符号帳と称する）を用いる。符号i（1≦ｉ≦MＥ）に対応する符号語は、全極パラメータ{ａi，σ２i}である。ただし、ａi＝[ａi，1，…，ai，P]Ｔとσ２iは、それぞれ線形予測係数{ai，ｋ}1≦ｋ≦Pのベクトルと残差パワーである。この符号帳を用いて雑音抑圧フィルタ{Ｇn,q}0≦ｑ≦Ｑ-1を式（２）にしたがって計算する。
【０００７】
【数２】
2011164126.tif 000004
【０００８】
ただし、ｆWiener（λ1，λ2）は式（３）で与えられるウィーナーフィルタである。
【０００９】
【数３】
2011164126.tif 000005
【００１０】
式（２）のλAP(q;a,σ2)は、線形予測係数ベクトルａ＝[ａ1，…，ａP]Ｔと残差パワーσ2が規定するパワースペクトルを表す（式（４））。線形予測分析の次数Pは通常小さな値に設定されるため、このパワースペクトルは微細な凹凸を持たない滑らかなものである。
【００１１】
【数４】
2011164126.tif 000006
【００１２】
また、p（i｜Xn）は観測短時間スペクトルがＸnであるときの符号ｉの事後確率を表し、式（５）で計算される。
【００１３】
【数５】
2011164126.tif 000007
【００１４】
ただし、ｐ（Xn｜i）は符号iから観測短時間スペクトルＸnが生起する確率密度関数を表し、式（６）で定義される。
【００１５】
【数６】
2011164126.tif 000008
【００１６】
ｆNC(ｘ；μ，σ2)は、式（７）で定義される平均μ、分散σ2の複素正規分布の確率密度関数である。
【００１７】
【数７】
2011164126.tif 000009
【００１８】
以上述べた計算を、図７に示す動作フローで説明する。まず、包絡符号帳を探索するインデックスｉが１に初期化される。次に、ウィーナーフィルタ算出過程において、全極パラメータ{ａi，σ２i}を用いて、ウィーナーフィルタの式（２）の右辺第２項fWiener(q;λAP(ai,σ2i)λVn(q))を計算する（ステップＳ７０）。そして、符号重み算出過程において、符号ｉの事後確率ｐ（i｜Xn）を計算する（ステップＳ７１）。もし、i＝MＥでなければ、ｉを１だけ増やして次の符号に対して同様の処理を行う。i＝MＥであれば、重みつき和計算過程において式（２）の計算を行う（ステップＳ７３）。
【先行技術文献】
【非特許文献】
【００１９】
【非特許文献１】Y.Ephraim, “A Bayesian estimation approach for speech enhancement using hidden Markov modeles,” IEEE Trans. Signal Process., vol. 40, no.4, pp.725-735, 1992.
【発明の概要】
【発明が解決しようとする課題】
【００２０】
一般的に、音声信号のパワースペクトルは包絡と微細構造とに分解することができる。従来技術は、線形予測係数と残差パワー（すなわち、全極パラメータ）に関する符号帳を用いるため、パワースペクトル包絡の回復に高い効果を持つ。しかし、パワースペクトルの微細構造は回復されない。
【００２１】
図１に音声信号の有声部におけるパワースペクトルを概念的に示す。図１の横軸は周波数（Hz）、縦軸はパワースペクトルの大きさ（dB）である。細い実線で示すパワースペクトル包絡はパワースペクトルの滑らかな傾斜であり、主に音韻を表現する。破線で示すパワースペクトル微細構造はパワースペクトルの微細な凹凸を表し、大局的には平坦である。パワースペクトル微細構造は主に声の高さを表現する。従来の雑音抑圧フィルタ算出方法で求めた雑音抑圧フィルタを用いると、このパワースペクトル包絡だけが用いられる。
【００２２】
しかし、音声は有声部においては調波構造を持つため、倍音周波数近傍にエネルギーが集中する微細構造を持つことが知られており、実際の音声信号のパワースペクトルは、パワースペクトル包絡にパワースペクトル微細構造を重畳した太い実線で示す形状をもつ。このようなスペクトル微細構造に着目することで、より高い雑音抑圧効果を発揮できることが期待される。
【００２３】
この発明は、このような考察に基づいてなされたものであり、従来から用いられているパワースペクトル包絡だけでなく、パワースペクトル微細構造も考慮するようにしてより高い雑音抑圧効果を発揮できる雑音抑圧フィルタ算出方法と、その装置と、プログラムを提供することを目的とする。
【課題を解決するための手段】
【００２４】
この発明の雑音抑圧フィルタ算出方法は、最大事後確率パラメータ推定過程と、ウィーナーフィルタ計算過程と、状態事後確率算出過程と、重みつき和計算過程と、を含み、パワースペクトル包絡を表す全極パラメータの事前分布である混合型包絡事前分布（混合型包絡事前分布は複数の包絡状態から成り立つ）と、取り得る全てのパワースペクトル微細構造が記録されている微細構造符号帳と、を用いる。最大事後確率パラメータ推定過程は、観測短時間スペクトルと外部から入力される雑音のパワースペクトルを入力として、混合型包絡事前分布の包絡状態と、微細構造符号帳の微細構造符号の各組み合わせについて、観測短時間スペクトルに最も合致する全極パラメータの推定値である全極パラメータ状態依存最大事後確率推定値を計算する。ウィーナーフィルタ計算過程は、雑音のパワースペクトルと、全極パラメータ状態依存最大事後確率推定値を入力として、包絡状態と微細構造符号の各組み合わせについて、状態依存ウィーナーフィルタを計算する。状態事後確率算出過程は、観測短時間スペクトルと雑音のパワースペクトルを入力として、包絡状態と微細構造符号の各組み合わせについて、結合状態事後確率を計算する。重みつき和計算過程は、状態依存ウィーナーフィルタと結合状態事後確率とを乗じた値を、包絡状態と微細構造符号の全ての組み合わせについて加算して雑音抑圧フィルタを計算する。
【発明の効果】
【００２５】
この発明の雑音抑圧フィルタ算出方法によれば、パワースペクトルの微細構造に着目することで、有声部の倍音周波数近傍にエネルギーが集中するという音声の微細構造の特徴も回復できる。したがって、この発明の雑音抑圧フィルタ算出方法は、より高い雑音抑圧効果を発揮する雑音抑圧フィルタの算出を可能にする。
【図面の簡単な説明】
【００２６】
【図１】音声信号のパワースペクトルの包絡と微細構造を概念的に示す図。
【図２】この発明の雑音抑圧フィルタ算出装置１００の機能構成例を示す図。
【図３】雑音抑圧フィルタ算出装置１００の動作フローを示す図。
【図４】実施例２の最大事後確率パラメータ推定過程の動作フローを示す図。
【図５】評価実験の結果を示す図。
【図６】従来の一般的な雑音抑圧方法の処理の流れを示す図。
【図７】従来の雑音抑圧フィルタ算出方法の処理の流れを示す図。
【発明を実施するための形態】
【００２７】
以下、この発明の実施の形態について、図面を参照して説明する。複数の図面中同一のものには同じ参照符号を付し、説明は繰り返さない。
【実施例１】
【００２８】
〔主要部分の説明〕
図２にこの発明の雑音抑圧フィルタ算出装置１００の機能構成例を示す。図３にその動作フローを示す。雑音抑圧フィルタ算出装置１００は、混合型包絡事前分布超パラメータ記録部１０と、微細構造符号帳１１と、最大事後確率パラメータ推定部１２と、ウィーナーフィルタ計算部１３と、状態事後確率算出部１４と、重みつき和計算部１５と、制御部１６と、を具備する。その各部の機能は、例えばＲＯＭ、ＲＡＭ、ＣＰＵ等で構成されるコンピュータに所定のプログラムが読み込まれて、ＣＰＵがそのプログラムを実行するこ９とで実現することができる。
【００２９】
この発明の新しい点は２つある。その１つは、パワースペクトル微細構造に関する予め作成された符号帳を用いる点である。この符号帳が微細構造符号帳１１である。微細構造符号帳１１の大きさをMＦとすると、符号ｊ（1≦ｊ≦MＦ）には、パワースペクトル微細構造{λＦ（q；j）}0≦ｑ≦Ｑ-1が対応付けられており、音声のパワースペクトル微細構造は{λＦ（q；1）}0≦ｑ≦Ｑ-1，…，{λＦ（q；MF）}0≦ｑ≦Ｑ-1の何れかしか取り得ないと仮定される。また、微細構造符号jが選択される事前確率p（j）も予め与えられており、微細構造符号帳１１に記録されている。
【００３０】
その２つは、パワースペクトル包絡、具体的には全極パラメータの事前知識を、符号帳ではなく混合連続確率分布の形式で表現することである。つまり、MＥ個の符号語の代わりにMＥ個の状態（包絡状態と称する）があって、各包絡状態ｉ（１≦ｉ≦MＥ）には、全極パラメータの事前分布ｐ（a,σ2｜i）と包絡状態iが選択される事前確率ｐ（i）が対応付けられている。よって、この発明によると、全極パラメータの事前分布は式（８）で与えられる。
【００３１】
【数８】
2011164126.tif 000010
【００３２】
式（８）で与えられるp(a,σ2)を混合型包絡事前分布と称し、p(a,σ2｜i)を状態依存包絡事前分布と称する。p(i)については、例えばp(i)＝1/MEと設定する方法がある。混合型包絡事前分布超パラメータ記録部１０には、p(i)とp(a,σ2｜i)を規定する超パラメータが記録されている。
【００３３】
最大事後確率パラメータ推定部１２は、観測短時間スペクトルXnと、外部から入力される雑音のパワースペクトル{λVn(q)} 0≦ｑ≦Ｑ-1と、各包絡状態iに対する状態依存包絡事前分布p(a,σ2｜i)の超パラメータと包絡状態事前確率ｐ(i)と、各微細構造符号jに対するパワースペクトル微細構造{λF(q;j)}0≦ｑ≦Ｑ-1と微細構造符号事前確率p(j)を入力として、包絡状態と微細構造符号の各組み合わせについて、観測短時間スペクトルXnに最も合致する全極パラメータの推定値である全極パラメータ状態依存最大事後確率推定値{a＾i,j,σ＾2 i,j }を計算する（ステップＳ１２）。
【００３４】
ウィーナーフィルタ計算部１３は、雑音のパワースペクトルλVn(q)と、全極パラメータ状態依存最大事後確率推定値{a＾i,j,σ＾2 i,j }を入力として、包絡状態と微細構造符号の各組み合わせについて、状態依存ウィーナーフィルタｆWiener(q; λS（a＾i,j,σ＾2 i,j,j）λVn(q))を計算する（ステップＳ１３）。
【００３５】
状態事後確率算出部１４は、観測短時間スペクトルXnと、雑音のパワースペクトル{λVn(q)} 0≦ｑ≦Ｑ-1と、各包絡状態iに対する状態依存包絡事前分布p(a,σ2｜i)の超パラメータと包絡状態事前確率p(i)と、各微細構造符号jに対するパワースペクトル微細構造{λF(q;j)}0≦ｑ≦Ｑ-1と微細構造符号事前確率p(j)を入力として、包絡状態と微細構造符号の各組み合わせについて、結合状態事後確率p（i,j｜Xn）を計算する（ステップＳ１４）。重みつき和計算部１５は、状態依存ウィーナーフィルタと結合状態事後確率p（i,j｜Xn）とを乗じた値を、包絡状態と微細構造符号の全ての組み合わせについて加算して雑音抑圧フィルタ{Gn,q}0≦ｑ≦Ｑ-1を計算する（ステップＳ１５）。
【００３６】
次に具体的な処理の内容を計算式を示して説明する。
【００３７】
この発明では、包絡符号帳を用いる場合と異なり、包絡状態に対して全極パラメータの値は予め一意に定められていない。そのため、先ず、最大事後確率パラメータ推定部１２で、包絡状態iと微細構造符号jの各組み合わせについて、最大事後確率推定により全極パラメータの状態依存最大事後確率推定値{a＾i,j,σ＾2 i,j }を式（９）で計算する（ステップＳ１２）。
【００３８】
【数９】
2011164126.tif 000011
【００３９】
ただし、全極パラメータの事後確率密度関数p（a,σ2｜Xn,i,j）は次式で与えられる。
【００４０】
【数１０】
2011164126.tif 000012
【００４１】
式（１０）に現れるp(a,σ2｜i)は、上記した状態依存包絡事前分布である。一方、式（１０）に現れる尤度関数p（Xn,｜a,σ2,j）は式（１１）によって与えられる。
【００４２】
【数１１】
2011164126.tif 000013
【００４３】
ただし、パワースペクトル{λS(q,a,σ2,j)}0≦ｑ≦Ｑ-1は、式（１２）に示すように全極パラメータ{ a,σ2}で決まるパワースペクトル包絡と微細構造符号jで決まるパワースペクトル微細構造の積で与えられる。
【００４４】
【数１２】
2011164126.tif 000014
【００４５】
この実施例では、重みつき和計算部１５が、こうして求められる全極パラメータ{a＾i,j,σ＾2 i,j }を用いて雑音抑圧フィルタ{Gn,q}0≦ｑ≦Ｑ-1を、式（１３）にしたがって計算する。
【００４６】
【数１３】
2011164126.tif 000015
【００４７】
ただし、ｆWiener(q;λS(a^i,j,σ^2i,j,j)λVn(q))は包絡状態iと微細構造符号jに依存して決まる状態依存ウィーナーフィルタである。また、p（i,j｜Xn）は包絡状態iと微細構造符号jの組み合わせに対して定義される結合状態事後確率である。これらはそれぞれ、ウィーナーフィルタ計算部１３と状態事後確率計算部１４で計算される。状態事後確率計算部１４は、p（i,j｜Xn）を式（１４）にしたがって計算する。
【００４８】
【数１４】
2011164126.tif 000016
【００４９】
なお、RPはP次元ユークリッド空間を表す。
【００５０】
以上、この発明の主要部である最大事後確率パラメータ推定部１２、状態事後確率算出部１４、ウィーナーフィルタ計算部１３、重みつき和計算部１５で処理される計算を説明した。これらの計算過程によって、図３に示す雑音抑圧フィルタ算出過程が構成される。
【００５１】
雑音抑圧フィルタ算出方法は、最大事後確率パラメータ推定過程（ステップＳ１２）と、ウィーナーフィルタ算出過程（ステップＳ１３）と、状態事後確率算出過程（ステップＳ１４）と、重みつき和計算過程（ステップＳ１５）とが、２重ループを含むように構成される。
【００５２】
先ず、混合型包絡事前分布を探索するインデックスiと微細構造符号帳を探索するインデックスjが、各々１に初期化される（ステップＳ１６０）。このようなインデックスの初期化と更新及び各部の制御は制御部１６が行う。最大事後確率パラメータ推定部１２は、全極パラメータの状態依存最大事後確率推定値{a＾i,j,σ＾2 i,j }を式（９）にしたがって計算する（ステップＳ１２）。ウィーナーフィルタ計算部１３は、外部から入力される雑音のパワースペクトル{λVn(q)} 0≦ｑ≦Ｑ-1と全極パラメータの状態依存最大事後確率推定値{a＾i,j,σ＾2 i,j }を入力として、状態依存ウィーナーフィルタｆWiener(q;λS(a^i,j,σ^2i,j,j)λVn(q))を式（１３）の右辺第２項で計算する（ステップＳ１３）。そして、状態事後確率算出部１４は、包絡状態iと微細構造符号jの結合状態事後確率p（i,j｜Xn）を式（１４）にしたがって計算する（ステップS１４）。以上の計算は、ｉ＝MＥで、かつj＝MＦになるまでｉとjの何れかを１ずつ増やして同様に計算される（ステップＳ１６１，Ｓ１６２）。重みつき和計算部１５は、各包絡状態iと各微細構造jの組み合わせについて計算された状態依存ウィーナーフィルタｆWiener(q;λS(a^i,j,σ^2i,j,j)λVn(q))と結合状態事後確率p（i,j｜Xn）とを乗じ、その全てを加算して雑音抑圧フィルタ{Gn,q}0≦ｑ≦Ｑ-1を計算する（式（１３）、ステップＳ１５）。
【００５３】
〔微細構造符号帳〕
この実施例では、微細構造符号帳１１を、その符号語が、相異なる調波パラメータを持つ調波混合ガウスモデルでモデル化されたパワースペクトルで規定する。調波混合ガウスモデルは次のように定義される。
【００５４】
モデルは２つのパラメータηとμを持つ。ηを発声指標と称し、ηは０か１の何れかの値を持ち、η＝１は有声音を、η＝０は無声音を表す。μは基本周波数（単位はラジアン）を表す。これら２つのパラメータをまとめて調波パラメータと称する。調波混合ガウスモデルは、調波パラメータを用いて式（１５）に示すように与えられる。
【００５５】
【数１５】
2011164126.tif 000017
【００５６】
ここでξは予め定められた定数であり、調波成分のスペクトルの広がりを表す。また、ｋはパワースペクトル微細構造のピーク（すなわち、高調波成分）のインデックスであり、そのピークは基本周波数μの間隔で発生する（図１参照）。なお、式（１５）の無限和は、実際には適当な有限の範囲で打ち切って計算する。
【００５７】
この調波混合ガウスモデルを用いて、微細構造符号帳を以下のように設定する。各微細構造符号j（1≦j≦MF）について、符号語として調波パラメータ{ηj, μｊ}を対応付ける。すなわち、微細構造符号jに対応するパワースペクトルは式（１６）で与えられる。
【００５８】
【数１６】
2011164126.tif 000018
【００５９】
MF個の符号語のうち、MF番目の符号語だけηj＝０（無声音）を満たし、残りのMF−１個の符号語についてはηj＝１（有声音）とする。MF−１個の有声音に対応する符号語について、基本周波数μｊが互いに異なるように設定する。
【００６０】
例えば、人間の声の基本周波数が概ね９０Ｈｚから３００Ｈｚ程度の範囲にあることを考慮して、この範囲で９０Ｈｚ，９２．６Ｈｚ，９５．４Ｈｚ，９８．１Ｈｚ，１０１．０Ｈｚ，１０４．０Ｈｚ，…のように、１/２４オクターブずつ対数軸上で等間隔にずらした周波数値をラジアンに変換したものをμ1，μ2，μ3，μ4，μ5，μ6，…とする方法がある。なお、符号事前確率p(j)は、p(j)＝1/MFとする。
【００６１】
〔混合型包絡事前分布〕
この実施例では、混合型包絡事前分布の各要素分布である状態依存包絡事前分布を相異なる正規−逆ガンマ分布で規定する。その正規−逆ガンマ分布は、第一信頼度、第二信頼度、線形予測係数平均ベクトル、線形予測係数共分散行列、残差パワースケールからなる全極超パラメータで規定される。これら全極超パラメータは、混合型包絡事前分布超パラメータ記録部１０に記録されている。その具体的な定義を以下に記述する。
【００６２】
各状態依存包絡事前分布p(a,σ2｜i)の関数形を式（１７）に示すように、正規−逆ガンマ分布を用いて定義する。
【００６３】
【数１７】
2011164126.tif 000019
【００６４】
ただし、ｆN(x;μ,Σ)は平均μ、共分散行列Σの多変量正規分布の確率密度関数（式（１８））である。ｆΓ-1(x;ρ,γ)は、形状パラメータρ、スケールパラメータγの逆ガンマ分布の確率密度関数（式（１９））である。
【００６５】
【数１８】
2011164126.tif 000020
【００６６】
ここで、Γ（・）はガンマ関数である。ζi,Θi,vi,ρi,γiは予め包絡状態i（１≦ｉ≦ME）毎に定められた定数であり、これら５種類のパラメータをまとめて全極超パラメータと称する。
【００６７】
ζiとρiはそれぞれ、線形予測係数ベクトルの事前分布（式（１７）の右辺第一項ｆN（・））と、残差パワーの事前分布（式（１７）の右辺第二項ｆΓ-1（・））の広がりの狭さを表す。確率密度関数が表す確率分布の広がりが狭いことは、その確率密度関数で与えられる情報の信頼度が高いことを意味するので、ζiとρiをそれぞれ第一信頼度、第二信頼度と称する。ｖiは線形予測係数ベクトルの事前分布の平均なので、これを線形予測係数平均ベクトルと称する。一方、Θiは線形予測係数ベクトルの事前分布の共分散行列であるから、線形予測係数共分散行列と称する。γiは逆ガンマ分布のスケールパラメータに相当するため、残差パワースケールと称する。
【００６８】
〔混合型包絡事前分布超パラメータ記録部〕
この実施例における混合型包絡事前分布超パラメータ記録１０を説明する。各包絡状態の状態依存包絡事前分布の全極超パラメータを、予め与えられた包絡符号帳を基に定めることができる（包絡符号帳は公知の方法で作成できる）。その場合、第一信頼度と第二信頼度は、予め定められた定数であり、線形予測係数平均ベクトルは、包絡符号帳の符号語が表す線形予測係数ベクトルに設定され、その線形予測係数ベクトルの共分散行列は、包絡符号帳の符号語が表す線形予測係数ベクトルと残差パワーを用いて、逆フーリエ変換とテプリッツ行列を介して設定され、残差パワースケールは、包絡符号帳の符号語が表す残差パワーと第二信頼度の積に設定される。なお、上記した通り、状態事前確率p(i)はp(i)＝1/MEと設定される。
【００６９】
具体的に説明する。先ず、全ての包絡状態i（１≦ｉ≦MF）について、第一信頼度ζiと第二信頼度ρiを予め定められた定数αに設定する。
【００７０】
【数１９】
2011164126.tif 000021
【００７１】
α→∞とすると、混合型包絡事前分布は、予め与えられた包絡符号帳と等価になる。
【００７２】
残りの全極パラメータは、包絡状態ごとに異なる値に設定される。包絡状態ｉの線形予測係数平均ベクトルは、予め与えられた包絡符号帳の包絡符号ｉの線形予測係数ベクトルａiに設定される。
【００７３】
【数２０】
2011164126.tif 000022
【００７４】
包絡状態ｉの残差パワースケールは、包絡符号ｉの残差パワーσ2iと第二信頼度ρiの積に設定される。
【００７５】
【数２１】
2011164126.tif 000023
【００７６】
最後に、線形予測係数共分散行列は次のようにして設定される。先ず、ｉ番目の包絡符号語{ai,σ2i}が表すパワースペクトル{λAP(q;ai,σ2i)}0≦ｑ≦Ｑ-1を計算する。このパワースペクトルに対して逆フーリエ変換することで、包絡符号語{ai,σ2i}に対応する自己相関係数列を求める。
【００７７】
【数２２】
2011164126.tif 000024
【００７８】
ただし、ｒi（ｔ）はラグｔの自己相関係数を表し、ＩＦＦＴ（・）は与えられたスペクトルをπ≦ω≦２πの部分に折り返してから逆フーリエ変換を計算する処理を表す。ラグP−１までの自己相関係数列からなるテプリッツ（Toplitz）行列を線形予測係数共分散行列とする。
【実施例２】
【００７９】
全極パラメータ状態依存最大事後確率推定値を計算する最大事後確率パラメータ推定過程の実施方法の一例を説明する。この方法では、混合型包絡事前分布を正規−逆ガンマ分布を用いて定義した場合において、未知のクリーン音声信号の短時間スペクトルＳn＝{sn,q}0≦ｑ≦Ｑ-1を潜在変数とするＥＭアルゴリズムを用いて、最大事後確率推定により全極パラメータを推定する。
【００８０】
その最大事後確率パラメータ推定過程は、Ｅ−ステップとＭステップの二つのステップで構成される。Ｅ−ステップは、観測短時間スペクトルと外部から入力される雑音のパワースペクトルと全極パラメータ状態依存最大事後確率推定値を入力として、クリーン音声信号の短時間スペクトルの最小平均二乗誤差推定値と最小平均二乗誤差を出力する。
【００８１】
Ｍ−ステップは、上記クリーン音声信号の短時間スペクトルの最小平均二乗誤差推定値と最小平均二乗誤差を入力として、上記全極パラメータ状態依存最大事後確率推定値を出力する。Ｅ−ステップはＥ−ステップ実行手段、Ｍ−ステップはＭ−ステップ実行手段で実行される。その動作について、図４を参照して説明する。
【００８２】
ＥＭアルゴリズムは繰り返し処理に基づく方法である。先ず、初期化ステップにおいて、全極パラメータ状態依存最大事後確率推定値{a^i,j,σ^2i,j}の初期値を定める（ステップＳ４０）。初期値の定め方としては、例えば、次のように混合型包絡事前分布の情報を用いる方法などが考えられる。  
【００８３】
【数２３】
2011164126.tif 000025
【００８４】
次に繰り返しのインデックスｋを１に初期化する。Ｅ−ステップでは、クリーン音声信号の短時間スペクトルの最小平均二乗誤差推定値s^n,q,i,jと最小平均二乗誤差εn,q,i,jを計算する（ステップＳ４１）。
【００８５】
【数２４】
2011164126.tif 000026
【００８６】
Ｍステップでは、全極パラメータ状態依存最大事後確率推定値を、次の最適化問題を解くことで更新する（ステップＳ４２）。
【００８７】
【数２５】
2011164126.tif 000027
【００８８】
ただし、Sn＝{ｓn,q}0≦ｑ≦Ｑ-1、CQはQ次元複素ベクトル空間を表し、ｗ（Sn）は式（３０）で与えられる確率密度関数である。
【００８９】
【数２６】
2011164126.tif 000028
【００９０】
式（２９）は、次の手順で計算される。先ず、エネルギースペクトル{Un,q,i,j}0≦ｑ≦Ｑ-1を式（３１）で計算する。
【００９１】
【数２７】
2011164126.tif 000029
【００９２】
定義から分かるように、{Un,q,i,j}0≦ｑ≦Ｑ-1は、パワースペクトル微細構造{λF（q,j）}0≦ｑ≦Ｑ-1の影響を取り除いたエネルギースペクトルである。このエネルギースペクトルに対して逆フーリエ変換して自己相関係数列を求める。
【００９３】
【数２８】
2011164126.tif 000030
【００９４】
こうして求められた自己相関係数列から、自己相関行列と自己相関ベクトルを構成する。
【００９５】
【数２９】
2011164126.tif 000031
【００９６】
Rn,i,jとrn,i,jを、状態依存包絡事前分布によって修正した修正自己相関行列と修正自己相関ベクトルを次式にしたがって計算する。
【００９７】
【数３０】
2011164126.tif 000032
【００９８】
更新後の全極パラメータ状態依存最大事後確率推定値は、次の線形方程式を解くことで求められる。
【００９９】
【数３１】
2011164126.tif 000033
【０１００】
ここで、０PはＰ次元零ベクトル、r〜n,i,j（０）はR〜n,i,jの対角成分である（R〜n,i,jはテプリッツ行列であるから、対角成分は全て等しい）。式（３７）はユール-ウォーカー（Yule-Walker）方程式の形式をしているため、レビンソン-ダービン（Levinson-Durbin）アルゴリズムによって高速に解ける。繰り返し回数ｋが予め定められた回数NRに達していれば終了し、そうでなければｋを１だけ増やしてＥ−ステップに戻る。
【実施例３】
【０１０１】
状態事後確率算出部１４で行う結合状態事後確率の計算（式（１４））を近似する方法を説明する。その算出過程は、包絡状態と微細構造符号の各組み合わせについて、観測短時間スペクトルと外部から入力される雑音のパワースペクトルだけでなく全極パラメータ状態依存最大事後確率推定値も用いて、結合状態事後確率を近似計算するものである。
【０１０２】
具体的には、式（１４）の右辺の分子を全極パラメータの状態依存最大事後確率推定値{a^i,j,σ＾2i,j}を用いて次式のように近似する。
【０１０３】
【数３２】
2011164126.tif 000034
【０１０４】
この近似を用いれば、式（１４）は次式によって計算される。
【０１０５】
【数３３】
2011164126.tif 000035
【０１０６】
以上説明した各実施例は、それぞれを組み合わせて実施することも可能である。
【０１０７】
〔評価実験〕
この発明による雑音抑圧フィルタ算出方法の効果を確認する目的で、実施例２で説明した雑音抑圧フィルタ算出方法を計算機上に実装して、Ａurora-２のテストセットＡを用いて評価実験を行った。実験条件としては、包絡状態数MEを６４、信頼度αの値を０．５，１．０，４．０の三種類に設定した。比較する従来技術については、包絡符号帳の大きさの設定を、ME＝６４とME＝５１２の２種類について実験した。
【０１０８】
図５に、この実験で得られた雑音抑圧処理前後における音声信号のケプストラム歪みを示す。図５の横方向は観測音声信号のＳＮ比を表し、縦方向は雑音抑圧に用いた方法である。入力音声信号は、観測信号でありフィルタに入力する前の未処理の信号である。
【０１０９】
平均ケプストラム歪みに着目するとこの発明の信頼度α＝０．５の時に最も良好（４．０３dB）な雑音抑圧効果を得た。この結果は、パワースペクトルの微細構造符号帳を用いることで、雑音抑圧効果が改善されることを示している。また、αの値を大きくすると雑音抑圧効果は低下した。式（２０）に示したように、αは、全極パラメータの事前分布の広がりの狭さ（すなわち、信頼度）を表す第一信頼度と第二信頼度の値である。式（２０）の直後で説明したように、α→∞とすると、混合型包絡事前分布は、包絡符号帳と等価になる。したがって、上記結果は、微細構造符号帳を用いる際には、全極パラメータの事前知識を、包絡符号帳ではなく混合型包絡事前分布で与える必要があることを示している。
【０１１０】
なお、上記装置における処理手段をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、各装置における処理手段がコンピュータ上で実現される。
【０１１１】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）/ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【０１１２】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記録装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【０１１３】
また、各装置の機能構成部は、コンピュータ上で所定のプログラムを実行させることにより構成することにしてもよいし、これらの処理内容の少なくとも一部をハードウェア的に実現することとしても良い。
【図１】
2011164126.tif 000036
【図２】
2011164126.tif 000037
【図３】
2011164126.tif 000038
【図４】
2011164126.tif 000039
【図５】
2011164126.tif 000040
【図６】
2011164126.tif 000041
【図７】
2011164126.tif 000042