ホーム > 特許ランキング > 学校法人東京理科大学
知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】5721098
(24)【登録日】2015年4月3日
(45)【発行日】2015年5月20日
(54)【発明の名称】雑音抑圧装置および雑音抑圧方法
(51)【国際特許分類】
G10L 21/0216 20130101AFI20150430BHJP
【ＦＩ】
G10L21/0216
【請求項の数】17
【全頁数】66
(21)【出願番号】特願2010-503784(P2010-503784)
(86)(22)【出願日】2009年3月18日
(86)【国際出願番号】JP2009001224
(87)【国際公開番号】WO2009116291
(87)【国際公開日】20090924
【審査請求日】2012年3月16日
(31)【優先権主張番号】特願2008-74691(P2008-74691)
(32)【優先日】2008年3月21日
(33)【優先権主張国】JP
(31)【優先権主張番号】特願2008-168835(P2008-168835)
(32)【優先日】2008年6月27日
(33)【優先権主張国】JP
(73)【特許権者】
【識別番号】803000115
【氏名又は名称】学校法人東京理科大学
(74)【代理人】
【識別番号】100105050
【弁理士】
【氏名又は名称】鷲田 公一
(72)【発明者】
【氏名】田邉 造
(72)【発明者】
【氏名】古川 利博
【審査官】
井上 健一
(56)【参考文献】
【文献】
特許第５１１５９５２（ＪＰ，Ｂ２）
【文献】
特開２００８−５８９０１（ＪＰ，Ａ）
【文献】
特開２００７−４１５０８（ＪＰ，Ａ）
【文献】
特開２００７−４１４９９（ＪＰ，Ａ）
【文献】
特開２００８−２９２８５８（ＪＰ，Ａ）
【文献】
田邉 造，有色性駆動雑音を考慮したカルマンフィルタを用いた雑音抑圧法，電子情報通信学会技術研究報告. EA, 応用音響 107(532)，２００８年 ２月２９日，ｐ．７９−８４
(58)【調査した分野】（Int.Cl.，ＤＢ名）
Ｇ１０Ｌ ２１／０２−２１／０２６４
(57)【特許請求の範囲】
【請求項1】
モデル化可能な観測対象のシステムにおいて、所望の情報を示す所望信号に不必要な情報である雑音が混在した観測信号のみから、前記所望信号を駆動源として含む状態空間モデルに基づく予測法を用いて、前記所望信号を推定する雑音抑圧装置であって、
前記観測信号を入力して各時刻の観測信号を生成する入力手段と、
時刻ｎまたは時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を算出する相関演算部と、
前記相関演算部によって算出された、時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数を算出する重み係数算出部と、
前記重み係数算出部によって算出された重み係数を用いて、時刻ｎ＋１での観測信号の推定誤差を補正し、補正した前記推定誤差と時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値とにより、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値を前記所望信号として算出する最適推定値算出部と、
を有し、
前記状態空間モデルは、任意の一時刻のみを抽出して構成された観測方程式を含み、前記相関演算部、前記重み係数算出部、および前記最適推定値算出部は、前記観測方程式に含まれる時刻ｎ＋１のみの観測信号に対し、反復的に、所定の処理を行い、
時刻ｎ＋１での観測信号の推定誤差は、時刻ｎ＋１での観測信号と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値から射影される観測信号の推定値との差であり、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値は、前記最適推定値算出部によって１つ前の時点での所望信号として算出された、時刻ｎまでの観測信号による時刻ｎでの前記状態量の最適推定値が、時刻ｎから１つの時刻だけ状態遷移したときの前記状態量の推定値である、
雑音抑圧装置。
前記観測信号を入力して各時刻の観測信号を生成する入力手段と、
時刻ｎまたは時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を算出する相関演算部と、
前記相関演算部によって算出された、時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数を算出する重み係数算出部と、
前記重み係数算出部によって算出された重み係数を用いて、時刻ｎ＋１での観測信号の推定誤差を補正し、補正した前記推定誤差と時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値とにより、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値を前記所望信号として算出する最適推定値算出部と、
を有し、
前記状態空間モデルは、任意の一時刻のみを抽出して構成された観測方程式を含み、前記相関演算部、前記重み係数算出部、および前記最適推定値算出部は、前記観測方程式に含まれる時刻ｎ＋１のみの観測信号に対し、反復的に、所定の処理を行い、
時刻ｎ＋１での観測信号の推定誤差は、時刻ｎ＋１での観測信号と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値から射影される観測信号の推定値との差であり、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値は、前記最適推定値算出部によって１つ前の時点での所望信号として算出された、時刻ｎまでの観測信号による時刻ｎでの前記状態量の最適推定値が、時刻ｎから１つの時刻だけ状態遷移したときの前記状態量の推定値である、
雑音抑圧装置。
【請求項2】
前記相関演算部は、
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値をスカラー量として算出し、
前記重み係数算出部は、
前記相関演算部によって算出された前記推定誤差の相関値のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をスカラー量として算出し、
前記最適推定値算出部は、
前記重み係数算出部によって算出された前記重み係数のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をスカラー量として算出する、
請求項１記載の雑音抑圧装置。
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値をスカラー量として算出し、
前記重み係数算出部は、
前記相関演算部によって算出された前記推定誤差の相関値のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をスカラー量として算出し、
前記最適推定値算出部は、
前記重み係数算出部によって算出された前記重み係数のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をスカラー量として算出する、
請求項１記載の雑音抑圧装置。
【請求項3】
前記相関演算部は、
与えられた駆動源の共分散のスカラー量を用いて、前記推定誤差の相関値のスカラー量の算出を行い、
前記重み係数算出部は、
与えられた雑音の共分散のスカラー量および前記相関演算部によって算出された前記推定誤差の相関値のスカラー量を用いて、前記重み係数のスカラー量の算出を行い、
前記最適推定値算出部は、
前記重み係数算出部によって算出された前記重み係数のスカラー量および時刻ｎ＋１のみの観測信号を用いて、前記状態量の最適推定値のスカラー量の算出を行う、
請求項２記載の雑音抑圧装置。
与えられた駆動源の共分散のスカラー量を用いて、前記推定誤差の相関値のスカラー量の算出を行い、
前記重み係数算出部は、
与えられた雑音の共分散のスカラー量および前記相関演算部によって算出された前記推定誤差の相関値のスカラー量を用いて、前記重み係数のスカラー量の算出を行い、
前記最適推定値算出部は、
前記重み係数算出部によって算出された前記重み係数のスカラー量および時刻ｎ＋１のみの観測信号を用いて、前記状態量の最適推定値のスカラー量の算出を行う、
請求項２記載の雑音抑圧装置。
【請求項4】
前記相関演算部は、第１の相関演算部および第２の相関演算部を有し、
前記最適推定値算出部は、第１の最適推定値算出部および第２の最適推定値算出部を有し、
前記第１の相関演算部は、
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を行列として算出し、
前記重み係数算出部は、
前記第１の相関演算部によって算出された前記推定誤差の相関値の行列を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をベクトル量として算出し、
前記第１の最適推定値算出部は、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の最適推定値算出部は、
前記重み係数算出部によって算出された前記重み係数のベクトル量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の相関演算部は、
時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記状態量を推定した場合の推定誤差の相関値を行列として算出する、
請求項１記載の雑音抑圧装置。
前記最適推定値算出部は、第１の最適推定値算出部および第２の最適推定値算出部を有し、
前記第１の相関演算部は、
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を行列として算出し、
前記重み係数算出部は、
前記第１の相関演算部によって算出された前記推定誤差の相関値の行列を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をベクトル量として算出し、
前記第１の最適推定値算出部は、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の最適推定値算出部は、
前記重み係数算出部によって算出された前記重み係数のベクトル量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の相関演算部は、
時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記状態量を推定した場合の推定誤差の相関値を行列として算出する、
請求項１記載の雑音抑圧装置。
【請求項5】
前記第１の相関演算部は、
所定の状態遷移行列、与えられた駆動源ベクトルの共分散の要素値、および与えられたまたは前回前記第２の相関演算部によって算出された前記推定誤差の相関値の行列を用いて、前記推定誤差の相関値の行列の算出を行い、
前記重み係数算出部は、
前記第１の相関演算部によって算出された前記推定誤差の相関値の行列、与えられた観測遷移ベクトル、および与えられた雑音の共分散のスカラー量を用いて、前記重み係数のベクトル量の算出を行い、
前記第１の最適推定値算出部は、
前記状態遷移行列、および、与えられたまたは前回前記第２の最適推定値算出部によって算出された前記状態量の最適推定値のベクトル量を用いて、前記状態量の最適推定値のベクトル量の算出を行い、
前記第２の最適推定値算出部は、
前記第１の最適推定値算出部によって算出された前記状態量の最適推定値のベクトル量、前記重み係数算出部によって算出された前記重み係数のベクトル量、前記観測遷移ベクトル、および時刻ｎ＋１のみの観測信号を用いて、前記状態量の最適推定値のベクトル量の算出を行い、
前記第２の相関演算部は、
前記重み係数算出部によって算出された前記重み係数のベクトル量、前記観測遷移ベクトル、および前記第１の相関演算部によって算出された前記推定誤差の相関値の行列を用いて、前記推定誤差の相関値の行列の算出を行う、
請求項４記載の雑音抑圧装置。
所定の状態遷移行列、与えられた駆動源ベクトルの共分散の要素値、および与えられたまたは前回前記第２の相関演算部によって算出された前記推定誤差の相関値の行列を用いて、前記推定誤差の相関値の行列の算出を行い、
前記重み係数算出部は、
前記第１の相関演算部によって算出された前記推定誤差の相関値の行列、与えられた観測遷移ベクトル、および与えられた雑音の共分散のスカラー量を用いて、前記重み係数のベクトル量の算出を行い、
前記第１の最適推定値算出部は、
前記状態遷移行列、および、与えられたまたは前回前記第２の最適推定値算出部によって算出された前記状態量の最適推定値のベクトル量を用いて、前記状態量の最適推定値のベクトル量の算出を行い、
前記第２の最適推定値算出部は、
前記第１の最適推定値算出部によって算出された前記状態量の最適推定値のベクトル量、前記重み係数算出部によって算出された前記重み係数のベクトル量、前記観測遷移ベクトル、および時刻ｎ＋１のみの観測信号を用いて、前記状態量の最適推定値のベクトル量の算出を行い、
前記第２の相関演算部は、
前記重み係数算出部によって算出された前記重み係数のベクトル量、前記観測遷移ベクトル、および前記第１の相関演算部によって算出された前記推定誤差の相関値の行列を用いて、前記推定誤差の相関値の行列の算出を行う、
請求項４記載の雑音抑圧装置。
【請求項6】
モデル化可能な観測対象のシステムにおいて、所望の情報を示す所望信号に不必要な情報である雑音が混在した観測信号のみから、前記所望信号を駆動源として含む状態空間モデルに基づく予測法を用いて、前記所望信号を推定する雑音抑圧方法であって、
前記観測信号を入力して各時刻の観測信号を生成する入力工程と、
時刻ｎまたは時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を算出する相関演算工程と、
前記相関演算工程で算出した、時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数を算出する重み係数算出工程と、
前記重み係数算出工程で算出した重み係数を用いて、時刻ｎ＋１での観測信号の推定誤差を補正し、補正した前記推定誤差と時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値とにより、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値を前記所望信号として算出する最適推定値算出工程と、
を有し、
前記状態空間モデルは、任意の一時刻のみを抽出して構成された観測方程式を含み、前記相関演算工程、前記重み係数算出工程、および前記最適推定値算出工程は、前記観測方程式に含まれる時刻ｎ＋１のみの観測信号に対し、反復的に、所定の処理を行い、
時刻ｎ＋１での観測信号の推定誤差は、時刻ｎ＋１での観測信号と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値から射影される観測信号の推定値との差であり、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値は、前記最適推定値算出工程で１つ前の時点での所望信号として算出した、時刻ｎまでの観測信号による時刻ｎでの前記状態量の最適推定値が、時刻ｎから１つの時刻だけ状態遷移したときの前記状態量の推定値である、
雑音抑圧方法。
前記観測信号を入力して各時刻の観測信号を生成する入力工程と、
時刻ｎまたは時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を算出する相関演算工程と、
前記相関演算工程で算出した、時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数を算出する重み係数算出工程と、
前記重み係数算出工程で算出した重み係数を用いて、時刻ｎ＋１での観測信号の推定誤差を補正し、補正した前記推定誤差と時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値とにより、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値を前記所望信号として算出する最適推定値算出工程と、
を有し、
前記状態空間モデルは、任意の一時刻のみを抽出して構成された観測方程式を含み、前記相関演算工程、前記重み係数算出工程、および前記最適推定値算出工程は、前記観測方程式に含まれる時刻ｎ＋１のみの観測信号に対し、反復的に、所定の処理を行い、
時刻ｎ＋１での観測信号の推定誤差は、時刻ｎ＋１での観測信号と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値から射影される観測信号の推定値との差であり、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値は、前記最適推定値算出工程で１つ前の時点での所望信号として算出した、時刻ｎまでの観測信号による時刻ｎでの前記状態量の最適推定値が、時刻ｎから１つの時刻だけ状態遷移したときの前記状態量の推定値である、
雑音抑圧方法。
【請求項7】
前記相関演算工程は、
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値をスカラー量として算出し、
前記重み係数算出工程は、
前記相関演算工程で算出した前記推定誤差の相関値のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をスカラー量として算出し、
前記最適推定値算出工程は、
前記重み係数算出工程で算出した前記重み係数のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をスカラー量として算出する、
請求項６記載の雑音抑圧方法。
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値をスカラー量として算出し、
前記重み係数算出工程は、
前記相関演算工程で算出した前記推定誤差の相関値のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をスカラー量として算出し、
前記最適推定値算出工程は、
前記重み係数算出工程で算出した前記重み係数のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をスカラー量として算出する、
請求項６記載の雑音抑圧方法。
【請求項8】
前記相関演算工程は、第１の相関演算工程および第２の相関演算工程を有し、
前記最適推定値算出工程は、第１の最適推定値算出工程および第２の最適推定値算出工程を有し、
前記第１の相関演算工程は、
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を行列として算出し、
前記重み係数算出工程は、
前記第１の相関演算工程で算出した前記推定誤差の相関値の行列を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をベクトル量として算出し、
前記第１の最適推定値算出工程は、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の最適推定値算出工程は、
前記重み係数算出工程で算出した前記重み係数のベクトル量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の相関演算工程は、
時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記状態量を推定した場合の推定誤差の相関値を行列として算出する、
請求項６記載の雑音抑圧方法。
前記最適推定値算出工程は、第１の最適推定値算出工程および第２の最適推定値算出工程を有し、
前記第１の相関演算工程は、
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を行列として算出し、
前記重み係数算出工程は、
前記第１の相関演算工程で算出した前記推定誤差の相関値の行列を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をベクトル量として算出し、
前記第１の最適推定値算出工程は、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の最適推定値算出工程は、
前記重み係数算出工程で算出した前記重み係数のベクトル量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の相関演算工程は、
時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記状態量を推定した場合の推定誤差の相関値を行列として算出する、
請求項６記載の雑音抑圧方法。
【請求項9】
モデル化可能な観測対象のシステムにおいて、所望の情報を示す所望信号に不必要な情報である雑音が混在した観測信号のみから、前記所望信号を駆動源として含む状態空間モデルに基づく予測法を用いて、前記所望信号を推定するための雑音抑圧プログラムであって、
コンピュータに、
各時刻の観測信号を入力する入力工程と、
時刻ｎまたは時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を算出する相関演算工程と、
前記相関演算工程で算出した、時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数を算出する重み係数算出工程と、
前記重み係数算出工程で算出した重み係数を用いて、時刻ｎ＋１での観測信号の推定誤差を補正し、補正した前記推定誤差と時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値とにより、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値を前記所望信号として算出する最適推定値算出工程と、
を実行させ、
前記状態空間モデルは、任意の一時刻のみを抽出して構成された観測方程式を含み、前記相関演算工程、前記重み係数算出工程、および前記最適推定値算出工程は、前記観測方程式に含まれる時刻ｎ＋１のみの観測信号に対し、反復的に、所定の処理を行い、
時刻ｎ＋１での観測信号の推定誤差は、時刻ｎ＋１での観測信号と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値から射影される観測信号の推定値との差であり、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値は、前記最適推定値算出工程で１つ前の時点での所望信号として算出した、時刻ｎまでの観測信号による時刻ｎでの前記状態量の最適推定値が、時刻ｎから１つの時刻だけ状態遷移したときの前記状態量の推定値である、
雑音抑圧プログラム。
コンピュータに、
各時刻の観測信号を入力する入力工程と、
時刻ｎまたは時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を算出する相関演算工程と、
前記相関演算工程で算出した、時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数を算出する重み係数算出工程と、
前記重み係数算出工程で算出した重み係数を用いて、時刻ｎ＋１での観測信号の推定誤差を補正し、補正した前記推定誤差と時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値とにより、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値を前記所望信号として算出する最適推定値算出工程と、
を実行させ、
前記状態空間モデルは、任意の一時刻のみを抽出して構成された観測方程式を含み、前記相関演算工程、前記重み係数算出工程、および前記最適推定値算出工程は、前記観測方程式に含まれる時刻ｎ＋１のみの観測信号に対し、反復的に、所定の処理を行い、
時刻ｎ＋１での観測信号の推定誤差は、時刻ｎ＋１での観測信号と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値から射影される観測信号の推定値との差であり、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値は、前記最適推定値算出工程で１つ前の時点での所望信号として算出した、時刻ｎまでの観測信号による時刻ｎでの前記状態量の最適推定値が、時刻ｎから１つの時刻だけ状態遷移したときの前記状態量の推定値である、
雑音抑圧プログラム。
【請求項10】
前記相関演算工程は、
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値をスカラー量として算出し、
前記重み係数算出工程は、
前記相関演算工程で算出した前記推定誤差の相関値のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をスカラー量として算出し、
前記最適推定値算出工程は、
前記重み係数算出工程で算出した前記重み係数のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をスカラー量として算出する、
請求項９記載の雑音抑圧プログラム。
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値をスカラー量として算出し、
前記重み係数算出工程は、
前記相関演算工程で算出した前記推定誤差の相関値のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をスカラー量として算出し、
前記最適推定値算出工程は、
前記重み係数算出工程で算出した前記重み係数のスカラー量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をスカラー量として算出する、
請求項９記載の雑音抑圧プログラム。
【請求項11】
前記相関演算工程は、第１の相関演算工程および第２の相関演算工程を有し、
前記最適推定値算出工程は、第１の最適推定値算出工程および第２の最適推定値算出工程を有し、
前記第１の相関演算工程は、
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を行列として算出し、
前記重み係数算出工程は、
前記第１の相関演算工程で算出した前記推定誤差の相関値の行列を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をベクトル量として算出し、
前記第１の最適推定値算出工程は、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の最適推定値算出工程は、
前記重み係数算出工程で算出した前記重み係数のベクトル量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の相関演算工程は、
時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記状態量を推定した場合の推定誤差の相関値を行列として算出する、
請求項９記載の雑音抑圧プログラム。
前記最適推定値算出工程は、第１の最適推定値算出工程および第２の最適推定値算出工程を有し、
前記第１の相関演算工程は、
時刻ｎまでの観測信号により時刻ｎ＋１での前記システムの状態量を推定した場合の推定誤差の相関値を行列として算出し、
前記重み係数算出工程は、
前記第１の相関演算工程で算出した前記推定誤差の相関値の行列を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値と、時刻ｎ＋１での観測信号の推定誤差と、の関係を規定するための重み係数をベクトル量として算出し、
前記第１の最適推定値算出工程は、
時刻ｎまでの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の最適推定値算出工程は、
前記重み係数算出工程で算出した前記重み係数のベクトル量を用いて、時刻ｎ＋１までの観測信号による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出し、
前記第２の相関演算工程は、
時刻ｎ＋１までの観測信号により時刻ｎ＋１での前記状態量を推定した場合の推定誤差の相関値を行列として算出する、
請求項９記載の雑音抑圧プログラム。
【請求項12】
胎児の心拍を検出して得られた信号を観測信号として、所望信号としての胎児の心拍を推定する、請求項１から請求項５のいずれかに記載の雑音抑圧装置と、
前記雑音抑圧装置によって推定された胎児の心拍を出力する出力手段と、
を有する胎児心音検出装置。
前記雑音抑圧装置によって推定された胎児の心拍を出力する出力手段と、
を有する胎児心音検出装置。
【請求項13】
使用者の音声を電気的な音声信号として検出するマイクロホンと、
無線送信された通信相手の音声信号を受信する受信手段と、
前記マイクロホンによって検出された音声信号または前記受信手段によって受信された音声信号を観測信号として、所望信号としての使用者または通信相手の音声信号を推定する、請求項１から請求項５のいずれかに記載の雑音抑圧装置と、
前記雑音抑圧装置によって推定された使用者の音声信号を無線送信する送信手段と、
前記雑音抑圧装置によって推定された通信相手の音声信号を出力するスピーカと、
を有する携帯端末装置。
無線送信された通信相手の音声信号を受信する受信手段と、
前記マイクロホンによって検出された音声信号または前記受信手段によって受信された音声信号を観測信号として、所望信号としての使用者または通信相手の音声信号を推定する、請求項１から請求項５のいずれかに記載の雑音抑圧装置と、
前記雑音抑圧装置によって推定された使用者の音声信号を無線送信する送信手段と、
前記雑音抑圧装置によって推定された通信相手の音声信号を出力するスピーカと、
を有する携帯端末装置。
【請求項14】
使用者の音声を電気的な音声信号として検出するマイクロホンと、
前記マイクロホンによって検出された音声信号を観測信号として、所望信号としての使用者の音声信号を推定する、請求項１から請求項５のいずれかに記載の雑音抑圧装置と、
前記雑音抑圧装置によって推定された使用者の音声信号を入力して音声認識処理を行う音声認識手段と、
前記音声認識手段の認識結果を操作指示として、カーナビゲーション処理を実行するカーナビゲーション処理手段と、
を有するカーナビゲーション装置。
前記マイクロホンによって検出された音声信号を観測信号として、所望信号としての使用者の音声信号を推定する、請求項１から請求項５のいずれかに記載の雑音抑圧装置と、
前記雑音抑圧装置によって推定された使用者の音声信号を入力して音声認識処理を行う音声認識手段と、
前記音声認識手段の認識結果を操作指示として、カーナビゲーション処理を実行するカーナビゲーション処理手段と、
を有するカーナビゲーション装置。
【請求項15】
使用者の音声を検出して得られた音声信号を観測信号として、所望信号としての使用者の音声信号を推定する、請求項１から請求項５のいずれかに記載の雑音抑圧装置と、
前記雑音抑圧装置によって推定された使用者の音声信号を入力して音声認識処理を行う音声認識手段と、
前記音声認識手段の認識結果を出力する出力手段と、
を有する音声認識装置。
前記雑音抑圧装置によって推定された使用者の音声信号を入力して音声認識処理を行う音声認識手段と、
前記音声認識手段の認識結果を出力する出力手段と、
を有する音声認識装置。
【請求項16】
検査対象が発する音を検出して得られた信号を観測信号として、所望信号を検査対象の信号として推定する、請求項１から請求項５のいずれかに記載の雑音抑圧装置と、
前記雑音抑圧装置によって推定された所望信号を正常音と比較して異常音を検出する異常音解析手段と、
前記異常音解析手段の解析結果を出力する出力手段と、
を有する異常検出装置。
前記雑音抑圧装置によって推定された所望信号を正常音と比較して異常音を検出する異常音解析手段と、
前記異常音解析手段の解析結果を出力する出力手段と、
を有する異常検出装置。
【請求項17】
画像信号を入力する画像入力手段と、
前記画像入力手段によって入力された画像信号を観測信号として、所望信号としての原画像を推定する、請求項１から請求項５のいずれかに記載の雑音抑圧装置と、
前記雑音抑圧装置によって推定された原画像を出力する出力手段と、
を有する画像処理装置。
前記画像入力手段によって入力された画像信号を観測信号として、所望信号としての原画像を推定する、請求項１から請求項５のいずれかに記載の雑音抑圧装置と、
前記雑音抑圧装置によって推定された原画像を出力する出力手段と、
を有する画像処理装置。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、雑音抑圧装置および雑音抑圧方法に関する。【背景技術】
【0002】
所望の情報（所望信号）に不必要な情報（雑音）が混在した観測情報（雑音などにより破損した情報）から不必要な情報を取り除き、所望情報のみを抽出することは、音声や無線通信、画像、姿勢制御、認識、産業用・福祉用・医療用ロボットなどの分野における重要な技術であり、近年盛んに研究開発が行われている。【0003】
例えば、音声分野における公知の雑音抑圧方法としては、単一のマイクロホンを用いた方法や、複数のマイクロホンから構成されるマイクロホンアレイを用いた方法が提案されている。【0004】
しかしながら、マイクロホンアレイを用いた方法では、雑音の音源数以上のマイクロホンが必要であるため、雑音の音源数が増大すると、マイクロホンの数も比例して増加することが避けられず、コストが増大する。また、携帯電話などの通信の小型化傾向における商品において、マイクロホンを装着できる数に制限がある場合や、各マイクロホンの特性の違いを制御する場合など、実用化が困難な場合も存在する。そのため、現在は、単一のマイクロホンを用いた雑音抑圧方法の開発が主流となっている。【0005】
単一のマイクロホンしか用いない従来の雑音抑圧方法のアルゴリズムとしては、以下のようなものが知られている。【0006】
非特許文献１記載のＡＮＣ（適応ノイズキャンセラ）アルゴリズムは、音声信号の周期性を利用してノイズ信号を低減する。【0007】
非特許文献２には、線形予測に基づく雑音抑圧アルゴリズムが記載されている。このアルゴリズムは、非特許文献１記載のＡＮＣで必要であったピッチ推定や、雑音パワースペクトラム、雑音の平均方向に関する事前知識を必要としない。【0008】
また、上記アルゴリズムとは別に、カルマンフィルタに基づく雑音抑圧アルゴリズムが、非特許文献３に提案されている。このアルゴリズムは、観測信号から音声信号の自己回帰（ＡＲ：Autoregressive）システムをモデル化する。さらに、このアルゴリズムは、ＡＲシステムのパラメータ（以下「ＡＲ係数」という）を推定し、推定したＡＲ係数を用いてカルマンフィルタに基づき雑音抑圧を実行する。【0009】
カルマンフィルタに基づくアルゴリズムの多くは、通常、２段階で動作する。すなわち、このようなアルゴリズムは、最初にＡＲ係数を推定し、次に推定したＡＲ係数を用いてカルマンフィルタに基づき雑音抑圧を行う。【非特許文献1】J.R. Deller, J.G. Proakis, J.H.L. Hansen, "Discrete-Time Processing of Speech Signals," Macmillan Press, 1993
【非特許文献2】A. Kawamura, K. Fujii, Y. Itoh and Y. Fukui, “A Noise Reduction Method Based on Linear Prediction Analysis,” IEICE Trans. Fundamentals, vol.J85-A, no.4, pp.415-423, May 2002
【非特許文献3】W. Kim and H. Ko, "Noise Variance Estimation for Kalman Filtering of Noise Speech," IEICE Trans. Inf. & syst., vol.E84-D, no.1, pp.155-160, Jan 2001
【非特許文献4】N. Tanabe, T. Inoue, K. Sueyoshi, T. Furukawa, H. Kubota, H. Matsue, and S. Tsujii, "Robust Noise Suppression Algorithm using Kalman Filter Theory with Colored Driving Source," IEICE Technical Report, EA2007-125, pp.79-84, Mar. 2008
【発明の開示】
【発明が解決しようとする課題】
【0010】
しかしながら、非特許文献１に記載された公知のアルゴリズムは、音声信号のピッチ周期の正確な推定を必要とする。そのため、このアルゴリズムは、その雑音抑圧能力が雑音によって劣化してしまうという問題点を有している。【0011】
この点、非特許文献２記載のアルゴリズムは、音声信号のピッチ周期の正確な推定を必要とせずに、雑音抑圧が可能である。さらに、このアルゴリズムは、その原理が単純であり、演算量を少なくすることができるといった長所を有している。しかし、このアルゴリズムは、その雑音抑圧能力が入力音声信号の周期性や線形性などの特性に依存している。言い換えると、このアルゴリズムは、アルゴリズムの中に音声信号に依存するパラメータが存在しているため、その実用には一定の限界がある。【0012】
非特許文献３記載のアルゴリズムは、高い雑音抑圧能力を有し、特に高い音質を得たい音響分野への応用に適した手法である。【0013】
しかしながら、一方で、このアルゴリズムは、ＡＲ係数を必要とするため、ＡＲ係数の推定精度に雑音抑圧能力（つまり、当該カルマンフィルタアルゴリズムの性能）が大きく依存してしまうという問題点を有している。すなわち、ＡＲ係数が正確に推定されない場合、雑音を抑圧し切れないのみならず、場合によっては雑音に加えて音声信号まで抑圧してしまう可能性がある。これらは、雑音が抑圧された音声信号の音質の劣化を引き起こす要因となりうる。【0014】
この点、一般には、ＡＲ係数の正確な推定は困難である。ＡＲ係数の正確な推定は、例えば、雑音抑圧であれば、クリアな信号つまり所望信号（例えば、音声信号）に依存しているからである。このことは、音声信号が既知でなければならないことを意味しているため、リアルタイム処理は困難となる。また、仮に何らかの手法でリアルタイムにＡＲ係数を正確に推定することが可能となったとしても、処理が増加するため演算量の問題は避けられない。また、そもそも、ＡＲ係数の推定はＡＲ係数の次数を決定した後に行われるが、ＡＲ係数の次数を決定するのは非常に困難であり、この点からもＡＲ係数の正確な推定は困難であるといえる。【0015】
そこで、本発明者は、カルマンフィルタに基づく従来の雑音抑圧方法（非特許文献３参照）の問題点を解決するために、非特許文献４に記載の雑音抑圧方法を提案した。具体的には、カルマンフィルタに基づく従来の雑音抑圧方法では、線形予測を用いてＡＲ係数を推定した後、その結果を用いてカルマンフィルタを実行することによって雑音抑圧を実現しているのに対し、この提案手法では、状態方程式と観測方程式から構成される新しい予測法によって雑音抑圧を実現している。そのため、この提案手法では、新しい状態空間モデル（状態方程式と観測方程式）を構成し、具体的には、情報源からのクリアな信号つまり所望信号（例えば、音声信号）のみを用いて状態方程式を構成し、そのクリアな信号と雑音とを用いて観測方程式を構成している。【0016】
しかしながら、この提案手法の状態空間モデルでは、観測信号をベクトル化することにより、多くの過去情報を用いて雑音抑圧を実行している。しかし、過去の観測信号にも雑音が混入するため（つまり、誤差あり）、多くの過去情報を用いた観測信号ベクトルは雑音（誤差）を含んだものとなる。すなわち、予測の推定精度を上げるために過去の観測信号が多く必要であるという命題は、必ずしも正しくない。したがって、この提案手法においては、所望信号の推定精度（つまり、雑音抑圧能力）の向上に一定の限界がある。また、この提案手法のアルゴリズムでは、逆行列の計算が必要であるため、カルマンフィルタに基づく従来の雑音抑圧方法に対する演算量の軽減効果にも一定の限界がある。なお、多くの過去情報を用いる点は、他の従来の雑音抑圧方法（非特許文献１〜３参照）も同様である。【0017】
本発明の目的は、よりシンプルな構成で、かつ、より少ない演算量で、所望情報の品質を劣化させることなくより高い雑音抑圧能力を実現することができる雑音抑圧装置および雑音抑圧方法を提供することである。【課題を解決するための手段】
【0018】
本発明の雑音抑圧装置は、所望の情報に雑音が混在した観測情報のみから前記所望情報を推定する雑音抑圧装置であって、時刻ｎのみの観測情報に対して、時刻ｎまたは時刻ｎ＋１までの情報により前記所望情報を含む時刻ｎ＋１のシステムの状態量を推定した場合の推定誤差の相関値を算出する相関演算部と、時刻ｎのみの観測情報に対して、前記相関演算部によって算出された相関値を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値と、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値と、前記観測情報を含む観測量の推定誤差と、の関係を規定するための重み係数を算出する重み係数算出部と、時刻ｎのみの観測情報に対して、前記重み係数算出部によって算出された重み係数を用いて、時刻ｎまたは時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値を算出する最適推定値算出部と、を有する構成を採る。【0019】
本発明の雑音抑圧装置は、好ましくは、所望の情報に雑音が混在した観測情報のみから前記所望情報を推定する雑音抑圧装置であって、時刻ｎのみの観測情報に対して、時刻ｎまでの情報により前記所望情報を含む時刻ｎ＋１のシステムの状態量を推定した場合の推定誤差の相関値をスカラー量として算出する相関演算部と、時刻ｎのみの観測情報に対して、前記相関演算部によって算出された前記推定誤差の相関値のスカラー量を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値と、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値と、前記観測情報を含む観測量の推定誤差と、の関係を規定するための重み係数をスカラー量として算出する重み係数算出部と、時刻ｎのみの観測情報に対して、前記重み係数算出部によって算出された前記重み係数のスカラー量を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値をスカラー量として算出する最適推定値算出部と、を有する構成を採る。【0020】
本発明の雑音抑圧装置は、好ましくは、所望の情報に雑音が混在した観測情報のみから前記所望情報を推定する雑音抑圧装置であって、時刻ｎのみの観測情報に対して、時刻ｎまでの情報により前記所望情報を含む時刻ｎ＋１のシステムの状態量を推定した場合の推定誤差の相関値を行列として算出する第１の相関演算部と、時刻ｎのみの観測情報に対して、前記第１の相関演算部によって算出された前記推定誤差の相関値の行列を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値と、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値と、前記観測情報を含む観測量の推定誤差と、の関係を規定するための重み係数をベクトル量として算出する重み係数算出部と、時刻ｎのみの観測情報に対して、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出する第１の最適推定値算出部と、時刻ｎのみの観測情報に対して、前記重み係数算出部によって算出された前記重み係数のベクトル量を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値をベクトル量として算出する第２の最適推定値算出部と、時刻ｎのみの観測情報に対して、時刻ｎ＋１までの情報により当該時刻の前記状態量を推定した場合の推定誤差の相関値を行列として算出する第２の相関演算部と、を有する構成を採る。【0021】
本発明の雑音抑圧方法は、所望の情報に雑音が混在した観測情報のみから前記所望情報を推定する雑音抑圧方法であって、時刻ｎのみの観測情報に対して、時刻ｎまたは時刻ｎ＋１までの情報により前記所望情報を含む時刻ｎ＋１のシステムの状態量を推定した場合の推定誤差の相関値を算出する相関演算工程と、時刻ｎのみの観測情報に対して、前記相関演算工程で算出した相関値を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値と、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値と、前記観測情報を含む観測量の推定誤差と、の関係を規定するための重み係数を算出する重み係数算出工程と、時刻ｎのみの観測情報に対して、前記重み係数算出工程で算出した重み係数を用いて、時刻ｎまたは時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値を算出する最適推定値算出工程と、を有するようにした。【0022】
本発明の雑音抑圧方法は、好ましくは、所望の情報に雑音が混在した観測情報のみから前記所望情報を推定する雑音抑圧方法であって、時刻ｎのみの観測情報に対して、時刻ｎまでの情報により前記所望情報を含む時刻ｎ＋１のシステムの状態量を推定した場合の推定誤差の相関値をスカラー量として算出する相関演算工程と、時刻ｎのみの観測情報に対して、前記相関演算工程で算出した前記推定誤差の相関値のスカラー量を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値と、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値と、前記観測情報を含む観測量の推定誤差と、の関係を規定するための重み係数をスカラー量として算出する重み係数算出工程と、時刻ｎのみの観測情報に対して、前記重み係数算出工程で算出した前記重み係数のスカラー量を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値をスカラー量として算出する最適推定値算出工程と、を有するようにした。【0023】
本発明の雑音抑圧方法は、好ましくは、所望の情報に雑音が混在した観測情報のみから前記所望情報を推定する雑音抑圧方法であって、時刻ｎのみの観測情報に対して、時刻ｎまでの情報により前記所望情報を含む時刻ｎ＋１のシステムの状態量を推定した場合の推定誤差の相関値を行列として算出する第１の相関演算工程と、時刻ｎのみの観測情報に対して、前記第１の相関演算工程で算出した前記推定誤差の相関値の行列を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値と、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値と、前記観測情報を含む観測量の推定誤差と、の関係を規定するための重み係数をベクトル量として算出する重み係数算出工程と、時刻ｎのみの観測情報に対して、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値をベクトル量として算出する第１の最適推定値算出工程と、時刻ｎのみの観測情報に対して、前記重み係数算出工程で算出した前記重み係数のベクトル量を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値をベクトル量として算出する第２の最適推定値算出工程と、時刻ｎのみの観測情報に対して、時刻ｎ＋１までの情報により当該時刻までの前記状態量を推定した場合の推定誤差の相関値を行列として算出する第２の相関演算工程と、を有するようにした。【発明の効果】
【0024】
本発明によれば、よりシンプルな構成で、かつ、より少ない演算量で、所望情報の品質を劣化させることなくより高い雑音抑圧能力を実現することができる。【図面の簡単な説明】
【0025】
【図1】本発明の実施の形態１に係る雑音抑圧装置の構成を示すブロック図
【図3】本発明が前提とするモデルを示す図
【図4】従来手法１の状態空間モデルのシステム構成を表すブロック線図
【図5】従来手法１のアルゴリズムを示す図
【図8】従来手法１とその問題点を視覚的にまとめた説明図
【図9】従来手法３の状態空間モデルのシステム構成を表すブロック線図
【図10】従来手法３のアルゴリズムを示す図
【図13】従来手法３を視覚的にまとめた説明図
【図14】従来手法３の問題点を視覚的にまとめた説明図
【図15】従来手法３の他の問題点を視覚的にまとめた説明図
【図16】本実施の形態の発明手法１の状態空間モデルを視覚的に説明するための概略図
【図17】発明手法１の状態空間モデルのシステム構成を表すブロック線図
【図18】発明手法１の状態方程式を視覚的に示す説明図
【図19】発明手法１における新しい観測方程式の構成を視覚的に示す説明図
【図20】発明手法１における新しい観測方程式の構成を視覚的に示す他の説明図
【図21】発明手法１のアルゴリズムを示す図
【図22】発明手法１のアルゴリズムの利点を視覚的にまとめた説明図
【図24】本実施の形態の発明手法２の導出方法を視覚的に示す段階別説明図
【図27】発明手法２のアルゴリズムを示す図
【図28】発明手法２のアルゴリズムの利点を視覚的にまとめた説明図
【図30】（Ａ）は、従来手法１の演算量を説明するための図、（Ｂ）は、従来手法３の演算量を説明するための図、（Ｃ）は、発明手法１の演算量を説明するための図、（Ｄ）は、発明手法２の演算量を説明するための図
【図31】主に従来手法１と従来手法３および発明手法１〜２との比較を示すグラフ図
【図32】従来手法３と発明手法１〜２との比較を示すグラフ図
【図33A-D】各手法の音声波形シミュレーションの第１の例を説明するための図
【図33E-H】各手法の音声波形シミュレーションの第１の例を説明するための図
【図34A-D】各手法の音声波形シミュレーションの第２の例を説明するための図
【図34E-H】各手法の音声波形シミュレーションの第２の例を説明するための図
【図35A-D】各手法の音声波形シミュレーションの第３の例を説明するための図
【図35E-H】各手法の音声波形シミュレーションの第３の例を説明するための図
【図36A-D】各手法の音声波形シミュレーションの第４の例を説明するための図
【図36E-H】各手法の音声波形シミュレーションの第４の例を説明するための図
【図37】各手法の雑音抑圧能力の数値シミュレーションの第１の例の結果を示す図
【図38】各手法の雑音抑圧能力の数値シミュレーションの第２の例の結果を示す図
【図39】各手法の雑音抑圧能力の数値シミュレーションの第３の例の結果を示す図
【図40】各手法の雑音抑圧能力の数値シミュレーションの第４の例の結果を示す図
【図41】各手法の主観的評価の第１の例の結果を示す図
【図42】各手法の主観的評価の第２の例の結果を示す図
【図43】各手法の主観的評価の第３の例の結果を示す図
【図44】各手法の主観的評価の第４の例の結果を示す図
【図45】本発明の実施の形態２に係る胎児心拍検出装置の構成を示すブロック図
【図46】本発明の実施の形態３に係る携帯端末装置の構成を示すブロック図
【図47】本発明の実施の形態４に係るカーナビゲーション装置の構成を示すブロック図
【図48】本発明の実施の形態５に係る音声認識装置の構成を示すブロック図
【図49】本発明の実施の形態６に係る異常検出装置の構成を示すブロック図
【図50】本発明の実施の形態７に係る画像復元装置の構成を示すブロック図
【発明を実施するための最良の形態】
【0026】
以下、本発明の実施の形態について、図面を参照して詳細に説明する。【0027】
なお、以下の説明においては、便宜上、非特許文献３に記載された、カルマンフィルタに基づく従来の雑音抑圧方法を「従来手法１」、カルマンフィルタ以外の、非特許文献２に記載された、線形予測に基づく従来の雑音抑圧方法を「従来手法２」、本発明者の提案に係る、非特許文献４に記載された雑音抑圧方法を「従来手法３」、とそれぞれ呼ぶことにする。【0028】
また、本明細書において、「雑音」とは、通常は、所望の情報以外の不必要な情報をすべて含む広義の意味で使用するが、特に画像の場合には、画像劣化の要因の１つである「ぼけ」と区別するため、「ぼけ」を除く狭義の意味で使用することにする。また、「音声」とは、人の声に限定されず、人の声を含む音一般を広く意味する。【0029】
（実施の形態１）図１は、本発明の実施の形態１に係る雑音抑圧装置の構成を示すブロック図である。
【0030】
図１に示す雑音抑圧装置１００は、入力部１１０、サンプリング部１２０、Ａ／Ｄ変換部１３０、バッファ１４０、雑音抑圧処理部１５０、および出力部１６０を有する。【0031】
入力部１１０は、観測情報または観測信号を入力する。観測信号は、情報源からのクリアな信号（所望信号）と、雑音とが合わさった（混合した）信号である。入力部１１０は、例えば、入力したアナログの観測信号を入力処理して、サンプリング部１２０に出力する。入力処理は、例えば、帯域制限処理や自動利得制御処理などである。【0032】
サンプリング部１２０は、所定のサンプリング周波数（例えば、１６ｋＨｚ）で、入力されたアナログの観測信号をサンプリング処理し、Ａ／Ｄ変換部１３０に出力する。サンプリング周波数は、検出対象（情報源）に応じて変更可能である。【0033】
Ａ／Ｄ変換部１３０は、サンプリングされた観測信号の振幅値を所定の分解能（例えば、８ｂｉｔ）でＡ／Ｄ変換処理し、バッファ１４０に送る。バッファ１４０は、所定のサンプリング数Ｎの信号フレーム（ブロック）を雑音抑圧処理部１５０に出力する。【0034】
雑音抑圧処理部１５０は、本発明の特徴的な構成要素であり、後述する雑音抑圧アルゴリズムを内蔵している。例えば、カルマンフィルタに基づく従来手法１では、線形予測を用いてＡＲ係数を推定した後、その結果を用いてカルマンフィルタを実行することにより雑音抑圧を実現しているのに対し、本発明の雑音抑圧方法（以下「発明手法」という）では、大きな基本的枠組に関しては、従来手法３と同様に、状態方程式と観測方程式から構成される新しい予測法によって雑音抑圧を実現している。すなわち、発明手法では、ＡＲ係数の推定を必要としない新しい状態空間モデル（状態方程式と観測方程式）を構成し、具体的には、情報源からのクリアな信号（所望信号）のみを用いて状態方程式を構成し、そのクリアな信号と雑音とを用いて観測方程式を構成している。しかし、発明手法では、その基本的枠組を実現するための具体的な構成に関しては、観測情報または観測信号について、従来手法３のように多くの過去情報を用いて雑音抑圧を実行するのではなく、ある１つの時刻の情報（瞬時情報）のみを用いて高性能な雑音抑圧を実行している。雑音抑圧処理部１５０は、内蔵する雑音抑圧アルゴリズムを用いて、ある１つの時刻の観測信号のみから所望信号（情報源からのクリアな信号）を推定する。雑音抑圧処理部１５０によって推定された推定信号は、雑音抑圧処理部１５０内に一時保存された後、出力部１６０に出力される。【0035】
このように、本明細書において「雑音抑圧」とは、観測信号のみから所望信号（例えば、音声信号）を推定することをいい、雑音を推定して引き算により所望信号を推定する「雑音除去」とは、引き算を行わない点で、全く異なる概念である。また、本明細書において「雑音」とは、通常、上記のように、所望の情報以外の不必要な情報すべて、換言すれば、観測信号に含まれるすべての信号のうち、所望信号以外のすべての信号を意味する。【0036】
出力部１６０は、雑音抑圧処理部１５０から入力した推定信号を所定の形態で外部に出力する。出力部１６０は、例えば、スピーカやディスプレイ、通信手段、記憶装置などで構成されている。【0037】
図２は、図１の雑音抑圧処理部１５０の構成を示すブロック図である。【0038】
図２に示すように、雑音抑圧処理部１５０は、初期設定部１５２、相関演算部１５４、重み係数算出部１５６、および最適推定値算出部１５８を有する。雑音抑圧処理部１５０は、各部１５２〜１５８の協働の下、内蔵する雑音抑圧アルゴリズム（発明手法）を実行することによって、ある１つの時刻の観測信号のみから所望信号（情報源からのクリアな信号）を推定する。このとき、初期設定部１５２は、発明手法のアルゴリズムの初期設定を行い、相関演算部１５４は、所望信号の推定誤差の相関演算を行い、重み係数算出部１５６は、所望信号の最適推定値の算出に必要な重み係数の算出を行い、最適推定値算出部１５８は、所望信号の最適推定値の算出を行う。なお、各部１５２〜１５８の具体的な処理内容は、後で詳細に説明する。【0039】
以下、雑音抑圧処理部１５０で行われる雑音抑圧処理動作について詳細に説明するが、ここでは、発明手法の特徴をより明確にするため、まず従来手法（具体的には、従来手法１および従来手法３）を説明し、この従来手法との対比において発明手法を詳細に説明することにする。なお、ここでは、情報源からのクリアな信号（所望信号）として、例えば、音声信号を例にとって説明する。【0040】
図３は、本発明が前提とするモデルを示す図である。雑音抑圧処理部１５０に入力される観測信号ｒ（ｎ）は、情報源からのクリアな信号（所望信号）（例えば、音声信号）ｄ（ｎ）以外に雑音ｖ（ｎ）を含んでおり、次の式（１）を満たす。【0041】
【数1】
【0042】
すなわち、音声信号を例にとると、一般に、人間の耳に聞こえる観測信号ｒ（ｎ）は、音声信号ｄ（ｎ）に雑音ｖ（ｎ）が加わった図３のモデルで表現されるとともに、数式では上記の式（１）のように書き表される。ただし、雑音ｖ（ｎ）は、分散σ２ｖとする。また、「ｎ」とは、装置の時刻ｎである。時刻ｎは、サンプリング部１２０で生成された離散的な時間系列において、処理開始時刻を時刻０と仮定したときに、そこからｎ番目の時刻のことを意味する。上記のように、本発明の目的は、ある１つの時刻の観測信号ｒ（ｎ）のみからその時刻の高音質な音声信号ｄ（ｎ）を復元（推定）することである。【0043】
［従来手法１］まず、従来手法１について説明する。
【0044】
従来のカルマンフィルタを用いた雑音抑圧方法は、まず第１段階（ステップ１）で、音声信号ｄ（ｎ）のＡＲ係数を推定した後、第２段階（ステップ２）で、ステップ１で推定したＡＲ係数を用いてカルマンフィルタアルゴリズムを実行することにより、観測信号ｒ（ｎ）から音声信号ｄ（ｎ）を抽出する。すなわち、従来手法１は、ステップ１で、線形予測により音声信号に対するＡＲシステムのパラメータ（ＡＲ係数）を推定し（ＡＲ係数の推定）、ステップ２で、ステップ１で推定したＡＲ係数を用いてカルマンフィルタにより雑音抑圧を行う（音声信号の復元）。【0045】
ステップ１（ＡＲ係数の推定）において、音声信号ｄ（ｎ）がＡＲ過程で表されると仮定すると、式（１）で与えられる音声信号ｄ（ｎ）は、次の式（２）のように表される。【0046】
【数2】
【0047】
ここで、αｌ（ｎ）は、時刻ｎでのＡＲ係数、Ｌｃ１は、ＡＲ係数の次数、ｅ（ｎ）は、駆動源、つまり、音声信号ｄ（ｎ）が式（２）に示すＬｃ１次のＡＲシステムでモデル化されるとした場合の予測誤差（モデリング誤差）である。また、駆動源は、ゼロ平均の白色ガウス過程とする。なお、添え字「ｃ１」は、従来技術１に係るものであることを示している。【0048】
公知のように、従来手法１では、雑音ｖ（ｎ）は、ゼロ平均であり白色雑音であることが前提条件である。言い換えると、従来手法１では、音声信号ｄ（ｎ）と雑音ｖ（ｎ）は無相関であると仮定されている。【0049】
すなわち、ステップ２（音声信号の復元）において、カルマンフィルタ理論に基づいて状態空間モデル（状態方程式と観測方程式）を構成するためには、駆動源δｃ１（ｎ）が白色信号であり、かつ、音声信号からなる状態量ｘｃ１（ｎ）と雑音ｖ（ｎ）とが無相関でなければならない。このことを踏まえた上で、ステップ１で求めたＡＲ係数αｌ（ｎ）を用いた従来手法の状態空間モデル（状態方程式と観測方程式）は、次の式（３）で書き表される。ただし、式（３）中のベクトルｘｃ１、δｃ１、ｍｃ１、スカラーｙｃ１、εｃ１、および行列Φｃ１は、次の式（４）でそれぞれ定義される。ベクトルｘｃ１は音声信号からなるＬｃ１×１次の状態ベクトル、ベクトルδｃ１はＬｃ１×１次の駆動源ベクトル、スカラーｙｃ１は観測信号、スカラーεｃ１は雑音、行列Φｃ１はＬｃ１×Ｌｃ１次の状態遷移行列、ベクトルｍｃ１はＬｃ１×１次の観測遷移ベクトルである。また、図４は、この状態空間モデルにおけるシステム構成図をブロック線図で表したものである。【0050】
【数3】
【0051】
【数4】
【0052】
図４において、「４０」は時刻ｎにおける状態ベクトルｘｃ１（ｎ）、「４１」は時刻ｎ＋１における状態ベクトルｘｃ１（ｎ＋１）、「４２」は時刻ｎにおける観測信号ｙｃ１（ｎ）、「４３」は時刻ｎにおける雑音εｃ１（ｎ）、「４４」は時刻ｎ＋１における駆動源ベクトルδｃ１（ｎ＋１）、「４５」は状態遷移行列Φｃ１、「４６」は観測遷移ベクトルｍｃ１である。式（３）における状態方程式は、推定対象（例えば、音声）のシステムを状態空間モデルで記述したものであり、内部状態つまり状態変数（ここでは、状態ベクトルｘｃ１（ｎ））の時間変化を表している。また、式（３）における観測方程式は、何らかの観測装置を通じて観測する過程を記述したものであり、観測結果（ここでは、観測信号ｙｃ１（ｎ））が、被観測量つまり入力（ここでは、状態ベクトルｘｃ１（ｎ））に依存して時間発展する様子を示している。なお、「時刻ｎにおける状態ベクトルｘｃ１（ｎ）」とは、時刻ｎまでの音声信号（所望信号）からなる状態ベクトルを意味する。【0053】
図５は、従来手法１のアルゴリズムを示す図である。従来手法１のアルゴリズムは、図５に示すように、初期設定（Initialization）の過程と反復（Iteration）の過程とに大別され、反復過程は、ステップ１とステップ２に分かれている。ステップ１では、線形予測アルゴリズムを用いて、クリアな信号（音声信号）ｄ（ｎ）に対するＡＲ係数｛αl（ｎ）｝の推定を行う。また、ステップ２では、１〜５の手順を逐次繰り返す。【0054】
なお、この従来手法１のアルゴリズムは、図６に示す雑音抑圧処理部５０によって実行される。この雑音抑圧処理部５０は、図２に示す本実施の形態における雑音抑圧処理部１５０に対応するものであり、初期設定部５１、ＡＲ係数処理部５２、相関行列演算部５３、カルマンゲインベクトル算出部５４、および最適推定値ベクトル算出部５５を有する。各部５１〜５５の具体的な処理内容は、後述する。【0055】
図７は、図５のアルゴリズムを実行する処理手順を示すフローチャートである。【0056】
まず、初期設定部５１で、初期設定を行う（ＳＴ１０）。具体的には、音声信号からなる状態ベクトルの最適推定値（以下「音声信号の最適推定値ベクトル」という）の初期値ｘｃ１（０｜０）、状態ベクトルの推定誤差（以下「音声信号の推定誤差ベクトル」という）の相関行列の初期値Ｐｃ１（０｜０）、雑音（スカラー）の共分散ｒεｃ１（ｎ）の値、観測遷移ベクトルｍｃ１、時刻ｎのカウンタの初期値、および駆動源ベクトルの共分散Ｒδｃ１（ｎ＋１）［ｉ，ｊ］の値を、次の式（５）に示すようにそれぞれ設定する。【0057】
【数5】
【0058】
ただし、行列Ｉは、単位行列である。また、σｖ２は、雑音εｃ１（ｎ）（＝ｖ（ｎ））の雑音分散であり、既知と仮定している。ここでいう「既知」とは、別の任意の方法（アルゴリズム）で求められて与えられることを意味する。もし雑音εｃ１（ｎ）が白色雑音でありゼロ平均であれば、σｖ２は、次の式（６）で与えられる。ここで、「Ｎ」は、所定のサンプル数である。【0059】
【数6】
【0060】
次に、ＡＲ係数処理部５２で、ＡＲ係数の推定を行う（ＳＴ１２）。具体的には、線形予測アルゴリズムを用いて、クリアな信号（音声信号）ｄ（ｎ＋１）に対するＡＲ係数｛αｌ（ｎ＋１）｝を推定する。【0061】
次に、ＡＲ係数処理部５２で、状態遷移行列の構成を行う（ＳＴ１４）。具体的には、ステップＳＴ１２で推定したＡＲ係数｛αｌ（ｎ＋１）｝を用いて、状態遷移行列Φｃ１（ｎ＋１）を、次の式（７）に示すように構成する。なお、ステップＳＴ１２とステップＳＴ１４が、従来手法１のステップ１に相当する。【0062】
【数7】
【0063】
次に、相関行列演算部５３で、ｎ→（ｎ＋１）の推定誤差の相関値（行列）を計算する（ＳＴ１６）。具体的には、時刻ｎまでの情報により時刻ｎ＋１の状態ベクトルを推定した場合の誤差（音声信号の推定誤差ベクトル）の相関行列Ｐｃ１（ｎ＋１｜ｎ）を計算する。この計算は、ステップＳＴ１０で設定した駆動源ベクトルの共分散Ｒδｃ１（ｎ＋１）［ｉ，ｊ］の値、ステップＳＴ１４で構成した状態遷移行列Φｃ１（ｎ＋１）、およびステップＳＴ１０で設定した（ｎ＝０の場合）または前回のステップＳＴ２６で計算した（ｎ≧１の場合）音声信号の推定誤差ベクトルの相関行列Ｐｃ１（ｎ｜ｎ）を用いて、次の式（８）により行う。なお、このステップＳＴ１６は、従来手法１のステップ２の手順１に相当する。【0064】
【数8】
【0065】
次に、カルマンゲインベクトル算出部５４で、カルマンゲイン（ベクトル）の計算を行う（ＳＴ１８）。具体的には、観測信号の推定誤差にカルマンゲイン（ベクトル）をかけて、時刻ｎまでの情報による時刻ｎ＋１での音声信号の最適推定値ベクトルｘｃ１（ｎ＋１｜ｎ）を加えたものが、時刻ｎ＋１までの情報による当該時刻での音声信号の最適推定値ベクトルｘｃ１（ｎ＋１｜ｎ＋１）になるような、そのカルマンゲインｋｃ１（ｎ＋１）を計算する。この計算は、ステップＳＴ１６で計算した音声信号の推定誤差ベクトルの相関行列Ｐｃ１（ｎ＋１｜ｎ）、ならびに、ステップＳＴ１０でそれぞれ設定した観測遷移ベクトルｍｃ１および雑音の共分散ｒεｃ１（ｎ）の値を用いて、次の式（９）により行う。なお、このステップＳＴ１８は、従来手法１のステップ２の手順２に相当する。【0066】
【数9】
【0067】
次に、最適推定値ベクトル算出部５５で、ｎ→（ｎ＋１）の音声信号の最適推定値（ベクトル）を計算する（ＳＴ２０）。具体的には、時刻ｎまでの情報による時刻ｎ＋１での音声信号の最適推定値ベクトルｘｃ１（ｎ＋１｜ｎ）を計算する。この計算は、前回のステップＳＴ１４で構成した状態遷移行列Φｃ１（ｎ）、および、前回のステップＳＴ２２で計算した音声信号の最適推定値ベクトルｘｃ１（ｎ｜ｎ）を用いて、次の式（１０）により行う。なお、このステップＳＴ２０は、従来手法のステップ２の手順３に相当する。【0068】
【数10】
【0069】
次に、最適推定値ベクトル算出部５５で、（ｎ＋１）→（ｎ＋１）の音声信号の最適推定値（ベクトル）を計算する（ＳＴ２２）。具体的には、時刻ｎ＋１までの情報による当該時刻での音声信号の最適推定値ベクトルｘｃ１（ｎ＋１｜ｎ＋１）を計算する。この計算は、ステップＳＴ２０で計算した音声信号の最適推定値ベクトルｘｃ１（ｎ＋１｜ｎ）、ステップＳＴ１８で計算したカルマンゲイン（ベクトル）ｋｃ１（ｎ＋１）、ステップＳＴ１０で設定した観測遷移ベクトルｍｃ１、および時刻ｎ＋１における観測信号ｙｃ１（ｎ＋１）を用いて、次の式（１１）により行う。なお、このステップＳＴ２２は、従来手法１のステップ２の手順４に相当する。【0070】
【数11】
【0071】
次に、処理を終了するか否かを判断する（ＳＴ２４）。この判断は、例えば、時刻ｎが所定のサンプル数Ｎに達したか否かを判定することによって行う。この判断の結果として、時刻ｎが所定のサンプル数Ｎに達していない場合は（ＳＴ２４：ＮＯ）、ステップＳＴ２６に進み、時刻ｎが所定のサンプル数Ｎに達した場合は（ＳＴ２４：ＹＥＳ）、ステップＳＴ３０に進む。なお、判断の基準は、上記の例に限定されない。例えば、リアルタイムで処理を行う場合は、時刻ｎが所定のサンプル数Ｎに達していなくても、サンプルがなくなった時点で処理を終了するようにしてもよい。【0072】
ステップＳＴ２６では、相関行列演算部５３で、（ｎ＋１）→（ｎ＋１）の推定誤差の相関値（行列）を計算する。具体的には、時刻ｎ＋１までの情報により当該時刻の状態ベクトルを推定した場合の誤差（音声信号の推定誤差ベクトル）の相関行列Ｐｃ１（ｎ＋１｜ｎ＋１）を計算する。この計算は、ステップＳＴ１８で計算したカルマンゲイン（ベクトル）ｋｃ１（ｎ＋１）、ステップＳＴ１０で設定した観測遷移ベクトルｍｃ１、およびステップＳＴ１６で計算した音声信号の推定誤差ベクトルの相関行列Ｐｃ１（ｎ＋１｜ｎ）を用いて、次の式（１２）により行う。なお、このステップＳＴ２６は、従来手法１のステップ２の手順５に相当する。【0073】
【数12】
【0074】
次に、ステップＳＴ２８では、時刻ｎのカウンタを１だけインクリメントして（ｎ＝ｎ＋１）、ステップＳＴ１２に戻る。【0075】
一方、ステップＳＴ３０では、本アルゴリズムの計算結果を出力値として一時保存する。具体的には、ステップＳＴ２２で計算した音声信号の最適推定値ベクトルｘｃ１（ｎ＋１｜ｎ＋１）を、本アルゴリズムの出力値として雑音抑圧処理部５０内に一時保存する。【0076】
図８は、従来手法１とその問題点を視覚的にまとめた説明図である。このように、従来手法１は、ステップ１でＡＲ係数を推定した後に、ステップ２で、ステップ１で推定したＡＲ係数を用いて状態空間モデルを構成することによって、雑音抑圧を実現している。しかし、一般に、実環境では、ステップ１のＡＲ係数の推定においては、ＡＲ係数の次数Ｌｃ１を決定することが非常に困難であるため、不十分な推定結果であるＡＲ係数を用いてカルマンフィルタアルゴリズムを実行しなければならない。このことが雑音抑圧能力の低下を引き起こすことは容易に想像できる。このことは、カルマンフィルタの雑音抑圧能力がＡＲ係数の推定精度に大きく依存していることを意味する。すなわち、従来手法１において最も重要な点は、カルマンフィルタを用いた高性能の雑音抑圧を達成するために、ＡＲ係数の正確な推定を必要とすることである。このことからも、カルマンフィルタの雑音抑圧能力がＡＲ係数の推定精度に大きく依存しており、ＡＲ係数の推定精度が不十分な場合に雑音抑圧能力が大きく劣化することは容易に想像可能である。また、音声信号の音質劣化についても容易に理解できる。【0077】
［従来手法３］次に、従来手法３について説明する。
【0078】
従来手法３では、上記のように、ＡＲ係数の推定を用いないようにするために、新しい状態空間モデルを構成する。すなわち、情報源からのクリアな信号（所望信号）のみを用いて状態方程式を構成し、かつ、そのクリアな信号（所望信号）および雑音を用いて観測方程式を構成する。具体的には、従来手法３では、新しい状態空間モデル（状態方程式と観測方程式）を構成し、この新しい状態空間モデルは、次の式（１３）で書き表される。ただし、式（１３）中のベクトルｘｃ３、δｃ３、ｙｃ３、εｃ３、および行列Φｃ３、Ｍｃ３は、次の式（１４）でそれぞれ定義される。ベクトルｘｃ３は所望信号からなるＬｃ３×１次の状態ベクトル、ベクトルδｃ３はＬｃ３×１次の駆動源ベクトル、ベクトルｙｃ３はＬｃ３×１次の観測信号ベクトル、ベクトルεｃ３はＬｃ３×１次の雑音ベクトル、行列Φｃ３はＬｃ３×Ｌｃ３次の状態遷移行列、行列Ｍｃ３はＬｃ３×Ｌｃ３次の観測遷移行列である。また、図９は、この状態空間モデルにおけるシステム構成図をブロック線図で表したものである。なお、添え字「ｃ３」は、従来手法３に係るものであることを示している。また、「Ｌｃ３」は、状態遷移行列のサイズである。【0079】
【数13】
【0080】
【数14】
【0081】
図９において、「６０」は時刻ｎにおける状態ベクトルｘｃ３（ｎ）、「６１」は時刻ｎ＋１における状態ベクトルｘｃ３（ｎ＋１）、「６２」は時刻ｎにおける観測信号ベクトルｙｃ３（ｎ）、「６３」は時刻ｎにおける雑音ベクトルεｃ３（ｎ）、「６４」は時刻ｎ＋１における駆動源ベクトルδｃ３（ｎ＋１）、「６５」は状態遷移行列Φｃ３、「６６」は観測遷移行列Ｍｃ３である。式（１３）における状態方程式は、推定対象（例えば、音声）のシステムを状態空間モデルで記述したものであり、内部状態つまり状態変数（ここでは、状態ベクトルｘｃ３（ｎ））の時間変化を表している。また、式（１３）における観測方程式は、何らかの観測装置を通じて観測する過程を記述したものであり、観測結果（ここでは、観測信号ベクトルｙｃ３（ｎ））が、被観測量つまり入力（ここでは、状態ベクトルｘｃ３（ｎ））に依存して時間発展する様子を示している。なお、「時刻ｎにおける状態ベクトルｘｃ３（ｎ）」とは、時刻ｎまでの所望信号（例えば、音声信号）からなる状態ベクトルを意味する。【0082】
図１０は、従来手法３のアルゴリズムを示す図である。従来手法３のアルゴリズムは、図１０に示すように、初期設定（Initialization）の過程と反復（Iteration）の過程とに大別され、反復過程は、従来手法１の場合（図５参照）と異なり、新しい状態空間モデル（状態方程式と観測方程式）を構成している。反復過程では、１〜５の手順を逐次繰り返す。【0083】
なお、この従来手法３のアルゴリズムは、図１１に示す雑音抑圧処理部７０によって実行される。この雑音抑圧処理部７０は、図２に示す本実施の形態における雑音抑圧処理部１５０に対応するものであり、初期設定部７２、相関行列演算部７４、重み係数行列算出部７６、および最適推定値ベクトル算出部７８を有する。各部７２〜７８の具体的な処理内容は、後述する。【0084】
図１２は、図１０のアルゴリズムを実行する処理手順を示すフローチャートである。【0085】
まず、初期設定部７２で、初期設定を行う（ＳＴ３０）。具体的には、初期設定部７２において、所望信号（例えば、音声信号）からなる状態ベクトルの最適推定値（以下「所望信号の最適推定値ベクトル」という）の初期値ｘｃ３（０｜０）、状態ベクトルの推定誤差（以下「所望信号の推定誤差ベクトル」という）の相関行列の初期値Ｐｃ３（０｜０）、時刻ｎのカウンタの初期値、状態遷移行列Φｃ３、観測遷移行列Ｍｃ３、雑音ベクトルの共分散Ｒεｃ３（ｎ＋１）［ｉ，ｊ］の値、および駆動源ベクトルの共分散Ｒδｃ３（ｎ＋１）［ｉ，ｊ］の値を、次の式（１５）に示すようにそれぞれ設定する。【0086】
【数15】
【0087】
ただし、行列Ｉは、単位行列である。また、σｖ２は、雑音εｃ３（ｎ）（＝ｖ（ｎ））の雑音分散であり、既知と仮定している。ここでいう「既知」とは、別の任意の方法（アルゴリズム）で求められて与えられることを意味する。もし雑音εｃ３（ｎ）が白色雑音でありゼロ平均であれば、σｖ２は、上記の式（６）で与えられる。【0088】
次に、相関行列演算部７４で、ｎ→（ｎ＋１）の推定誤差の相関値（行列）を計算する（ＳＴ３２）。具体的には、相関行列演算部７４において、時刻ｎまでの情報により時刻ｎ＋１の状態ベクトルを推定した場合の誤差（所望信号の推定誤差ベクトル）の相関行列Ｐｃ３（ｎ＋１｜ｎ）を計算する。この計算は、ステップＳＴ３０でそれぞれ設定した状態遷移行列Φｃ３および駆動源ベクトルの共分散Ｒδｃ３（ｎ＋１）［ｉ，ｊ］の値、ならびに、ステップＳＴ３０で設定した（ｎ＝０の場合）または前回のステップＳＴ４２で計算した（ｎ≧１の場合）所望信号の推定誤差ベクトルの相関行列Ｐｃ３（ｎ｜ｎ）を用いて、次の式（１６）により行う。なお、このステップＳＴ３２は、図１０の反復過程の手順１に相当する。【0089】
【数16】
【0090】
次に、重み係数行列算出部７６で、重み係数（行列）の計算を行う（ＳＴ３４）。具体的には、重み係数行列算出部７６において、観測量である観測信号ベクトルの推定誤差（以下「観測信号の推定誤差ベクトル」という）に重み係数（行列）をかけて、時刻ｎまでの情報による時刻ｎ＋１での所望信号の最適推定値ベクトルｘｃ３（ｎ＋１｜ｎ）を加えたものが、時刻ｎ＋１までの情報による当該時刻での所望信号の最適推定値ベクトルｘｃ３（ｎ＋１｜ｎ＋１）になるような、その重み係数行列Ｋｃ３（ｎ＋１）を計算する。この計算は、ステップＳＴ３２で計算した所望信号の推定誤差ベクトルの相関行列Ｐｃ３（ｎ＋１｜ｎ）、ならびに、ステップＳＴ３０でそれぞれ設定した観測遷移行列Ｍｃ３および雑音ベクトルの共分散Ｒεｃ３（ｎ＋１）［ｉ，ｊ］の値を用いて、次の式（１７）により行う。なお、このステップＳＴ３４は、図１０の反復過程の手順２に相当する。【0091】
【数17】
【0092】
次に、最適推定値ベクトル算出部７８で、ｎ→（ｎ＋１）の状態量（所望信号）の最適推定値（ベクトル）を計算する（ＳＴ３６）。具体的には、最適推定値ベクトル算出部７８において、時刻ｎまでの情報による時刻ｎ＋１での所望信号の最適推定値ベクトルｘｃ３（ｎ＋１｜ｎ）を計算する。この計算は、ステップＳＴ３０で設定した状態遷移行列Φｃ３、および、前回のステップＳＴ３８で計算した所望信号の最適推定値ベクトルｘｃ３（ｎ｜ｎ）を用いて、次の式（１８）により行う。なお、このステップＳＴ３６は、図１０の反復過程の手順３に相当する。【0093】
【数18】
【0094】
次に、同じく最適推定値ベクトル算出部７８で、（ｎ＋１）→（ｎ＋１）の状態量（所望信号）の最適推定値（ベクトル）を計算する（ＳＴ３８）。具体的には、最適推定値ベクトル算出部７８において、時刻ｎ＋１までの情報による当該時刻での所望信号の最適推定値ベクトルｘｃ３（ｎ＋１｜ｎ＋１）を計算する。この計算は、ステップＳＴ３６で計算した所望信号の最適推定値ベクトルｘｃ３（ｎ＋１｜ｎ）、ステップＳＴ３４で計算した重み係数行列Ｋｃ３（ｎ＋１）、ステップＳＴ３０で設定した観測遷移行列Ｍｃ３、および時刻ｎ＋１における観測信号ｙｃ３（ｎ＋１）を用いて、次の式（１９）により行う。なお、このステップＳＴ３８は、図１０の反復過程の手順４に相当する。【0095】
【数19】
【0096】
次に、処理を終了するか否かを判断する（ＳＴ４０）。この判断は、例えば、時刻ｎが所定のサンプル数Ｎに達したか否かを判定することによって行う。この判断の結果として、時刻ｎが所定のサンプル数Ｎに達していない場合は（ＳＴ４０：ＮＯ）、ステップＳＴ４２に進み、時刻ｎが所定のサンプル数Ｎに達した場合は（ＳＴ４０：ＹＥＳ）、ステップＳＴ４６に進む。なお、判断の基準は、上記の例に限定されない。例えば、リアルタイムで処理を行う場合は、時刻ｎが所定のサンプル数Ｎに達していなくても、サンプルがなくなった時点で処理を終了するようにしてもよい。【0097】
ステップＳＴ４２では、相関行列演算部７４で、（ｎ＋１）→（ｎ＋１）の推定誤差の相関値（行列）を計算する。具体的には、相関行列演算部７４において、時刻ｎ＋１までの情報により当該時刻の状態ベクトルを推定した場合の誤差（所望信号の推定誤差ベクトル）の相関行列Ｐｃ３（ｎ＋１｜ｎ＋１）を計算する。この計算は、ステップＳＴ３４で計算した重み係数行列Ｋｃ３（ｎ＋１）、ステップＳＴ３０で設定した観測遷移行列Ｍｃ３、およびステップＳＴ３２で計算した所望信号の推定誤差ベクトルの相関行列Ｐｃ３（ｎ＋１｜ｎ）を用いて、次の式（２０）により行う。なお、このステップＳＴ４２は、図１０の反復過程の手順５に相当する。【0098】
【数20】
【0099】
次に、ステップＳＴ４４では、時刻ｎのカウンタを１だけインクリメントして（ｎ＝ｎ＋１）、ステップＳＴ３２に戻る。【0100】
一方、ステップＳＴ４６では、本アルゴリズムの計算結果を出力値として一時保存する。具体的には、ステップＳＴ３８で計算した所望信号の最適推定値ベクトルｘｃ３（ｎ＋１｜ｎ＋１）を、本アルゴリズムの出力値として雑音抑圧処理部７０内に一時保存する。【0101】
図１３は、従来手法３を視覚的にまとめた説明図である。このように、従来手法３では、新しい状態空間モデルを構成しているため、１段階処理で雑音抑圧が可能となる。なお、雑音抑圧方法の基本的枠組に関する限り、この点は、後述する発明手法１および発明手法２でも同様である。【0102】
ところで、従来手法３のアルゴリズムは、従来手法１とは異なり、駆動源が有色であっても実行することができる。すなわち、従来手法１で述べたように、カルマンフィルタ理論を用いるためには、駆動源δｃ３（ｎ＋１）が白色性であり、かつ、音声信号からなる状態量ｘｃ３（ｎ＋１）と雑音ｖ（ｎ）とが無相関でなければならない。しかし、従来手法３の状態空間モデルにおける状態方程式の駆動源δｃ３（ｎ＋１）は、上記の式（１４）に示すように、有色信号である音声信号ｄ（ｎ＋１）を含んでいる。したがって、一般にカルマンフィルタ理論を適用することはできないものの、従来手法３のアルゴリズムは、駆動源が有色であるにもかかわらず、実行することができる。【0103】
以下では、駆動源が有色信号である従来手法３の有効性、つまり、駆動源が有色であるにもかかわらず従来手法３のアルゴリズムを実行できる理由を説明しておく。もちろん、この理由は、後述する発明手法１および発明手法２にも当てはまる。なお、以下の表記において、行列内を灰色に塗り潰した部分は、駆動源の影響がある部分とし、塗り潰していない部分は、駆動源の影響がない部分とする。【0104】
駆動源δｃ３が有色信号という条件の下で、時刻ｎまでの情報により時刻ｎ＋１の状態ベクトルｘｃ３（ｎ＋１｜ｎ）を推定した場合の誤差（所望信号の推定誤差ベクトル）の相関行列Ｐｃ３（ｎ＋１｜ｎ）は、次の式（２１）と書き表される。【0105】
【数21】
【0106】
ただし、Ｌｃ３×Ｌｃ３次の行列Ｒδｃ３（ｎ＋１）とＬｃ３×Ｌｃ３次の行列Ｑｃ３（ｎ＋１）は、次の式（２２）のように置くことにする。【0107】
【数22】
【0108】
ここで、もし駆動源が白色信号であれば、所望信号の推定誤差ベクトルの相関行列Ｐｃ３（ｎ＋１｜ｎ）は、Ｐｃ３（ｎ＋１｜ｎ）＝Φｃ３Ｐ（ｎ｜ｎ）Φｃ３Ｔ＋Ｒδｃ３（ｎ＋１）となる。このことは、駆動源ベクトルδｃ３（ｎ＋１）と状態ベクトルｘｃ３（ｎ｜ｎ）とは無相関であることを意味する。すなわち、Ｑｃ３Ｔ（ｎ＋１）＝Ｏ（Ｏはゼロ行列）。一方、駆動源が有色信号の場合、駆動源ベクトルδｃ３（ｎ＋１）は、所望信号の推定誤差ベクトルと相関を持つ。すなわち、Ｑｃ３Ｔ（ｎ＋１）≠Ｏとなる。【0109】
以上のことから、次の式（２３）の関係を用いて、行列Ｑｃ３Ｔ（ｎ＋１）の各要素について解析を行う。【0110】
【数23】
【0111】
ただし、Ｌｃ３×Ｌｃ３次の行列Γは、次の式（２４）となる。【0112】
【数24】
【0113】
上記の式（２３）を用いて、行列Ｑｃ３（ｎ＋１）を変形すると、次の式（２５）のように書き表される。【0114】
【数25】
【0115】
行列Ｑｃ３（ｎ＋１）の各要素を明らかにするために、次の式（２６）を定義すると、行列Ｑｃ３（ｎ＋１）の各要素は、次の式（２７）および式（２８）となる。【0116】
【数26】
【0117】
【数27】
【0118】
【数28】
【0119】
ただし、行列Ｑｃ３（ｎ＋１）の１列目の要素｛ｅｉ（ｎ）｝は、次の式（２９）となる。【0120】
【数29】
【0121】
ここで、行列Ｑｃ３（ｎ＋１）の１列目の要素がすべてゼロ、つまり、｛ｅｉ（ｎ）｝＝０であれば、駆動源の影響を受けていない。しかし、上記の式（２８）の行列Ｑｃ３（ｎ＋１）の１列目の各要素はゼロではない、つまり、｛ｅｉ（ｎ）｝≠０であるため、駆動源の影響を受けてしまっている。【0122】
一方、もし上記の式（２８）の行列Ｑｃ３（ｎ＋１）の１列目の各要素｛ｅｉ（ｎ）｝を何らかの方法で得ることができれば、｛ｅｉ（ｎ）｝の位置は、更新と無関係に常に一定であるため、何らかの方法で得られた｛ｅｉ（ｎ）｝を差し引きすることで、駆動源の影響を取り除くことが可能となる。このことは、行列Ｑｃ３（ｎ＋１）をゼロ行列と置くことに等しい。したがって、従来手法３の状態空間モデルであれば、駆動源が有色であっても実行可能であるということができる。すなわち、従来手法３のアルゴリズムは、駆動源が有色であるにもかかわらず、実行することができる。【0123】
［発明手法］次に、発明手法について説明する。
【0124】
発明手法では、演算量のさらなる低減と雑音抑圧能力のさらなる向上とを図るために、従来手法３の状態空間モデルを基礎にして、さらに新しい状態空間モデルを構成する。すなわち、発明手法では、情報源からのクリアな信号（所望信号）のみを用いて状態方程式を構成し、かつ、ある１つの時刻のみの観測信号に対して、そのクリアな信号（所望信号）および雑音を用いて観測方程式を構成する。なお、以下では、発明手法として２つの具体的な手法を提示する。ここでは、便宜上、第１の発明手法および第２の発明手法を、発明手法１および発明手法２とそれぞれ称し、順次説明する。【0125】
＜発明手法１＞図１４および図１５は、従来手法３の問題点を視覚的にまとめた説明図である。
【0126】
まず、上記のように、従来手法３の状態空間モデルでは、観測信号ｒ（ｎ）をベクトル化｛ｙｃ３（ｎ）｝することにより、多くの過去情報を用いて雑音抑圧を実行している。例えば、図１４Ａに示すように、従来手法３では、「おはよう」における未来の「う」を推定する場合のように、予測の推定精度の観点から観測信号ｒ（ｎ）のデータ量を多く使用している。具体的には、図１４Ａに示すように、「おはよう」における未来の「う」を推定する場合、１時刻前の「よ」を知っていても未来の「う」は推測できない。また、２時刻前までの「は」と「よ」を知っていても「はよ」から未来の「う」の推定は困難である。一方、３時刻前までの「お」と「は」と「よ」を知っていれば「おはよ」から未来の「う」は推測できる。それゆえ、予測の推定精度の観点から、観測信号ｒ（ｎ）のデータ量は多くなる。すなわち、過去のデータを多く使用したほうが、未来を推定しやすくなる。【0127】
しかし、図１４Ｂに示すように、過去の観測信号ｒ（ｎ）にも雑音が混入するため（つまり、誤差あり）、多くの過去情報を用いた観測信号ベクトルは雑音（誤差）を含んだものとなる。この場合、予測の推定精度を上げるために過去の観測信号ｒ（ｎ）が多く必要であるという命題は、必ずしも正しくない。したがって、従来手法３においては、所望信号の推定精度（つまり、雑音抑圧能力）の向上に一定の限界がある。【0128】
また、図１５に強調的に示すように、従来手法３のアルゴリズムでは、逆行列の計算が必要であるため（図１０の反復過程の手順２、つまり、上記の式（１７）参照）、従来手法１に対する演算量の軽減効果にも一定の限界がある。【0129】
そこで、発明手法２では、従来手法３に対してさらに演算量の軽減および雑音抑圧能力の向上を図るために、さらに新しい状態空間モデル（状態方程式と観測方程式）を、次の式（３０）のように構成する。ただし、式（３０）中のベクトルｘｐ１、δｐ１、ｍｐ１、スカラーｙｐ１、εｐ１、および行列Φｐ１は、次の式（３１）でそれぞれ定義される。ベクトルｘｐ１は所望信号からなるＬｐ１×１次の状態ベクトル、ベクトルδｐ１はＬｐ１×１次の駆動源ベクトル、スカラーｙｐ１は観測信号、スカラーεｐ１は雑音、行列Φｐ１はＬｐ１×Ｌｐ１次の状態遷移行列、ベクトルｍｐ１はＬｐ１×１次の観測遷移ベクトルである。なお、添え字「ｐ１」は、発明手法１に係るものであることを示している。また、「Ｌｐ１」は、状態遷移行列のサイズである。【0130】
【数30】
【0131】
【数31】
【0132】
図１６は、このさらに新しい状態空間モデルを視覚的に説明するための概略図である。図１６に示すように、例えば、音声信号の場合、ある１つの時刻ｎにおける観測信号ｒ（ｎ）は、その時刻ｎにおける音声信号ｄ（ｎ）と雑音ｖ（ｎ）を加算した信号である（図３および式（１）参照）。音声信号の場合、状態方程式は、気道１７０の構造（つまり、音声の構造）を記述したものであり、観測方程式は、任意の観測装置を通じて、ある１つの時刻ｎにおいて、気道１７０から発せられた音声信号ｄ（ｎ）と、この音声信号ｄ（ｎ）以外の雑音ｖ（ｎ）とが混ざり合って、観測される過程を記述したものである。なお、音声信号の場合には、状態方程式を「音声構造方程式」、観測方程式を「音声観測方程式」とそれぞれ呼ぶこともできる。【0133】
図１７は、この状態空間モデルにおけるシステム構成図をブロック線図で表したものである。図１７において、「２００」は時刻ｎにおける状態ベクトルｘｐ１（ｎ）、「２０１」は時刻ｎ＋１における状態ベクトルｘｐ１（ｎ＋１）、「２０２」は時刻ｎにおける観測信号ｙｐ１（ｎ）、「２０３」は時刻ｎにおける雑音εｐ１（ｎ）、「２０４」は時刻ｎ＋１における駆動源ベクトルδｐ１（ｎ＋１）、「２０５」は状態遷移行列Φｐ１、「２０６」は観測遷移ベクトルｍｐ１である。式（３０）における状態方程式は、推定対象（例えば、音声）のシステムを状態空間モデルで記述したものであり、内部状態つまり状態変数（ここでは、状態ベクトルｘｐ１（ｎ））の時間変化を表している。また、式（３０）における観測方程式は、何らかの観測装置を通じて観測する過程を記述したものであり、観測結果（ここでは、観測信号ｙｐ１（ｎ））が、被観測量つまり入力（ここでは、状態ベクトルｘｐ１（ｎ））に依存して時間発展する様子を示している。なお、「時刻ｎにおける状態ベクトルｘｐ１（ｎ）」とは、時刻ｎまでの所望信号（例えば、音声信号）からなる状態ベクトルを意味する。【0134】
図１８は、発明手法１の状態方程式を視覚的に示す説明図であり、図１９および図２０は、発明手法１における新しい観測方程式の構成を視覚的に示す説明図である。すなわち、発明手法１は、従来手法３と比較して、図１８に示すように、状態方程式は、状態量の個数の点を除き、同じであるが（図１３参照）、図１９および図２０に示すように、観測方程式を新しく構成して演算量の軽減を図っている。具体的には、発明手法１では、従来手法３の観測方程式から、状態方程式の駆動源の影響がないように、任意の一時刻の部分のみを抽出して、従来手法３と異なる新しい観測方程式を構成している（特に図１９参照）。これにより、発明手法１では、例えば、従来手法３の観測方程式における観測信号ベクトルｙｃ３（ｎ＋１）、観測遷移行列Ｍｃ３、および雑音ベクトルεｃ３（ｎ＋１）が、それぞれ、スカラー化、ベクトル化、スカラー化されて、観測信号ｙｐ１（ｎ＋１）、観測遷移ベクトルｍｐ１、および雑音εｐ１（ｎ＋１）となっている（特に図２０参照）。この結果、発明手法１は、従来手法３よりも演算量が大幅に軽減される。なお、演算量の軽減効果については、後で詳述する。【0135】
なお、状態方程式に関して、発明手法１では、従来手法３よりも、状態量の個数、つまり、状態遷移行列のサイズを大きくとる。すなわち、発明手法１の状態遷移行列のサイズＬｐ１は、従来手法３の状態遷移行列のサイズＬｃ３よりも大きい（Ｌｐ１＞Ｌｃ３）。これは、例えば、音声信号の場合、気道１７０の構造を明らかにするためには、気道１７０のデータは多ければ多いほど良いためである（図１６参照）。一方、観測方程式に関しては、発明手法１では、上記のように、従来手法３の観測方程式から任意の一時刻の部分のみを抽出して従来手法３の観測方程式と異なる新しい観測方程式を構成している。すなわち、発明手法１では、例えば、観測信号として、過去のデータを使用する従来手法３とは異なり、現在のデータ（瞬時データ）のみを使用する。この点で、発明手法１は、本質的に、従来手法３とは異なる構成を採っている。【0136】
図２１は、発明手法１のアルゴリズムを示す図である。発明手法１のアルゴリズムは、図２１に示すように、初期設定（Initialization）の過程と反復（Iteration）の過程とに大別され、反復過程は、従来手法３の場合（図１０参照）よりも演算量が軽減されるように、さらに新しい状態空間モデル（状態方程式と観測方程式）を構成している。反復過程では、１〜５の手順を逐次繰り返す。【0137】
図２２は、発明手法１のアルゴリズムの利点を視覚的にまとめた説明図である。発明手法１では、図２２に示すように、従来手法３のアルゴリズムにおける重み係数行列Ｋｃ３（ｎ＋１）、観測遷移行列Ｍｃ３、および観測信号ベクトルｙｃ３（ｎ＋１）が、それぞれ、ベクトル化、ベクトル化、スカラー化されて、重み係数ベクトルｋｐ１（ｎ＋１）、観測遷移ベクトルｍｐ１、および観測信号ｙｐ１（ｎ＋１）となっている。これにより、発明手法１では、従来手法３における逆行列の計算（図１０の反復過程の手順２参照）が、スカラーの逆数の計算（図２１の反復過程の手順２参照）となる。一般に、逆行列の計算は複雑であるが、スカラーの逆数の計算は簡単である。したがって、発明手法１は、従来手法３よりも演算量が大幅に軽減される。【0138】
図２３は、図２１のアルゴリズムを実行する処理手順を示すフローチャートである。【0139】
まず、初期設定部１５２で、初期設定を行う（ＳＴ１０００）。具体的には、初期設定部１５２において、所望信号（例えば、音声信号）からなる状態ベクトルの最適推定値（以下「所望信号の最適推定値ベクトル」という）の初期値ｘｐ１（０｜０）、状態ベクトルの推定誤差（以下「所望信号の推定誤差ベクトル」という）の相関行列の初期値Ｐｐ１（０｜０）、雑音（スカラー）の共分散ｒεｐ１（ｎ＋１）の値、観測遷移ベクトルｍｐ１、時刻ｎのカウンタの初期値、状態遷移行列Φｐ１、および駆動源ベクトルの共分散Ｒδｐ１（ｎ＋１）［ｉ，ｊ］の値を、次の式（３２）に示すようにそれぞれ設定する。
【0140】
【数32】
【0141】
ただし、行列Ｉは、単位行列である。また、σｖ２は、雑音εｐ１（ｎ）（＝ｖ（ｎ））の雑音分散であり、既知と仮定している。ここでいう「既知」とは、別の任意の方法（アルゴリズム）で求められて与えられることを意味する。もし雑音εｐ１（ｎ）が白色雑音でありゼロ平均であれば、σｖ２は、上記の式（６）で与えられる。【0142】
次に、第１の相関演算部としての相関演算部１５４で、ｎ→（ｎ＋１）の推定誤差の相関値（ベクトル）を計算する（ＳＴ１１００）。具体的には、相関演算部１５４において、時刻ｎのみの観測信号に対して、時刻ｎまでの情報により時刻ｎ＋１の状態ベクトルを推定した場合の誤差（所望信号の推定誤差ベクトル）の相関行列Ｐｐ１（ｎ＋１｜ｎ）を計算する。この計算は、ステップＳＴ１０００でそれぞれ設定した状態遷移行列Φｐ１および駆動源ベクトルの共分散Ｒδｐ１（ｎ＋１）［ｉ，ｊ］の値、ならびに、ステップＳＴ１０００で設定した（ｎ＝０の場合）または前回のステップＳＴ１６００で計算した（ｎ≧１の場合）所望信号の推定誤差ベクトルの相関行列Ｐｐ１（ｎ｜ｎ）を用いて、次の式（３３）により行う。なお、このステップＳＴ１１００は、図２１の反復過程の手順１に相当する。【0143】
【数33】
【0144】
次に、重み係数算出部１５６で、重み係数（ベクトル）の計算を行う（ＳＴ１２００）。具体的には、重み係数算出部１５６において、時刻ｎのみの観測信号に対して、観測量である観測信号の推定誤差（スカラー）に重み係数（行列）をかけて、時刻ｎまでの情報による時刻ｎ＋１での所望信号の最適推定値ベクトルｘｐ１（ｎ＋１｜ｎ）を加えたものが、時刻ｎ＋１までの情報による当該時刻での所望信号の最適推定値ベクトルｘｐ１（ｎ＋１｜ｎ＋１）になるような、その重み係数ベクトルｋｐ１（ｎ＋１）を計算する。この計算は、ステップＳＴ１１００で計算した所望信号の推定誤差ベクトルの相関行列Ｐｐ１（ｎ＋１｜ｎ）、ならびに、ステップＳＴ１０００でそれぞれ設定した観測遷移ベクトルｍｐ１および雑音の共分散ｒεｐ１（ｎ＋１）の値を用いて、次の式（３４）により行う。なお、このステップＳＴ１２００は、図２１の反復過程の手順２に相当する。【0145】
【数34】
【0146】
次に、第１の最適推定値算出部としての最適推定値算出部１５８で、ｎ→（ｎ＋１）の状態量（所望信号）の最適推定値（ベクトル）を計算する（ＳＴ１３００）。具体的には、最適推定値算出部１５８において、時刻ｎのみの観測信号に対して、時刻ｎまでの情報による時刻ｎ＋１での所望信号の最適推定値ベクトルｘｐ１（ｎ＋１｜ｎ）を計算する。この計算は、ステップＳＴ１０００で設定した状態遷移行列Φｐ１、および、前回のステップＳＴ１４００で計算した所望信号の最適推定値ベクトルｘｐ１（ｎ｜ｎ）を用いて、次の式（３５）により行う。なお、このステップＳＴ１３００は、図２１の反復過程の手順３に相当する。【0147】
【数35】
【0148】
次に、第２の最適推定値算出部としての最適推定値算出部１５８で、（ｎ＋１）→（ｎ＋１）の状態量（所望信号）の最適推定値（ベクトル）を計算する（ＳＴ１４００）。具体的には、最適推定値算出部１５８において、時刻ｎのみの観測信号に対して、時刻ｎ＋１までの情報による当該時刻での所望信号の最適推定値ベクトルｘｐ１（ｎ＋１｜ｎ＋１）を計算する。この計算は、ステップＳＴ１３００で計算した所望信号の最適推定値ベクトルｘｐ１（ｎ＋１｜ｎ）、ステップＳＴ１２００で計算した重み係数ベクトルｋｐ１（ｎ＋１）、ステップＳＴ１０００で設定した観測遷移ベクトルｍｐ１、および時刻ｎ＋１における観測信号ｙｐ１（ｎ＋１）を用いて、次の式（３６）により行う。なお、このステップＳＴ１４００は、図２１の反復過程の手順４に相当する。【0149】
【数36】
【0150】
次に、処理を終了するか否かを判断する（ＳＴ１５００）。この判断は、例えば、時刻ｎが所定のサンプル数Ｎに達したか否かを判定することによって行う。この判断の結果として、時刻ｎが所定のサンプル数Ｎに達していない場合は（ＳＴ１５００：ＮＯ）、ステップＳＴ１６００に進み、時刻ｎが所定のサンプル数Ｎに達した場合は（ＳＴ１５００：ＹＥＳ）、ステップＳＴ１８００に進む。なお、判断の基準は、上記の例に限定されない。例えば、リアルタイムで処理を行う場合は、時刻ｎが所定のサンプル数Ｎに達していなくても、サンプルがなくなった時点で処理を終了するようにしてもよい。【0151】
ステップＳＴ１６００では、第２の相関演算部としての相関演算部１５４で、（ｎ＋１）→（ｎ＋１）の推定誤差の相関値（ベクトル）を計算する。具体的には、相関演算部１５４において、時刻ｎのみの観測信号に対して、時刻ｎ＋１までの情報により当該時刻の状態ベクトルを推定した場合の誤差（所望信号の推定誤差ベクトル）の相関行列Ｐｐ１（ｎ＋１｜ｎ＋１）を計算する。この計算は、ステップＳＴ１２００で計算した重み係数ベクトルｋｐ１（ｎ＋１）、ステップＳＴ１０００で設定した観測遷移ベクトルｍｐ１、およびステップＳＴ１１００で計算した所望信号の推定誤差ベクトルの相関行列Ｐｐ１（ｎ＋１｜ｎ）を用いて、次の式（３７）により行う。なお、このステップＳＴ１６００は、図１９の反復過程の手順５に相当する。【0152】
【数37】
【0153】
次に、ステップＳＴ１７００では、時刻ｎのカウンタを１だけインクリメントして（ｎ＝ｎ＋１）、ステップＳＴ１１００に戻る。【0154】
一方、ステップＳＴ１８００では、本アルゴリズムの計算結果を出力値として一時保存する。具体的には、ステップＳＴ１４００で計算した所望信号の最適推定値ベクトルｘｐ１（ｎ＋１｜ｎ＋１）を、本アルゴリズムの出力値として雑音抑圧処理部１５０内に一時保存する。【0155】
このように、発明手法１においても、従来手法３と同様に、ＡＲ係数の推定を必要としない新しい状態空間モデルを構成しているため、１段階処理で雑音抑圧が可能となる。【0156】
また、従来手法３と同様に、発明手法１のアルゴリズムは、従来手法１とは異なり、駆動源が有色であっても実行することができる。すなわち、従来手法１で述べたように、カルマンフィルタ理論を用いるためには、駆動源δｐ１（ｎ＋１）が白色性であり、かつ、音声信号（所望信号）からなる状態量ｘｐ１（ｎ＋１）と雑音ｖ（ｎ）とが無相関でなければならない。しかし、発明手法１の状態空間モデルにおける状態方程式の駆動源δｐ１（ｎ＋１）は、上記の式（３１）に示すように、有色信号である音声信号ｄ（ｎ＋１）を含んでいる。したがって、一般にカルマンフィルタ理論を適用することはできないものの、発明手法１のアルゴリズムは、駆動源が有色であるにもかかわらず、実行することができる。【0157】
以下では、駆動源が有色信号である発明手法１の有効性、つまり、駆動源が有色であるにもかかわらず発明手法１のアルゴリズムを実行できる理由を説明する。なお、以下の表記において、行列内を灰色に塗り潰した部分は、駆動源の影響がある部分とし、塗り潰していない部分は、駆動源の影響がない部分とする。また、更新回数ｎにおける各手順において、Ｐｐ１（ｎ｜ｎ）、ｘｐ１（ｎ｜ｎ）は、Ｑｐ１（ｎ＋１）の影響を最大限に見積るために、これらのすべての要素に行列Ｑｐ１（ｎ＋１）の影響が含まれていることを仮定する。すなわち、Ｐｐ１（ｎ｜ｎ）、ｘｐ１（ｎ｜ｎ）のすべての要素は灰色で表される。【0158】
駆動源δｐ１が有色信号という条件の下で、時刻ｎまでの情報により時刻ｎ＋１の状態ベクトルｘｐ１（ｎ＋１｜ｎ）を推定した場合の誤差（所望信号の推定誤差ベクトル）の相関行列Ｐｐ１（ｎ＋１｜ｎ）は、次の式（３８）と書き表される。【0159】
【数38】
【0160】
ただし、Ｌｐ１×Ｌｐ１次の行列Ｒδｐ１（ｎ＋１）とＬｐ１×Ｌｐ１次の行列Ｑｐ１（ｎ＋１）は、次の式（３９）のように置くことにする。【0161】
【数39】
【0162】
行列Ｐｐ１（ｎ＋１｜ｎ）に対して、駆動源の影響を受ける部分は、次の式（４０）に示すものとなる。【0163】
【数40】
【0164】
ただし、｛ｅｉ（ｎ）｝は、次の式（４１）である。【0165】
【数41】
【0166】
行列Ｐｐ１（ｎ＋１｜ｎ）の結果を用いて、ベクトルｋｐ１（ｎ＋１）に対する駆動源の影響部分は、次の式（４２）に示すものとなる。【0167】
【数42】
【0168】
また、次の式（４３）を定義すると、最適推定値ベクトルｘｐ１（ｎ＋１｜ｎ）に対する駆動源の影響部分は、次の式（４４）に示すものとなる。【0169】
【数43】
【0170】
【数44】
【0171】
したがって、最適推定値ベクトルｘｐ１（ｎ＋１｜ｎ＋１）に対する駆動源の影響部分は、次の式（４５）のように書き表される。【0172】
【数45】
【0173】
復元された信号は、最適推定値ベクトルｘｐ１（ｎ＋１｜ｎ＋１）の１行１列目の要素、つまり、最適推定値ｘｐ１（ｎ＋１｜ｎ＋１）［１，１］であるため、発明手法１は、有色性の駆動源に対して影響はないといえる。【0174】
同様に、所望信号の誤差推定ベクトルの相関行列Ｐｐ１（ｎ＋１｜ｎ＋１）に対する駆動源の影響部分は、次の式（４６）のように書き表される。【0175】
【数46】
【0176】
以上のことから、復元された信号は、有色駆動源の影響を受けないことが明らかとなった。また、これらの議論は、ｎ回目の更新に限らず、すべての更新回数においても成立する。それゆえ、新に発明した有色駆動源を含んだ状態空間モデルの有効性は明らかである。【0177】
＜発明手法２＞上記のように、発明手法１は、従来手法３と比べて、非常に少ない演算量で雑音抑圧を実現することができる。しかし、発明手法１のアルゴリズムを精査してみると、発明手法１のアルゴリズムは、今まだ余分な計算を含んでいる、つまり、計算する必要のない部分についても演算を行っている。このことは、演算量の増加や不必要な演算誤差などを引き起こす原因になる可能性がある。
【0178】
そこで、発明手法２では、所望信号の推定に必要な部分のみ、つまり、取り出す推定信号のみに着目して、発明手法１のアルゴリズムの一部の演算を省略するようにしている。具体的には、発明手法２では、復元された所望信号（例えば、音声信号）、つまり、所望信号の最適推定値ｘｐ１（ｎ＋１｜ｎ＋１）［１，１］のみに着目することにより、発明手法１に対して、さらなる演算量の軽減を図るようにしている。【0179】
すなわち、復元された所望信号ｘｐ１（ｎ＋１｜ｎ＋１）［１，１］に着目すると、これは、次の式（４７）および式（４８）と書き表される。【0180】
【数47】
【0181】
【数48】
【0182】
このとき、ｙｐ１（ｎ＋１）は観測信号であり既知であるため、復元された所望信号ｘｐ１（ｎ＋１｜ｎ＋１）［１，１］を得るためには、重み係数ベクトルｋｐ１（ｎ＋１）の１行１列目の要素のみ必要となる。【0183】
重み係数ベクトルｋｐ１（ｎ＋１）は、下記の式（４９）と書き表されるため、重み係数ベクトルｋｐ１（ｎ＋１）の１行１列目の要素は、次の式（５０）で得ることが可能である。【0184】
【数49】
【0185】
【数50】
【0186】
同様に、重み係数ｋｐ１（ｎ＋１）［１，１］を得るためには、所望信号の推定誤差ベクトルの相関行列Ｐｐ１（ｎ＋１｜ｎ）の１行１列目の要素が必要となる。したがって、この行列Ｐｐ１（ｎ＋１｜ｎ）は、下記の式（５１）と書き表されるため、この行列Ｐｐ１（ｎ＋１｜ｎ）の１行１列目の要素は、次の式（５２）により得られる。【0187】
【数51】
【0188】
【数52】
【0189】
図２４〜図２６は、発明手法２の導出方法を視覚的に示す段階別説明図である。上記した発明手法２の導出方法を視覚的に説明すると、図２４〜図２６に示す通りである。【0190】
以上のことから、発明手法１から余分な演算を取り除いた発明手法２のアルゴリズムは、図２７で与えられる。なお、添え字「ｐ２」は、発明手法２に係るものであることを示している。また、図２８は、発明手法２のアルゴリズムの利点を視覚的にまとめた説明図である。【0191】
発明手法２のアルゴリズムは、図２７に示すように、初期設定（Initialization）の過程と反復（Iteration）の過程とに大別され、反復過程は、従来手法３や発明手法１の場合（図１０、図２１参照）と同様にＡＲ係数の推定を用いないように、かつ、従来手法３や発明手法１の場合（図１０、図２１参照）よりも演算量が軽減されるように、処理手順が構成されている。反復過程では、１〜３の手順を逐次繰り返す。すなわち、発明手法２の反復過程の手順の数は、従来手法３や発明手法１と比較して、５ステップから３ステップに低減されている。なお、発明手法２は、発明手法１の演算の一部を省略したものであり、その状態空間モデル（状態方程式と観測方程式）は、発明手法１と同様である。【0192】
図２９は、図２７のアルゴリズムを実行する処理手順を示すフローチャートである。【0193】
まず、初期設定部１５２で、初期設定を行う（ＳＴ２０００）。具体的には、初期設定部１５２において、雑音（スカラー）の共分散ｒεｐ２（ｎ＋１）の値、時刻ｎのカウンタの初期値、および駆動源（スカラー）の共分散ｒδｐ２（ｎ＋１）の値を、次の式（５３）に示すようにそれぞれ設定する。【0194】
【数53】
【0195】
ただし、σｖ２は、雑音εｐ２（ｎ）（＝ｖ（ｎ））の雑音分散であり、既知と仮定している。ここでいう「既知」とは、別の任意の方法（アルゴリズム）で求められて与えられることを意味する。もし雑音εｐ２（ｎ）が白色雑音でありゼロ平均であれば、σｖ２は、上記の式（６）で与えられる。【0196】
次に、相関演算部１５４で、ｎ→（ｎ＋１）の推定誤差の相関値（スカラー）を計算する（ＳＴ２１００）。具体的には、相関演算部１５４において、時刻ｎのみの観測信号に対して、時刻ｎまでの情報により時刻ｎ＋１の状態ベクトルを推定した場合の誤差（所望信号の推定誤差ベクトル）の相関行列Ｐｐ２（ｎ＋１｜ｎ）の１行１列目の要素、つまり、スカラーであるＰｐ２（ｎ＋１｜ｎ）［１，１］を計算する。この計算は、ステップＳＴ２０００で設定した駆動源（スカラー）の共分散ｒδｐ２（ｎ＋１）の値を用いて、次の式（５４）により行う。なお、このステップＳＴ２１００は、図２７の反復過程の手順１に相当する。【0197】
【数54】
【0198】
次に、重み係数算出部１５６で、重み係数（スカラー）の計算を行う（ＳＴ２２００）。具体的には、重み係数算出部１５６において、時刻ｎのみの観測信号に対して、観測量である観測信号の推定誤差（スカラー）に重み係数（スカラー）をかけて、時刻ｎまでの情報による時刻ｎ＋１での所望信号の最適推定値ベクトルｘｐ２（ｎ＋１｜ｎ）を加えたものが、時刻ｎ＋１までの情報による当該時刻での所望信号の最適推定値ベクトルｘｐ２（ｎ＋１｜ｎ＋１）になるような、その重み係数ベクトルｋｐ２（ｎ＋１）の１行１列目の要素、つまり、スカラーであるｋｐ２（ｎ＋１）［１，１］を計算する。この計算は、ステップＳＴ２１００で計算した所望信号の推定誤差ベクトルの相関行列Ｐｐ２（ｎ＋１｜ｎ）の１行１列目の要素であるＰｐ２（ｎ＋１｜ｎ）［１，１］、および、ステップＳＴ２０００で設定した雑音の共分散ｒεｐ２（ｎ＋１）の値を用いて、次の式（５５）により行う。なお、このステップＳＴ２２００は、図２５の反復過程の手順２に相当する。【0199】
【数55】
【0200】
次に、最適推定値算出部１５８で、（ｎ＋１）→（ｎ＋１）の状態量（所望信号）の最適推定値（スカラー）を計算する（ＳＴ２３００）。具体的には、最適推定値算出部１５８において、時刻ｎのみの観測信号に対して、時刻ｎ＋１までの情報による当該時刻での所望信号の最適推定値ベクトルのｘｐ２（ｎ＋１｜ｎ＋１）の１行１列目の要素、つまり、所望信号の最適推定値ｘｐ２（ｎ＋１｜ｎ＋１）［１，１］を計算する。この計算は、ステップＳＴ２２００で計算した重み係数ｋｐ２（ｎ＋１）［１，１］、および、時刻ｎ＋１における観測信号ｙｐ２（ｎ＋１）を用いて、次の式（５６）により行う。なお、このステップＳＴ２３００は、図２７の反復過程の手順３に相当する。【0201】
【数56】
【0202】
次に、処理を終了するか否かを判断する（ＳＴ２４００）。この判断は、例えば、時刻ｎが所定のサンプル数Ｎに達したか否かを判定することによって行う。この判断の結果として、時刻ｎが所定のサンプル数Ｎに達していない場合は（ＳＴ２４００：ＮＯ）、ステップＳＴ２５００に進み、時刻ｎが所定のサンプル数Ｎに達した場合は（ＳＴ２４００：ＹＥＳ）、ステップＳＴ２６００に進む。なお、判断の基準は、上記の例に限定されない。例えば、リアルタイムで処理を行う場合は、時刻ｎが所定のサンプル数Ｎに達していなくても、サンプルがなくなった時点で処理を終了するようにしてもよい。【0203】
ステップＳＴ２５００では、時刻ｎのカウンタを１だけインクリメントして（ｎ＝ｎ＋１）、ステップＳＴ２１００に戻る。【0204】
一方、ステップＳＴ２６００では、本アルゴリズムの計算結果を出力値として一時保存する。具体的には、ステップＳＴ２３００で計算した所望信号の最適推定値ｘｐ２（ｎ＋１｜ｎ＋１）［１，１］を、本アルゴリズムの出力値として雑音抑圧処理部１５０内に一時保存する。【0205】
このように、発明手法２においても、従来手法３や発明手法１と同様に、ＡＲ係数の推定を必要としない新しい状態空間モデルを構成しているため、１段階処理で雑音抑圧が可能となる。【0206】
また、発明手法２では、発明手法１と比較して、さらなるスカラー化が図られ、反復過程の手順の数が５ステップから３ステップに低減されているため、演算量のさらなる軽減を実現することができる（図２８参照）。【0207】
また、発明手法２では、発明手法１と比較して、観測遷移ベクトルｍや状態遷移行列Φを設定する必要がないため（図２７の初期設定、式（５３）参照）、反復過程の手順の数の低減に伴う演算量の軽減と相俟って、メモリ容量の大幅な削減を図ることができる。【0208】
次に、発明手法１および発明手法２の演算量の削減効果について説明する。【0209】
図３０は、上記各手法の演算量を説明するための図であり、特に、図３０Ａは、従来手法１の場合、図３０Ｂは、従来手法３の場合、図３０Ｃは、発明手法１の場合、図３０Ｄは、発明手法２の場合をそれぞれ示している。また、図３１および図３２は、次数と演算量との関係を示すグラフであり、特に、図３１は、主に従来手法１と従来手法３および発明手法１〜２との比較を示すグラフ図であり、図３２は、従来手法３と発明手法１〜２との比較を示すグラフ図である。【0210】
図３０〜図３２において、演算量は、乗算の回数によって表される。図中の「Ｌｃ１」は、ＡＲ係数の次数であり、「Ｌｃ３」および「Ｌｐ」は、状態遷移行列のサイズである。ここでは、便宜上、Ｌｃ１とＬｃ３をまとめて「Ｌｃ」と表記する。また、「Ｌｐ」は、発明手法１および発明手法２における状態遷移行列のサイズＬｐ１、Ｌｐ２をまとめて表記したものである。例えば、従来手法１の場合、図３０Ａに示すように、アルゴリズムの反復過程の手順１〜５の演算量は、順に、２Ｌｃ１２、Ｌｃ１、Ｌｃ１、Ｌｃ１、Ｌｃ１２であり、その合計は、３Ｌｃ１２＋３Ｌｃ１である。また、従来手法３の場合、図３０Ｂに示すように、アルゴリズムの反復過程の手順１〜５の演算量は、順に、０、Ｌｃ３−１、０、Ｌｃ３−１、Ｌｃ３−１であり、その合計は、３Ｌｃ３−３である。また、発明手法１の場合、図３０Ｃに示すように、アルゴリズムの反復過程の手順１〜５の演算量は、順に、０、１、０、１、１であり、その合計は、３である。また、発明手法２の場合、図３０Ｄに示すように、アルゴリズムの反復過程の手順１〜３の演算量は、順に、０、１、１であり、その合計は、２である。【0211】
したがって、従来手法１では、ＡＲ係数の次数Ｌｃ１の２乗に比例して演算量が著しく増大し（図３１参照）、従来手法３では、状態遷移行列のサイズＬｃ３に比例して演算量が増大するのに対し（図３２参照）、発明手法１および発明手法２では、状態遷移行列のサイズＬｐに関係なくそれぞれ演算量が一定であることがわかる（図３２参照）。すなわち、演算量の大小に関して、従来手法１＞従来手法３＞発明手法１＞発明手法２である。特に、従来手法１と従来手法３および発明手法１〜２とを比較した場合、図３１に示すように、従来手法３および発明手法１〜２は、従来手法１に比べて演算量が大幅に軽減されている。さらに、従来手法３と発明手法１〜２とを比較した場合、図３２に示すように、発明手法１および発明手法２は、従来手法３に比べて演算量がさらに大幅に軽減されている。また、図示しないが、このような演算量の削減によって、上記のように、メモリ容量も大幅に削減される。【0212】
また、本発明者は、本発明の効果（発明手法１〜２の有効性）を実証するために実験を行った。具体的には、発明手法１〜２の雑音抑圧能力を評価するために、音声信号を用いて、客観的評価と主観的評価を行った。客観的評価は、（１）波形による評価（音声波形のシミュレーション）と（２）数値による評価である。前者は、いわば見た目による評価であり、後者は、能力（雑音抑圧能力）の評価である。また、主観的評価は、（３）聞き取り調査である。これは、いわば音質の評価である。これらの実験では、本発明が、カルマンフィルタに基づく従来の雑音抑圧方法のみならずカルマンフィルタ以外の他の原理に基づく従来の雑音抑圧方法に対しても格別の作用効果を有することを実証するために、従来手法１〜３と発明手法１〜２とを比較した。以下、順に説明する。【0213】
（１）客観的評価（波形による評価）図３３は、各手法の音声波形シミュレーションの第１の例を説明するための図、図３４は、各手法の音声波形シミュレーションの第２の例を説明するための図、図３５は、各手法の音声波形シミュレーションの第３の例を説明するための図、図３６は、各手法の音声波形シミュレーションの第４の例を説明するための図である。
【0214】
シミュレーション条件は、次の通りである。【0215】
本シミュレーションでは、２つの音声信号と２つの雑音を用いる。一方の音声信号は、クリアな信号（音声信号）として、成人男性の音声を無音室で録音したものであり、音声（Ａ−１）と称する。もう一方の音声信号は、クリアな信号（音声信号）として、成人女性の音声を無音室で録音したものであり、音声（Ａ−２）と称する。また、一方の雑音は、ガウス性白色雑音（つまり、白色ガウス雑音）であり、雑音（Ｂ−１）と称する。もう一方の雑音は、バブル雑音（有色雑音）であり、雑音（Ｂ−２）と称する。２つの雑音は、いずれも、分散σｖ２は既知である、つまり、別の任意の方法（アルゴリズム）で求められて与えられるとする。【0216】
また、信号雑音比ＳＮＲｉｎを、次の式（５７）で定義する。【0217】
【数57】
【0218】
図３３は、音声（Ａ−１）と雑音（Ｂ−１）の組み合わせにおける音声波形シミュレーションの結果を示す図であり、例えば、ＳＮＲｉｎ＝１０ｄＢ、Ｌｃ＝Ｌｐ＝５０（Ｌは状態遷移行列またはベクトルのサイズ）の条件における各種波形を示している。具体的には、図３３Ａは、音声（Ａ−１）の波形を示し、図３３Ｂは、雑音（Ｂ−１）の波形を示し、図３３Ｃは、音声（Ａ−１）と雑音（Ｂ−１）の合成波形（観測信号）を示している。また、図３３Ｄは、音声（Ａ−１）と雑音（Ｂ−１）の合成波形に対して従来手法１による雑音抑圧を行った場合の推定音声信号の波形を示し、図３３Ｅは、音声（Ａ−１）と雑音（Ｂ−１）の合成波形に対して従来手法２による雑音抑圧を行った場合の推定音声信号の波形を示し、図３３Ｆは、音声（Ａ−１）と雑音（Ｂ−１）の合成波形に対して従来手法３による雑音抑圧を行った場合の推定音声信号の波形を示し、図３３Ｇは、音声（Ａ−１）と雑音（Ｂ−１）の合成波形に対して発明手法１による雑音抑圧を行った場合の推定音声信号の波形を示し、図３３Ｈは、音声（Ａ−１）と雑音（Ｂ−１）の合成波形に対して発明手法２による雑音抑圧を行った場合の推定音声信号の波形を示している。【0219】
図３４は、音声（Ａ−２）と雑音（Ｂ−１）の組み合わせにおける音声波形シミュレーションの結果を示す図であり、例えば、上記と同様に、ＳＮＲｉｎ＝１０ｄＢ、Ｌｃ＝Ｌｐ＝５０の条件における各種波形を示している。具体的には、図３４Ａは、音声（Ａ−２）の波形を示し、図３４Ｂは、雑音（Ｂ−１）の波形を示し、図３４Ｃは、音声（Ａ−２）と雑音（Ｂ−１）の合成波形（観測信号）を示している。また、図３４Ｄは、音声（Ａ−２）と雑音（Ｂ−１）の合成波形に対して従来手法１による雑音抑圧を行った場合の推定音声信号の波形を示し、図３４Ｅは、音声（Ａ−２）と雑音（Ｂ−１）の合成波形に対して従来手法２による雑音抑圧を行った場合の推定音声信号の波形を示し、図３４Ｆは、音声（Ａ−２）と雑音（Ｂ−１）の合成波形に対して従来手法３による雑音抑圧を行った場合の推定音声信号の波形を示し、図３４Ｇは、音声（Ａ−２）と雑音（Ｂ−１）の合成波形に対して発明手法１による雑音抑圧を行った場合の推定音声信号の波形を示し、図３４Ｈは、音声（Ａ−２）と雑音（Ｂ−１）の合成波形に対して発明手法２による雑音抑圧を行った場合の推定音声信号の波形を示している。
【0220】
図３５は、音声（Ａ−１）と雑音（Ｂ−２）の組み合わせにおける音声波形シミュレーションの結果を示す図であり、例えば、上記と同様に、ＳＮＲｉｎ＝１０ｄＢ、Ｌｃ＝Ｌｐ＝５０の条件における各種波形を示している。具体的には、図３５Ａは、音声（Ａ−１）の波形を示し、図３５Ｂは、雑音（Ｂ−２）の波形を示し、図３５Ｃは、音声（Ａ−１）と雑音（Ｂ−２）の合成波形（観測信号）を示している。また、図３５Ｄは、音声（Ａ−１）と雑音（Ｂ−２）の合成波形に対して従来手法１による雑音抑圧を行った場合の推定音声信号の波形を示し、図３５Ｅは、音声（Ａ−１）と雑音（Ｂ−２）の合成波形に対して従来手法２による雑音抑圧を行った場合の推定音声信号の波形を示し、図３５Ｆは、音声（Ａ−１）と雑音（Ｂ−２）の合成波形に対して従来手法３による雑音抑圧を行った場合の推定音声信号の波形を示し、図３５Ｇは、音声（Ａ−１）と雑音（Ｂ−２）の合成波形に対して発明手法１による雑音抑圧を行った場合の推定音声信号の波形を示し、図３５Ｈは、音声（Ａ−１）と雑音（Ｂ−２）の合成波形に対して発明手法２による雑音抑圧を行った場合の推定音声信号の波形を示している。【0221】
図３６は、音声（Ａ−２）と雑音（Ｂ−２）の組み合わせにおける音声波形シミュレーションの結果を示す図であり、例えば、上記と同様に、ＳＮＲｉｎ＝１０ｄＢ、Ｌｃ＝Ｌｐ＝５０の条件における各種波形を示している。具体的には、図３６Ａは、音声（Ａ−２）の波形を示し、図３６Ｂは、雑音（Ｂ−２）の波形を示し、図３６Ｃは、音声（Ａ−２）と雑音（Ｂ−２）の合成波形（観測信号）を示している。また、図３６Ｄは、音声（Ａ−２）と雑音（Ｂ−２）の合成波形に対して従来手法１による雑音抑圧を行った場合の推定音声信号の波形を示し、図３６Ｅは、音声（Ａ−２）と雑音（Ｂ−２）の合成波形に対して従来手法２による雑音抑圧を行った場合の推定音声信号の波形を示し、図３６Ｆは、音声（Ａ−２）と雑音（Ｂ−２）の合成波形に対して従来手法３による雑音抑圧を行った場合の推定音声信号の波形を示し、図３６Ｇは、音声（Ａ−２）と雑音（Ｂ−２）の合成波形に対して発明手法１による雑音抑圧を行った場合の推定音声信号の波形を示し、図３６Ｈは、音声（Ａ−２）と雑音（Ｂ−２）の合成波形に対して発明手法２による雑音抑圧を行った場合の推定音声信号の波形を示している。【0222】
まず、図３３Ａと図３３Ｄ、および、図３４Ａと図３４Ｄをそれぞれ比較することにより、従来手法１による雑音抑圧では、雑音抑圧後に推定音声信号の振幅が小さくなっており、クリアな信号（音声信号）が抑圧されていることがわかる。また、従来手法１による雑音抑圧では、サンプル数の増加と共に、雑音抑圧後の推定音声信号の波形がクリアな信号（音声信号）の波形から変形していることがわかる。【0223】
また、特に図３４Ａと図３４Ｄの比較から、従来手法１の雑音抑圧では、無声区間を有する音声（Ａ−２）に対して、推定音声信号が抑圧されるだけでなく、無声区間においてオリジナルの雑音とは異なる雑音が観察されている。これは、従来手法１では、無声区間において音声信号ｄ（ｎ）は０であるにもかかわらず、上記の式（２）でＡＲ係数を求めようとするためにＡＲ係数の値が発散し、不安定な状態を与えるからであると推測される。【0224】
また、このことから、雑音が有色の場合、従来手法１の適用は困難であろうことは容易に推測される。【0225】
これに対し、発明手法１〜２による雑音抑圧では、従来手法１と対照的に、雑音抑圧後の推定音声信号の波形は、いずれの場合においても、クリアな信号（音声信号）の波形と非常に似ている。【0226】
次に、図３５Ａと図３５Ｄ、および、図３６Ａと図３６Ｄをそれぞれ比較することにより、従来手法１による雑音抑圧では、雑音（Ｂ−２）を含む観測音声信号に対して、非常に劣った結果を与えていることがわかる。これは、従来手法１では、有色雑音である雑音（Ｂ−２）を含んだ観測音声信号に対してＡＲ係数を正確に推定することが困難であるためである。【0227】
これに対し、発明手法１〜２による雑音抑圧法では、従来手法１と対照的に、雑音（Ｂ−２）の場合も、雑音（Ｂ−１）の場合と同程度の雑音抑圧が達成されている。【0228】
このように、本発明の雑音抑制方法（発明手法１〜２）は、白色雑音か有色雑音かにかかわらず、また、無声区間の有無にかかわらず、有効である。これは、本発明の雑音抑制方法の大きな特徴の１つである。【0229】
（２）客観的評価（数値による評価）図３７は、各手法の雑音抑圧能力の数値シミュレーションの第１の例の結果を示す図、図３８は、各手法の雑音抑圧能力の数値シミュレーションの第２の例の結果を示す図、図３９は、各手法の雑音抑圧能力の数値シミュレーションの第３の例の結果を示す図、図４０は、各手法の雑音抑圧能力の数値シミュレーションの第４の例の結果を示す図である。
【0230】
本シミュレーションでは、雑音抑制能力を数値により評価するため、次の式（５８）で表されるＳＮＲｏｕｔ［ｄＢ］を用いて雑音抑圧量を評価した。なお、ＳＮＲは、信号に対する雑音の比であり、数値が大きいほど雑音が少なく音声として良好であるといえる。また、図３７〜図４０において、「Ｌ」は、上記の「Ｌｃ」と「Ｌｐ」をまとめて表記したものである。【0231】
【数58】
【0232】
すなわち、図３７は、音声（Ａ−１）と雑音（Ｂ−１）の組み合わせにおける雑音抑圧能力の数値シミュレーションの結果を示す図であり、成人男性の音声、白色雑音の条件における雑音抑圧量ＳＮＲｏｕｔを示している。図３８は、音声（Ａ−２）と雑音（Ｂ−１）の組み合わせにおける雑音抑圧能力の数値シミュレーションの結果を示す図であり、成人女性の音声、白色雑音の条件における雑音抑圧量ＳＮＲｏｕｔを示している。図３９は、音声（Ａ−１）と雑音（Ｂ−２）の組み合わせにおける雑音抑圧能力の数値シミュレーションの結果を示す図であり、成人男性の音声、有色雑音の条件における雑音抑圧量ＳＮＲｏｕｔを示している。図４０は、音声（Ａ−２）と雑音（Ｂ−２）の組み合わせにおける雑音抑圧能力の数値シミュレーションの結果を示す図であり、成人女性の音声、有色雑音の条件における雑音抑圧量ＳＮＲｏｕｔを示している。また、４つのいずれの場合においても、ＳＮＲｉｎと状態遷移行列またはベクトルのサイズＬとのいくつかの値の組み合わせに対して、従来手法１〜３と発明手法１〜２によるＳＮＲｏｕｔの値を比較して示している。なお、発明手法１と発明手法２の結果はほぼ同じであるため、ここでは、便宜上、両者をまとめて「発明手法１＆２」と表記している。【0233】
図３７〜図４０を参照すると、発明手法１〜２は、すべてのＳＮＲｉｎ、Ｌｐの値において、従来手法１〜３に比べて雑音抑圧能力を改善していることがわかる。【0234】
特に、図３９および図４０に示す有色雑音の場合には、従来手法１は非常に劣った結果を与えているのに対し、発明手法１〜２は、図３７および図３８に示す白色雑音の場合と同程度の結果を示している。すなわち、発明手法１〜２は、白色雑音と有色雑音の両者に効果的で、雑音の性質に堅牢な雑音抑圧方法であるといえる。【0235】
また、図３７〜図４０に見られるように、発明手法１〜２では、Ｌｐの値に対して雑音抑圧能力ＳＮＲｏｕｔは安定であり、Ｌｐの値の増加に伴い増加する傾向にある。これと対照的に、従来手法１では、図３７〜図４０に見られるように、Ｌｃの値に対して雑音抑圧能力ＳＮＲｏｕｔは不安定である。これは、従来手法１では、最適なＬｃの値、つまりＡＲ係数の次数を決定することが困難であることを意味している。【0236】
ＡＲ係数の推定を必要とする従来手法１において最も問題になることは、一般に、ＡＲ係数の次数の決定は非常に困難であるということである。なぜなら、ＡＲ係数の次数の正確な推定は、例えば、雑音抑圧であれば、クリアな信号（音声信号）に依存しているからである。【0237】
このことは、クリアな信号（音声信号）が既知でなければならないことを意味しているため、リアルタイム処理は困難となる。ＡＲ係数の次数が適切でない場合には、カルマンフィルタアルゴリズムの性能が劣化することは容易に想像可能である。また、何らかの手法でリアルタイムに推定することが可能となったとしても、処理が増加することより演算量などの問題を避けることは不可能である。【0238】
さらに、図３７〜図４０に見られるように、従来手法３と発明手法１〜２とを比較した場合、発明手法１および発明手法２は、４つのいずれの場合においても、従来手法３に比べて雑音抑圧能力が高いといえる。特に、ＳＮＲｉｎ、Ｌｐの値の値が小さいほど、発明手法１および発明手法２は、従来手法３よりも高い雑音抑圧能力を発揮することがわかる。このような雑音抑圧能力の違いは、演算量の違いに起因するものと考えられる。すなわち、演算量が少ないほど、誤差の発生が少なくなり、雑音抑圧能力の向上に資するものと考えられる。【0239】
また、図３７〜図４０に見られるように、従来手法２と発明手法１〜２とを比較した場合、発明手法１および発明手法２は、４つのいずれの場合においても、従来手法３との比較の場合にも増して、従来手法２に比べて雑音抑圧能力が高いといえる。従来手法２は、音声に特化したパラメータの設定を必要とするが、この実験の結果は、実環境では、そのパラメータの設定が非常に困難であることを示している。【0240】
（３）主観的評価（聞き取り調査）図４１は、各手法の主観的評価の第１の例の結果を示す図、図４２は、各手法の主観的評価の第２の例の結果を示す図、図４３は、各手法の主観的評価の第３の例の結果を示す図、図４４は、各手法の主観的評価の第４の例の結果を示す図である。
【0241】
ここでは、推定音声信号の音声品質を評価するために、リスニングテスト（聞き取り調査）による主観的評価を行った。音声品質評価に用いた音声信号と雑音は、上記のシミュレーション（客観的評価）に用いたものと同一である。雑音は、異なるＳＮＲｉｎ（＝０、５［ｄＢ］）で音声信号に加えた。また、音声品質評価は、ＡＣＲ（絶対範疇評価）に基づく５段階ＭＯＳ（平均オピニオン値）を用いたリスニングテストにより行った。５０人の聴取者が雑音抑圧により得られた推定音声信号のうちいくつかを評価した。各々の聴取者は、ポイント１からポイント５を決定する。ポイント５が最良である。【0242】
図４１は、音声（Ａ−１）と雑音（Ｂ−１）の組み合わせにおける雑音抑圧後の音声品質の主観的評価結果の一例を示す図であり、例えば、Ｌｃ＝Ｌｐ＝５０の条件における各手法のリスニングテストの結果を示している。図４２は、音声（Ａ−２）と雑音（Ｂ−１）の組み合わせにおける雑音抑圧後の音声品質の主観的評価結果の一例を示す図であり、例えば、Ｌｃ＝Ｌｐ＝５０の条件における各手法のリスニングテストの結果を示している。図４３は、音声（Ａ−１）と雑音（Ｂ−２）の組み合わせにおける雑音抑圧後の音声品質の主観的評価結果の一例を示す図であり、例えば、Ｌｃ＝Ｌｐ＝５０の条件における各手法のリスニングテストの結果を示している。図４４は、音声（Ａ−２）と雑音（Ｂ−２）の組み合わせにおける雑音抑圧後の音声品質の主観的評価結果の一例を示す図であり、例えば、Ｌｃ＝Ｌｐ＝５０の条件における各手法のリスニングテストの結果を示している。なお、客観的評価のみならず主観的評価においても発明手法１と発明手法２の評価結果はほぼ同じであるため、ここでも、便宜上、両者をまとめて「発明手法１＆２」と表記している。【0243】
図４１〜図４４から、発明手法１〜２で推定した音声信号のスコアは、すべてのＳＮＲｉｎの値において、従来手法１〜３のスコアよりも高いことがわかる。特に、従来手法１〜２との比較において、発明手法１〜２と従来手法１〜２との差は、音声（Ａ−２）と雑音（Ｂ−１）の組み合わせ、音声（Ａ−１）と雑音（Ｂ−２）の組み合わせ、および音声（Ａ−２）と雑音（Ｂ−２）の組み合わせに対して大きい。すなわち、発明手法１〜２では、従来手法１〜２に比べて音質が大幅に改善されている。【0244】
また、従来手法３との比較においても、発明手法１および発明手法２は、４つのいずれの場合においても、従来手法３に比べて主観的評価も高いといえる。【0245】
以上の実験結果により、本発明の雑音抑圧方法（発明手法１〜２）は、音声信号の音声品質を犠牲にすることのない、白色雑音および有色雑音に効果的な優れた雑音抑圧方法であるといえる。特に、従来手法３との比較において、発明手法１および発明手法２は、従来手法３に比べて、数値による客観的評価が高く、より高い雑音抑圧能力を実現しているといえる。【0246】
以上をまとめると、本発明に係る雑音抑圧方法（発明手法１〜２）は、新しい状態空間モデル（特に観測方程式）を構成することによって、従来手法１〜３に比べて、演算量を大幅に軽減することが可能である。具体的には、発明手法１〜２は、まず、ＡＲ係数の推定を必要としないため、従来手法１で必要であったＡＲ係数を推定するステップを削減することができ、従来手法１に比べて演算量を大幅に軽減することができる（図３０、図３１参照）。また、上記のように、発明手法１〜２は、従来手法１と異なり、白色雑音か有色雑音かにかかわらず、また、無声区間の有無にかかわらず、有効であり、従来手法１に比べて、客観的にも主観的にも高い雑音抑圧能力を実現することができる（図３３〜図４４参照）。また、発明手法１〜２は、従来手法１で必要であったＡＲ係数を推定するステップを必要としないため、ハードウエアとソフトウエアのいずれかで構成するとしても、従来手法１よりも単純化されたシンプルな構成で実現することができる。したがって、発明手法１〜２によれば、よりシンプルな構成で、かつ、より少ない演算量で、所望情報の品質を劣化させることなくより高い雑音抑圧能力を実現することができる。さらに、演算量の大幅な軽減によって、メモリ容量も大幅に削減することができる。また、音声の場合には、音質も大幅に改善することができる。【0247】
特に、従来手法３との比較においても、発明手法１および発明手法２は、上記のように、従来手法３に比べて、さらに演算量の軽減（図３０、図３２参照）および雑音抑圧能力の向上（特に図３７〜図４０参照）を図ることができ、従来技術３に対しても顕著な効果を有する。さらに、演算量の軽減によってメモリ容量も削減することができ、また、音声の場合には音質も改善することができる。【0248】
さらに、発明手法２は、上記のように、発明手法１に比べて、さらなるスカラー化が図られ、反復過程の手順の数が５ステップから３ステップに低減されているため、さらに演算量を軽減することができ、また、状態遷移行列Φおよび観測遷移ベクトルｍの設定ならびに多くの計算結果の保存が不必要になるため、より一層メモリ容量の削減を図ることができる。【0249】
また、発明手法１〜２を、例えば、半導体集積回路や半導体回路などのハードウエアとして実施する場合や、パーソナルコンピュータなどで実行可能なソフトウエアとして実施する場合のいずれにおいても、その構成は、従来手法よりも単純化される。したがって、発明手法１〜２を用いれば、回路規模やプログラム量を大幅に低減できるであろうことは明らかである。【0250】
なお、本発明に係る雑音抑圧装置および雑音抑圧方法は、いろいろな技術分野に適用可能である。【0251】
例えば、本発明の雑音抑圧装置は、雑音を含む音声信号（観測信号）からクリアな信号（所望信号）としての音声信号を取得することが可能である。この分野の適用例としては、例えば、携帯電話やカーナビゲーション、対話型ロボットなどにおける音声情報処理を挙げることができる。具体的には、一例として、例えば、カーナビゲーションの場合は、カーナビゲーションに必要不可欠な音声認識装置の前処理雑音抑圧装置への適用が考えられる。【0252】
また、画像処理の分野においては、本発明によれば、何らかの原因でぼけと雑音が入ってしまった劣化画像（観測信号）からぼけと雑音を取り除いた原画像をクリアな信号（所望信号）として取得することが可能であり、画像処理装置として活用可能である。【0253】
さらに、従来手法が適用されていた通信や信号処理の分野全般にわたって、本発明が適応可能であることはいうまでもない。【0254】
また、医療分野では、従来、妊婦の胎児の状況を検査するには、個人が購入できない高価な装置と高い専門知識とが必要であったが、本発明によれば、妊婦の体から得られる観測信号（母胎の心拍その他の雑音を含む）から不必要な音（雑音）を抑圧し、胎児の心拍（所望信号）を取得することが可能になり、通院せずとも自宅で胎児の健康状態を、その心拍から容易に確認することが可能となる。また、出産前の胎児のみならず出産後の赤ちゃんの管理にも役立つことが考えられる（周産期医療）。なお、ここで「心拍」とは、広く心臓の動きを意味し、例えば、心音や心電図なども含まれる。【0255】
また、本実施の形態の説明に用いた各機能要素は、例えば、集積回路として実現される。これらは、個別に１チップ化されてもよいし、一部またはすべてを含むように１チップ化されてもよい。また、集積回路製造後にプログラムすることが可能なＦＰＧＡ（Ｆield Programmable Gate Array）や、回路を構成可能なリコンフィギュラブル・プロセッサを利用してもよい。【0256】
さらに、本実施の形態は、ハードウエアに限定されるものではなく、ソフトウエアによってもよい。その逆も真である。また、それらの組み合わせであってもよい。【0257】
上記のように、本発明に係る雑音抑圧装置および雑音抑圧方法は、いろいろな技術分野に適用可能である。そこで、以下では、本発明に係る雑音抑圧装置および雑音抑圧方法の各種分野への具体的な適用例について説明する。ここでは、本発明に係る雑音抑圧装置および雑音抑圧方法を、例えば、音響機器関連（胎児心拍検出装置や携帯電話など）、音声認識装置関連（カーナビゲーションなど）、検査装置関連（異常音検出装置など）、および画像処理装置関連（画像復元装置など）にそれぞれ適用した場合について、具体的に説明する。なお、各適用例において、雑音抑圧処理部は、上記の発明手法１〜２のいずれの手法をも任意に実行することができる。【0258】
（実施の形態２）実施の形態２は、実施の形態１に係る雑音抑圧装置を胎児心拍検出装置に適用した場合である。なお、ここで「心拍」とは、上記のように、広く心臓の動きを意味し、例えば、心音や心電図なども含まれる。
【0259】
図４５は、本発明の実施の形態２に係る胎児心拍検出装置の構成を示すブロック図である。【0260】
図４５に示す胎児心拍検出装置４００は、実施の形態１の雑音抑圧処理（発明手法１〜２）を実行可能なコンピュータ本体４１０と、マイクロホン４２０と、信号入力部４３０と、操作部４４０と、スピーカ４５０と、ディスプレイ４６０とを有する。【0261】
コンピュータ本体４１０は、インタフェース部４１１、記憶部４１２（記録装置４１３と主記憶メモリ４１４からなる）、モデム４１５、Ｄ／Ａ変換器４１６、雑音抑圧処理部４１７、胎児心拍解析処理部４１８、および正常胎児心拍情報記憶部４１９を有する。雑音抑圧処理部４１７および胎児心拍解析処理部４１８は、中央演算処理装置（ＣＰＵ：Central Processing Unit）で構成されている。コンピュータ本体４１０は、モデム４１５を介して外部の通信ネットワーク（例えば、電話回線やＬＡＮ、インターネットなど）と接続されている。また、コンピュータ本体４１０は、Ｄ／Ａ変換器４１６を介してスピーカ４５０と接続されている。スピーカ４５０は、例えば、聴診器のスピーカであってもよい。なお、図示しないが、コンピュータ本体４１０にプリンタを接続することも可能である。【0262】
コンピュータ本体４１０において、実施の形態１の雑音抑圧処理を実行させるプログラム（雑音抑圧アルゴリズム）は、記録装置４１３に格納されていてもよいし、モデム４１５およびインタフェース部４１１を介して外部からダウンロードされてきてもよい。記録装置４１３は、典型的にはハードディスク装置であるが、ＣＤ−ＲＯＭ装置やＤＶＤ装置、フラッシュメモリなどの可搬性のあるものであってもよいし、また、それらの組み合わせであってもよい。雑音抑圧処理部４１７は、そのプログラム（雑音抑圧アルゴリズム）を実行することによって、マイクロホン４２０により検出された胎児の心拍（雑音を含む）に対して実施の形態１の雑音抑圧処理を実行して胎児の心拍を取得する。【0263】
信号入力部４３０は、図１に示すサンプリング部１２０およびＡ／Ｄ変換部１３０を有する。信号入力部４３０は、コンピュータ本体４１０の内部に格納された内蔵カード（ボード）であってもよいし、インタフェース部４１１を経由して接続された外部設置型機器であってもよい。【0264】
操作部４４０は、典型的にはキーボートやマウス、タッチパネルなどであるが、音声認識装置などを用いてもよい。使用者は、操作部４４０を用い、ディスプレイ４６０で確認をしながらコンピュータを操作することができる。また、操作部４４０は、パラメータ設定部４４１を有する。パラメータ設定部４４１は、使用者の入力操作により、実施の形態１の雑音抑圧処理に必要な各種パラメータの値を設定し、コンピュータ本体４１０に出力する。【0265】
胎児心拍解析処理部４１８は、雑音抑圧処理部４１７によって取得された胎児の心拍を解析する。例えば、胎児心拍解析処理部４１８は、与えられた正常胎児心拍情報（胎児の正常心拍）と、雑音抑圧処理部４１７によって取得された胎児の心拍とを比較して、異常心拍を識別したり症状の診断を行ったりする。この場合、雑音抑圧処理部４１７は、胎児心拍解析処理部４１８の前処理としての機能を有する。ここでは、正常胎児心拍情報は、正常胎児心拍情報記憶部４１９に格納されている。正常胎児心拍情報記憶部４１９は、例えば、検診対象となる胎児の胎児月齢を操作部４４０から入力することにより、胎児月齢ごとに正常胎児心拍情報を収集し蓄積する。正常胎児心拍情報は、この正常胎児心拍情報記憶部４１９から読み出されて胎児心拍解析処理部４１８に与えられる。【0266】
マイクロホン４２０からの観測音声信号は、信号入力部４３０のサンプリング部１２０に入力される。サンプリング部１２０は、所定のサンプリング周波数（例えば、１６ｋＨｚ）で、入力されたアナログの観測音声信号をサンプリング処理し、Ａ／Ｄ変換部１３０に出力する。Ａ／Ｄ変換部１３０は、サンプリングされた観測音声信号の振幅値を所定の分解能（例えば、８ｂｉｔ）でＡ／Ｄ変換処理し、一時格納する。Ａ／Ｄ変換部１３０は、所定のサンプリング数Ｎの音声フレーム単位で、デジタル化した観測音声信号をコンピュータ本体４１０のインタフェース部４１１に出力する。【0267】
コンピュータ本体４１０は、インタフェース部４１１に出力された観測音声信号を一時、記憶部４１２の主記憶メモリ４１４に格納し、その後、所定の音声フレーム（サンプリング数）単位で、雑音抑圧処理を施した上で、再度主記憶メモリ４１４に格納する。雑音抑圧処理は、主記憶メモリ４１４や記録装置４１３に格納されたソフトウエアをインタフェース部４１１経由で雑音抑圧処理部４１７に呼び出し、実行させることによって行われる。【0268】
コンピュータ本体４１０は、使用者の操作により、処理を実行したり、中断、終了させたりする。また、コンピュータ本体４１０は、使用者の操作により、雑音抑圧処理部４１７で取得した推定音声信号（胎児の心拍）を、胎児心拍解析処理部４１８に出力したり、モデム４１５やスピーカ４５０、ディスプレイ４６０などを介して外部に出力してもよい。【0269】
このように構成された胎児心拍検出装置４００は、例えば、検出した胎児の心拍を聴診器のスピーカに出力したり、あるいは、正常心拍との比較で解析した結果をモデム４１５を介して特定の医療センターに送信したりすることができる。送信を受けた医療センターは、受信した解析結果に基づいて総合判定を行うことができる。また、解析結果については、単独の情報としてディスプレイ４６０に表示したり、あるいは、前回検診時の心拍データを記憶部４１２から読み出して両者を対比する形でディスプレイ４６０に表示させたりすることもできる。また、解析結果を図示しないプリンタに出力して視認することも可能である。【0270】
このように、本実施の形態によれば、妊婦の体から得られる観測音声信号（母胎の心拍その他の雑音を含む）から不必要な音（雑音）を抑圧し、胎児の心拍（所望信号）を取得することができるため、胎児の健康状態を、その心拍から容易にかつ正確に確認することができる。また、出産前の胎児のみならず出産後の赤ちゃんの管理にも役立つことが考えられる（周産期医療・生体観測）。【0271】
なお、本実施の形態では、胎児心拍検出装置４００は、出力手段として、スピーカ、ディスプレイ、通信手段、記憶装置、およびプリンタ（図示せず）を有するが、これに限定されない。出力手段については、用途や機能などに応じて、適宜、取捨選択すればよい。また、胎児の心拍さえ検出できればよい場合は、胎児心拍解析処理部４１８および正常胎児心拍情報記憶部４１９は省略してもよい。【0272】
また、本実施の形態では、コンピュータ本体４１０に胎児心拍解析処理部４１８を設けているが（内蔵タイプ）、もちろんこれに限定されない。胎児心拍解析処理部４１８は、コンピュータ本体４１０の外部装置（胎児心拍解析装置）として構成することも可能である（外付けタイプ）。胎児心拍解析処理部４１８をコンピュータ本体４１０に組み込むか外部装置（胎児心拍解析装置）として構成するかは、用途やデータ処理量などに応じて任意に決定すればよい。また、この点は、正常胎児心拍情報記憶部４１９についても同様である。【0273】
（実施の形態３）実施の形態３は、実施の形態１に係る雑音抑圧装置を携帯電話などの携帯端末装置に適用した場合である。
【0274】
図４６は、本発明の実施の形態３に係る携帯端末装置の構成を示すブロック図である。なお、この携帯端末装置５００は、図４５に示す胎児心拍検出装置４００と同様の基本的構成を有しており、同一の構成要素には同一の符号を付し、オプションを含めてその説明を省略する。【0275】
図４６に示す携帯端末装置５００は、送受信用のアンテナ５１０と、送受信部５２０とを有する。送受信部５２０は、アンテナ５１０で送受信した音声信号をベースバンド処理する。【0276】
この携帯端末装置５００では、マイクロホン４２０からの観測音声信号（使用者の音声信号）は、信号入力部４３０に入力され、ここでサンプリング部１２０およびＡ／Ｄ変換部１３０によりデジタル化された後、インタフェース部４１１に出力される。インタフェース部４１１に出力された観測音声信号は、記憶部４１２に一時格納され、その後、所定の音声フレーム（サンプリング数）単位で、雑音抑圧処理が施された後、記憶部４１２に再度格納される。雑音抑圧処理は、記憶部４１２に格納されたプログラム（雑音抑圧アルゴリズム）をインタフェース部４１１経由で雑音抑圧処理部４１７に呼び出し、実行させることによって行われる。雑音抑圧処理後のクリアな信号（音声信号）は、送受信部５２０で送信用のベースバンド処理が施された後、アンテナ５１０から無線送信される。【0277】
一方、アンテナ５１０で受信した観測音声信号（通信相手の音声信号）は、送受信部５２０で受信用のベースバンド処理が施された後、デジタル信号としてインタフェース部４１１に出力される。インタフェース部４１１に出力された観測音声信号は、記憶部４１２に一時格納され、その後、所定の音声フレーム（サンプリング数）単位で、雑音抑圧処理が施された後、記憶部４１２に再度格納される。雑音抑圧処理後のクリアな信号（音声信号）は、Ｄ／Ａ変換器４１６を介してスピーカ４５０に出力される。【0278】
このように、本実施の形態によれば、マイクロホン４２０からの観測音声信号（雑音を含む）およびアンテナ５１０で受信した観測音声信号（雑音を含む）からそれぞれ不必要な音（雑音）を抑圧し、音声信号をクリアな信号（所望信号）として取得することができるため、簡単な構成で高音質化を図ることができる。【0279】
（実施の形態４）実施の形態４は、実施の形態１に係る雑音抑圧装置をカーナビゲーション装置に適用した場合である。カーナビゲーション装置として、ドライバが運転に集中しながら情報を入力できるよう、つまり、音声で操作できるよう、音声認識機能を搭載したものがある。本実施の形態は、このような音声認識機能を有するカーナビゲーション装置において、当該カーナビゲーションに必要不可欠な音声認識装置の前処理としての雑音抑圧装置への適用例である。
【0280】
図４７は、本発明の実施の形態４に係るカーナビゲーション装置の構成を示すブロック図である。なお、このカーナビゲーション装置６００は、図４５に示す胎児心拍検出装置４００と同様の基本的構成を有しており、同一の構成要素には同一の符号を付し、オプションを含めてその説明を省略する。【0281】
図４７に示すカーナビゲーション装置６００は、音声認識処理部６１０、行先設定部６２０、現在位置検出部６３０、およびルート演算部６４０を有する。この構成により、カーナビゲーション装置６００は、音声認識処理部６１０で、検出した音声（使用者の指示）を認識し、行先設定部６２０で、その認識した指示に基づいて行先を設定し、現在位置検出部６３０で、ＧＰＳ（Global Positioning System）からの信号を受信して現在位置を把握し、ルート演算部６４０で、音声により指示された行先から最適なルートを演算し、ディスプレイ４６０に、演算されたルートや地図などを表示する。なお、このカーナビゲーション処理は、単なる一例であって、音声操作できる処理であれば任意のカーナビゲーション処理に適用可能である。【0282】
このとき、このカーナビゲーション装置６００では、マイクロホン４２０からの観測音声信号（使用者の指示）は、信号入力部４３０に入力され、ここでサンプリング部１２０およびＡ／Ｄ変換部１３０によりデジタル化された後、インタフェース部４１１に出力される。インタフェース部４１１に出力された観測音声信号は、記憶部４１２の主記憶メモリ４１４に一時格納され、その後、所定の音声フレーム（サンプリング数）単位で、雑音抑圧処理が施された後、主記憶メモリ４１４に再度格納される。雑音抑圧処理は、記憶部４１２（記録装置４１３または主記憶メモリ４１４）に格納されたプログラム（雑音抑圧アルゴリズム）をインタフェース部４１１経由で雑音抑圧処理部４１７に呼び出し、実行させることによって行われる。雑音抑圧処理後のクリアな信号（音声信号）は、音声認識処理部６１０に出力される。【0283】
このように、本実施の形態によれば、音声認識処理部６１０の前処理として、マイクロホン４２０からの観測音声信号（雑音を含む）から不必要な音（雑音）を抑圧し、音声信号をクリアな信号（所望信号）として取得することができるため、音声認識処理部６１０の音声認識能力を最大限に発揮させることができ、カーナビゲーションを音声認識で確実に操作することができる。特に、運転中は音声指示以外の雑音がかなり多くかつ大きいため、この効果は非常に顕著である。【0284】
（実施の形態５）実施の形態５は、実施の形態１に係る雑音抑圧装置を音声認識装置に適用した場合である。
【0285】
図４８は、本発明の実施の形態５に係る音声認識装置の構成を示すブロック図である。なお、この音声認識装置７００は、図４７に示すカーナビゲーション装置６００と同様の基本的構成を有しており、同一の構成要素には同一の符号を付し、オプションを含めてその説明を省略する。【0286】
図４８に示す音声認識装置７００は、図４７に示すカーナビゲーション装置６００から、カーナビゲーションに特有の構成要素、つまり、行先設定部６２０、現在位置検出部６３０、およびルート演算部６４０を削除した構成を有する。この音声認識装置７００においても、雑音抑圧処理部４１７は、図４７に示すカーナビゲーション装置６００の場合と同様に、音声認識処理部６１０の前処理としての機能を有する。【0287】
すなわち、この音声認識装置７００では、マイクロホン４２０からの観測音声信号は、信号入力部４３０に入力され、ここでサンプリング部１２０およびＡ／Ｄ変換部１３０によりデジタル化された後、インタフェース部４１１に出力される。インタフェース部４１１に出力された観測音声信号は、記憶部４１２の主記憶メモリ４１４に一時格納され、その後、所定の音声フレーム（サンプリング数）単位で、雑音抑圧処理が施された後、主記憶メモリ４１４に再度格納される。雑音抑圧処理は、記憶部４１２（記録装置４１３または主記憶メモリ４１４）に格納されたプログラム（雑音抑圧アルゴリズム）をインタフェース部４１１経由で雑音抑圧処理部４１７に呼び出し、実行させることによって行われる。雑音抑圧処理後のクリアな信号（音声信号）は、音声認識処理部６１０に出力される。【0288】
このように、本実施の形態によれば、音声認識処理部６１０の前処理として、マイクロホン４２０からの観測音声信号（雑音を含む）から不必要な音（雑音）を抑圧し、音声信号をクリアな信号（所望信号）として取得することができるため、音声認識処理部６１０の音声認識能力を最大限に発揮させることができ、非常に高精度の音声認識を実現することができる。【0289】
（実施の形態６）実施の形態６は、実施の形態１に係る雑音抑圧装置を異常検出装置に適用した場合である。
【0290】
図４９は、本発明の実施の形態６に係る異常検出装置の構成を示すブロック図である。なお、この異常検出装置８００は、図４５に示す胎児心拍検出装置４００と同様の基本的構成を有しており、同一の構成要素には同一の符号を付し、オプションを含めてその説明を省略する。【0291】
図４９に示す異常検出装置８００は、検査対象が発する異常音を検出する装置であって、異常音解析処理部８１０および正常音情報記憶部８２０を有する。異常音解析処理部８１０は、雑音抑圧処理部４１７によって取得された音を解析する。例えば、異常音解析処理部８１０は、与えられた正常音情報と、雑音抑圧処理部４１７によって取得された音とを比較して、検査対象の異常音を検出する。この場合、雑音抑圧処理部４１７は、異常音解析処理部８１０の前処理としての機能を有する。ここでは、正常音情報は、正常音情報記憶部８２０に格納されている。本実施の形態では、例えば、検査対象が発生する音を複数の所定箇所から収集するように構成されており、正常音情報記憶部８２０は、検査対象の複数の所定箇所における正常動作音をそれぞれ収集し記憶する。正常音情報は、検査対象の所定箇所ごとに正常音情報記憶部８２０から読み出されて異常音解析処理部８１０に与えられる。また、正常音情報記憶部８２０には、直前に検出された正常動作音をも蓄積するようにしてもよい。この場合、正常音情報記憶部８２０は、正常音情報を、直近の過去分の正常動作音として出力することができる。【0292】
この異常検出装置８００では、マイクロホン４２０からの観測音声信号は、信号入力部４３０に入力され、ここでサンプリング部１２０およびＡ／Ｄ変換部１３０によりデジタル化された後、インタフェース部４１１に出力される。インタフェース部４１１に出力された観測音声信号は、記憶部４１２の主記憶メモリ４１４に一時格納され、その後、所定の音声フレーム（サンプリング数）単位で、雑音抑圧処理が施された後、主記憶メモリ４１４に再度格納される。雑音抑圧処理は、記憶部４１２（記録装置４１３または主記憶メモリ４１４）に格納されたプログラム（雑音抑圧アルゴリズム）をインタフェース部４１１経由で雑音抑圧処理部４１７に呼び出し、実行させることによって行われる。雑音抑圧処理後のクリアな信号（音声信号）は、異常音解析処理部８１０に出力される。【0293】
この異常検出装置８００は、検査対象から検出した音を正常音と比較して解析した結果を、ディスプレイ４６０に表示したり、あるいは、異常音を検出した場合はスピーカ４５０から警報を発したりすることができる。また、解析結果は、モデム４１５を介して特定の監視センターなどに送信することができる。この場合、異常音の検出状況などを、リモートで監視センターなどに連絡することができる。なお、異常音をどのようなタイミングで検出するかは、各装置に依存する。【0294】
このように、本実施の形態によれば、異常音解析処理部８１０の前処理として、マイクロホン４２０からの観測音声信号（雑音を含む）から不必要な音（雑音）を抑圧し、音声信号をクリアな信号（所望信号）として取得することができるため、異常音解析処理部８１０の異常音解析能力を最大限に発揮させることができ、非常に高精度の異常音検出を実現することができる。【0295】
（実施の形態７）実施の形態７は、実施の形態１に係る雑音抑圧装置を画像処理装置、特に画像復元装置に適用した場合である。
【0296】
図５０は、本発明の実施の形態７に係る画像復元装置の構成を示すブロック図である。なお、この画像復元装置９００は、図４５に示す胎児心拍検出装置４００と同様の基本的構成を有しており、同一の構成要素には同一の符号を付し、オプションを含めてその説明を省略する。【0297】
図５０に示す画像復元装置９００は、例えば、ぼけと雑音により劣化した画像を復元する装置であって、画像復元処理部９１０を有する。ここでは、一例として、例えば、スキャナ９２０で読み取った画像を復元する場合を例にとって説明する。すなわち、この画像復元装置９００は、デジタルスキャナ（以下単に「スキャナ」という）９２０で読み取った画像に含まれるぼけと雑音を抑圧してクリーンな画像を出力する構成を有する。この場合、雑音抑圧処理部４１７は、画像復元処理部９１０の前処理としての機能を有する。【0298】
この画像復元装置９００では、スキャナ９２０からの観測画像信号は、信号入力部４３０ａに入力され、ここでサンプリング部１２０によりサンプリング処理された後、インタフェース部４１１に出力される。インタフェース部４１１に出力された観測画像信号は、記憶部４１２の主記憶メモリ４１４に一時格納され、その後、所定の画像フレーム（サンプリング数）単位で、雑音抑圧処理が施された後、主記憶メモリ４１４に再度格納される。雑音抑圧処理は、記憶部４１２（記録装置４１３または主記憶メモリ４１４）に格納されたプログラム（雑音抑圧アルゴリズム）をインタフェース部４１１経由で雑音抑圧処理部４１７に呼び出し、実行させることによって行われる。雑音抑圧処理後のクリーンな画像信号は、画像復元処理部９１０に出力される。なお、画像復元処理部９１０によって復元された画像は、プリンタ９３０またはディスプレイ４６０に出力される。【0299】
このとき、観測画像信号に対する雑音抑圧処理は、スキャナ９２０で読み取った画像の特定領域を指定することにより、指定された領域の画像に対してのみぼけと雑音の抑圧を行うことができる。画像の特定領域の指定は、操作部４４０を介して行われる。これにより、スキャナ９２０で読み取った画像の一部分に対してのみ、ぼけと雑音の抑圧を行い、その部分の画像を復元することができる。【0300】
このように、本実施の形態によれば、画像復元処理部９１０の前処理として、スキャナ９２０からの観測画像信号（ぼけと雑音を含む）からぼけと雑音を抑圧し、クリーンな画像（所望信号）を取得することができるため、画像復元処理部９１０の画像復元能力を最大限に発揮させることができ、非常に高精度の画像復元を実現することができる。【0301】
なお、本実施の形態では、スキャナ９２０で読み取った画像を復元する場合を例にとって説明したが、もちろんこれに限定されない。例えば、スキャナ９２０に代えてデジタルカメラやデジタルビデオカメラなどで撮影した画像を復元する場合にも適用可能である。さらには、すでにある画像情報を取り込んで復元する場合にも適用可能である。【0302】
また、本実施の形態では、画像を復元する場合を例にとって説明したが、もちろんこれに限定されない。本発明は、広く画像処理装置において、何らかの原因でぼけと雑音が入ってしまった劣化画像（観測信号）からぼけと雑音を取り除いた原画像をクリアな信号（所望信号）として取得し、得られた原画像を画像処理する場合に適用可能である。【0303】
２００８年３月２１日出願の特願２００８−０７４６９１の日本出願および２００８年６月２７日出願の特願２００８−１６８８３５の日本出願にそれぞれ含まれる明細書、図面および要約書の開示内容は、すべて本願に援用される。【産業上の利用可能性】
【0304】
本発明に係る雑音抑圧装置および雑音抑圧方法は、よりシンプルな構成で、かつ、より少ない演算量で、所望情報の品質を劣化させることなくより高い雑音抑圧能力を実現することができる雑音抑圧装置および雑音抑圧方法として有用である。【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】
【図25】
【図26】
【図27】
【図28】
【図29】
【図30】
【図31】
【図32】
【図33A-D】
【図33E-H】
【図34A-D】
【図34E-H】
【図35A-D】
【図35E-H】
【図36A-D】
【図36E-H】
【図37】
【図38】
【図39】
【図40】
【図41】
【図42】
【図43】
【図44】
【図45】
【図46】
【図47】
【図48】
【図49】
【図50】