(57)【要約】【課題】ＡＲ係数の推定を必要としないシンプルで雑音抑圧能力が高い、カルマンフィルタに基づく雑音抑圧装置および雑音抑圧方法を提供する。【解決手段】雑音抑圧処理部１５０は、ＡＲ係数を用いないカルマンフィルタを内蔵している。従来手法が線形予測を用いてＡＲ係数を推定した後、その結果を用いてカルマンフィルタを実行することで雑音抑圧を実現しているのに対し、本手法では、カルマンフィルタのみを用いて雑音抑圧を実現するべく、信号源からのクリアな信号のみを用いて状態方程式を構成し、そのクリアな信号と付加雑音信号を用いて観測方程式を構成することで、高い雑音抑圧を実現している。【選択図】図１
2008236270.tif 000002
【特許請求の範囲】
【請求項１】
所望の情報に不必要な情報が混在した情報を取得する取得手段と、
カルマンフィルタのみを用いて、取得された情報から前記不必要な情報を除去して前記所望情報を抽出する抽出手段と、を有し、
前記カルマンフィルタは、
状態空間モデルの状態方程式において自己回帰モデルの係数を使用しないように構成されている、
雑音抑制装置。
【請求項２】
前記カルマンフィルタは、
時刻ｎの所望情報ｄ（ｎ）、不必要情報ｖ（ｎ）、および取得情報ｒ（ｎ）に関して、状態空間モデルの状態ベクトルｘｐを記述する状態方程式および観測ベクトルｙｐを記述する観測方程式が、次の式（１）および式（２）でそれぞれ記述され、
【数１】
2008236270.tif 000053
ここで、Ｋ×Ｋ次の状態遷移行列Φｐおよび観測状態遷移行列Ｍｐは、次の式（３）で表され、
【数２】
2008236270.tif 000054
また、状態ベクトルｘｐ、駆動雑音ベクトルδｐ、および付加雑音ベクトルεｐは、それぞれ、次の式（４）、式（５）、および式（６）で定義される、
【数３】
2008236270.tif 000055
請求項１記載の雑音抑圧装置。
【請求項３】
前記カルマンフィルタは、
時刻ｎの所望情報ｄ（ｎ）、不必要情報ｖ（ｎ）、および取得情報ｒ（ｎ）が、次の式（７）に示す関係を満たし、
【数４】
2008236270.tif 000056
ここで、ΦはＫ×Ｋ次の状態遷移行列であり、次の式（８）で表される、
【数５】
2008236270.tif 000057
請求項１記載の雑音抑圧装置。
【請求項４】
所望の信号に雑音が混在した観測信号を入力する入力手段と、
カルマンフィルタのみを用いて、入力された観測信号から前記雑音を除去して前記所望信号を抽出する抽出手段と、を有し、
前記カルマンフィルタは、
状態空間モデルの状態方程式において自己回帰モデルの係数を使用しないように構成されている、
雑音抑制装置。
【請求項５】
前記カルマンフィルタは、
時刻ｎの所望情報ｄ（ｎ）、不必要情報ｖ（ｎ）、および取得情報ｒ（ｎ）に関して、状態空間モデルの状態ベクトルｘｐを記述する状態方程式および観測ベクトルｙｐを記述する観測方程式が、次の式（１）および式（２）でそれぞれ記述され、
【数６】
2008236270.tif 000058
ここで、Ｋ×Ｋ次の状態遷移行列Φｐおよび観測状態遷移行列Ｍｐは、次の式（３）で表され、
【数７】
2008236270.tif 000059
また、状態ベクトルｘｐ、駆動雑音ベクトルδｐ、および付加雑音ベクトルεｐは、それぞれ、次の式（４）、式（５）、および式（６）で定義される、
【数８】
2008236270.tif 000060
請求項４記載の雑音抑圧装置。
【請求項６】
前記カルマンフィルタの時刻ｎおよび時刻ｎ＋１での推定誤差の共分散Ｐ（ｎ＋１｜ｎ）およびカルマンゲインＫｐが、それぞれ、次の式（９）および式（１０）で定義される、
【数９】
2008236270.tif 000061
ただし、Ｒδｐ、Ｒεｐは、それぞれ、駆動雑音ベクトルおよび付加雑音ベクトルの共分散である、
請求項５記載の雑音抑圧装置。
【請求項７】
前記抽出手段は、
時刻ｎ＋１での駆動ベクトルの共分散を算出する駆動ベクトル共分散算出部と、
状態遷移行列、時刻ｎでの駆動ベクトルの共分散、および前記駆動ベクトル共分散算出部によって算出された時刻ｎ＋１での駆動ベクトルの共分散から、時刻ｎおよび時刻ｎ＋１での推定誤差の共分散行列を算出する推定誤差共分散行列算出部と、
前記推定誤差共分散行列算出部によって算出された時刻ｎおよび時刻ｎ＋１での推定誤差の共分散行列、観測状態遷移行列、ならびに時刻ｎ＋１での付加雑音ベクトルから、時刻ｎ＋１でのカルマンゲインを算出するカルマンゲイン算出部と、
前回算出した時刻ｎでの推定信号ベクトルおよび状態遷移行列から、時刻ｎから時刻ｎ＋１への推定信号ベクトルを算出する第１の推定信号ベクトル算出部と、
前記第１の推定信号ベクトル算出部によって算出された時刻ｎから時刻ｎ＋１への推定信号ベクトル、前記カルマンゲイン算出部によって算出された時刻ｎ＋１でのカルマンゲイン、時刻ｎ＋１での観測ベクトル、および前記観測状態遷移行列から、時刻ｎ＋１での推定信号ベクトルを算出する第２の推定信号ベクトル算出部と、
前記第２の推定信号ベクトル算出部によって算出された時刻ｎから時刻ｎ＋１への推定信号ベクトル、前記カルマンゲイン算出部によって算出された時刻ｎ＋１でのカルマンゲイン、前記観測状態遷移行列、ならびに前記推定誤差共分散行列算出部によって算出された時刻ｎおよび時刻ｎ＋１での推定誤差の共分散行列から、時刻ｎ＋１での推定誤差の共分散行列を算出し更新する推定誤差共分散行列更新部と、
を有する請求項４記載の雑音抑圧装置。
【請求項８】
前記抽出手段は、
時刻ｎ＋１での駆動ベクトルの共分散を算出する駆動ベクトル共分散算出部と、
状態遷移行列、時刻ｎでの駆動ベクトルの共分散、前記駆動ベクトル共分散算出部によって算出された時刻ｎ＋１での駆動ベクトルの共分散から、時刻ｎおよび時刻ｎ＋１での推定誤差の共分散行列を算出する推定誤差共分散行列算出部と、
前記推定誤差共分散行列算出部によって算出された時刻ｎおよび時刻ｎ＋１での推定誤差の共分散行列、観測状態遷移ベクトル、ならびに時刻ｎ＋１でのスカラー量である付加雑音信号から、時刻ｎ＋１でのカルマンゲインを算出するカルマンゲイン算出部と、
前回算出した時刻ｎでの推定信号ベクトルおよび状態遷移行列から、時刻ｎから時刻ｎ＋１への推定信号ベクトルを算出する第１の推定信号ベクトル算出部と、
前記第１の推定信号ベクトル算出部によって算出された時刻ｎから時刻ｎ＋１への信号ベクトル、前記カルマンゲイン算出部によって算出された時刻ｎ＋１でのカルマンゲイン、スカラー量である時刻ｎ＋１での観測信号、および前記観測状態遷移ベクトルから、時刻ｎ＋１での推定信号ベクトルを算出する第２の推定信号ベクトル算出部と、
前記第２の推定信号ベクトル算出部によって算出された時刻ｎから時刻ｎ＋１への推定信号ベクトル、前記カルマンゲイン算出部によって算出された時刻ｎ＋１でのカルマンゲイン、前記観測状態遷移ベクトル、ならびに前記推定誤差共分散行列算出部によって算出された時刻ｎおよび時刻ｎ＋１での推定誤差の共分散行列から、時刻ｎ＋１での推定誤差の共分散行列を算出し更新する推定誤差共分散行列更新部と、
を有する請求項４記載の雑音抑圧装置。
【請求項９】
前記入力手段によって入力された観測信号をサンプリングするサンプリング手段と、
前記サンプリング手段によってサンプリングされた観測信号の振幅強度をディジタル値に変換する変換手段と、
前記変換手段によって変換された観測信号を一時格納し、要求に応じて前記推定手段に出力するバッファ手段と、
をさらに有する請求項４記載の雑音抑圧装置。
【請求項１０】
前記所望信号は、音声信号である、請求項１から請求項９のいずれかに記載の雑音抑圧装置。
【請求項１１】
所望の情報に不必要な情報が混在した情報を取得する取得ステップと、
カルマンフィルタのみを用いて、取得した情報から前記不必要な情報を除去して前記所望情報を抽出する抽出ステップと、を有し、
前記カルマンフィルタは、
状態空間モデルの状態方程式において自己回帰モデルの係数を使用しないように構成されている、
雑音抑制方法。
【請求項１２】
前記カルマンフィルタは、
時刻ｎの所望情報ｄ（ｎ）、不必要情報ｖ（ｎ）、および取得情報ｒ（ｎ）に関して、状態空間モデルの状態ベクトルｘｐを記述する状態方程式および観測ベクトルｙｐを記述する観測方程式が、次の式（１）および式（２）でそれぞれ記述され、
【数１０】
2008236270.tif 000062
ここで、Ｋ×Ｋ次の状態遷移行列Φｐおよび観測状態遷移行列Ｍｐは、次の式（３）で表され、
【数１１】
2008236270.tif 000063
また、状態ベクトルｘｐ、駆動雑音ベクトルδｐ、および付加雑音ベクトルεｐは、それぞれ、次の式（４）、式（５）、および式（６）で定義される、
【数１２】
2008236270.tif 000064
請求項１１記載の雑音抑圧方法。
【請求項１３】
前記カルマンフィルタは、
時刻ｎの所望情報ｄ（ｎ）、不必要情報ｖ（ｎ）、および取得情報ｒ（ｎ）が、次の式（７）に示す関係を満たし、
【数１３】
2008236270.tif 000065
ここで、ΦはＫ×Ｋ次の状態遷移行列であり、次の式（８）で表される、
【数１４】
2008236270.tif 000066
請求項１１記載の雑音抑圧方法。
【請求項１４】
所望の情報に不必要な情報が混在した情報を取得する取得ステップと、
カルマンフィルタのみを用いて、取得した情報から前記不必要な情報を除去して前記所望情報を抽出する抽出ステップと、前記カルマンフィルタは、状態空間モデルの状態方程式において自己回帰モデルの係数を使用しないように構成されている、
をコンピュータに実行させるための雑音抑制プログラム。
【請求項１５】
前記カルマンフィルタは、
時刻ｎの所望情報ｄ（ｎ）、不必要情報ｖ（ｎ）、および取得情報ｒ（ｎ）に関して、状態空間モデルの状態ベクトルｘｐを記述する状態方程式および観測ベクトルｙｐを記述する観測方程式が、次の式（１）および式（２）でそれぞれ記述され、
【数１５】
2008236270.tif 000067
ここで、Ｋ×Ｋ次の状態遷移行列Φｐおよび観測状態遷移行列Ｍｐは、次の式（３）で表され、
【数１６】
2008236270.tif 000068
また、状態ベクトルｘｐ、駆動雑音ベクトルδｐ、および付加雑音ベクトルεｐは、それぞれ、次の式（４）、式（５）、および式（６）で定義される、
【数１７】
2008236270.tif 000069
請求項１４記載の雑音抑圧プログラム。
【請求項１６】
前記カルマンフィルタは、
時刻ｎの所望情報ｄ（ｎ）、不必要情報ｖ（ｎ）、および取得情報ｒ（ｎ）が、次の式（７）に示す関係を満たし、
【数１８】
2008236270.tif 000070
ここで、ΦはＫ×Ｋ次の状態遷移行列であり、次の式（８）で表される、
【数１９】
2008236270.tif 000071
請求項１４記載の雑音抑圧プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、カルマンフィルタに基づく雑音抑圧装置および雑音抑圧方法に関する。
【背景技術】
【０００２】
所望の情報に不必要な情報が混在した受信情報（付加雑音などにより破損した情報）から不必要な情報を取り除き、所望情報のみを抽出することは、音声や無線通信、画像、姿勢制御などの分野における重要な技術であり、近年盛んに研究開発が行われている。
【０００３】
例えば、音声分野における公知の雑音抑圧方法としては、単一のマイクロホンを用いた方法や、複数のマイクロホンから構成されるマイクロホンアレイを用いた方法が提案されている。
【０００４】
しかしながら、マイクロホンアレイを用いた方法では、雑音信号の数が増大すると、マイクロホンの数も比例して増加することが避けらず、コストが増大する。そのため、現在は、単一のマイクロホンを用いた雑音抑圧方法の開発が主流となっている。
【０００５】
単一のマイクロホンしか用いない従来の雑音抑圧方法のアルゴリズムとしては、以下のようなものが知られている。
【０００６】
非特許文献１記載のＡＮＣ（適応ノイズキャンセラ）アルゴリズムは、音声信号の周期性を利用してノイズ信号を低減する。
【０００７】
非特許文献２には、線形予測に基づいた雑音抑圧アルゴリズムが記載されている。このアルゴリズムは、ピッチ推定や、雑音パワースペクトラム、雑音の平均方向に関する事前知識を必要としない。
【０００８】
また、上記アルゴリズムとは別に、カルマンフィルタに基づいた雑音抑圧アルゴリズムが、非特許文献３に提案されている。このアルゴリズムは、音声信号を自己回帰（ＡＲ：AutoRegressive）過程でモデル化する。さらに、このアルゴリズムは、自己回帰モデルの係数（以下「ＡＲ係数」という）を推定し、推定したＡＲ係数を用いたカルマンフィルタに基づいて雑音抑圧を実行する。このアルゴリズムは、ＡＲ係数を推定する必要があるため、上記他のアルゴリズムとは本質的に異なっている。
【０００９】
カルマンフィルタに基づいたアルゴリズムの多くは、通常、２段階で動作する。すなわち、このようなアルゴリズムは、最初にＡＲ係数を推定し、次に推定したＡＲ係数を用いたカルマンフィルタに基づいて雑音抑圧を行う。
【非特許文献１】J.R. Deller, J.G. Proakis, J.H.L. Hansen, "Discrete-Time Processing of Speech Signals," Macmillan Press, 1993
【非特許文献２】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
【非特許文献３】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
【発明の開示】
【発明が解決しようとする課題】
【００１０】
しかしながら、非特許文献１に記載された公知のアルゴリズムは、音声信号のピッチ周期の正確な推定を必要とする。そのため、このアルゴリズムは、その雑音抑圧能力が付加雑音によって劣化してしまうという問題点を有している。
【００１１】
この点、非特許文献２記載のアルゴリズムは、音声信号のピッチ周期の正確な推定を必要とせずに、雑音抑圧が可能である。さらに、このアルゴリズムは、その原理が単純であり、演算量を少なくすることができるといった長所を有している。しかし、このアルゴリズムは、その雑音抑圧能力が入力音声信号の周期性や線形性などの特性に依存している。言い換えると、このアルゴリズムは、雑音抑圧能力の中に音声信号に依存するパラメータが存在しているため、その実用には一定の限界がある。
【００１２】
非特許文献３記載のアルゴリズムは、強力な雑音抑圧能力を有し、特に高い音質を得たい音響分野への応用に適した手法である。
【００１３】
しかしながら、一方で、このアルゴリズムは、ＡＲ係数を必要とするため、ＡＲ係数の推定精度に雑音抑圧性能（つまり、当該カルマンフィルタアルゴリズムの性能）が大きく依存してしまうという問題点を有している。すなわち、ＡＲ係数が正確に推定されない場合、雑音を除去し切れないのみならず、場合によっては雑音に加えて音声信号まで除去してしまう可能性がある。これらは、雑音除去された音声信号の音質の劣化を引き起こす要因となりうる。
【００１４】
この点、一般には、ＡＲ係数の正確な推定は困難である。ＡＲ係数の正確な推定は、例えば、雑音除去であれば、クリアな音声信号に依存しているからである。このことは、音声信号が既知でなければならないことを意味しているため、リアルタイム処理は困難となる。また、仮に何らかの手法でリアルタイムにＡＲ係数を正確に推定することが可能となったとしても、処理が増加するため演算量の問題は避けられない。
【００１５】
本発明は、かかる点に鑑みてなされたものであり、ＡＲ係数の推定を必要としないシンプルで雑音抑圧能力が高い、カルマンフィルタに基づく雑音抑圧装置および雑音抑圧方法を提供することを目的とする。
【課題を解決するための手段】
【００１６】
本発明の雑音抑圧装置は、所望の情報に不必要な情報が混在した情報を取得する取得手段と、カルマンフィルタのみを用いて、取得された情報から前記不必要な情報を除去して前記所望情報を抽出する抽出手段と、を有し、前記カルマンフィルタは、状態空間モデルの観測方程式において自己回帰モデルの係数を使用しないように構成されている、構成を採る。
【００１７】
本発明の雑音抑圧方法は、所望の情報に不必要な情報が混在した情報を取得する取得ステップと、カルマンフィルタのみを用いて、取得した情報から前記不必要な情報を除去して前記所望情報を抽出する抽出ステップと、を有し、前記カルマンフィルタは、状態空間モデルの観測方程式において自己回帰モデルの係数を使用しないように構成されている、ようにした。
【発明の効果】
【００１８】
本発明によれば、カルマンフィルタを用いつつ、ＡＲ係数の推定を必要とすることなく、シンプルな構成で、雑音抑圧能力を向上することができる。
【発明を実施するための最良の形態】
【００１９】
以下、本発明の実施の形態について、図面を参照して詳細に説明する。
【００２０】
（実施の形態１）
図１は、本発明の実施の形態１に係る雑音抑圧装置の構成を示すブロック図である。
【００２１】
図１に示す雑音抑圧装置１００は、入力部１１０、サンプリング部１２０、Ａ／Ｄ変換部１３０、バッファ１４０、雑音抑圧処理部１５０、および出力部１６０を有する。
【００２２】
入力部１１０は、観測信号を入力する。観測信号は、信号源からのクリアな信号と、付加雑音信号とが合わさった信号である。入力部１１０は、例えば、入力したアナログの観測信号を入力処理して、サンプリング部１２０に出力する。入力処理は、例えば、帯域制限処理や自動利得制御処理などである。
【００２３】
サンプリング部１２０は、所定のサンプリング周波数（例えば、１６ｋＨｚ）で、入力されたアナログの観測信号をサンプリング処理し、Ａ／Ｄ変換部１３０に出力する。Ａ／Ｄ変換部１３０は、サンプリングされた観測信号の振幅値を所定の分解能（例えば、８ｂｉｔ）でＡ／Ｄ変換処理し、バッファ１４０に送る。バッファ１４０は、所定のサンプリング数Ｎの信号フレーム（ブロック）を雑音抑圧処理部１５０に出力する。
【００２４】
雑音抑圧処理部１４０は、本発明の特徴的な構成要素であり、ＡＲ係数を用いないカルマンフィルタを内蔵している。すなわち、カルマンフィルタに基づく従来の雑音抑圧方法（非特許文献３参照）（以下「従来手法」という）では、線形予測を用いてＡＲ係数を推定した後、その結果を用いてカルマンフィルタを実行することで雑音抑圧を実現しているのに対し、本発明の雑音抑圧方法（以下「本手法」という）では、カルマンフィルタのみを用いて雑音抑圧を実現している。そのため、本手法では、信号源からのクリアな信号のみを用いて状態方程式を構成し、そのクリアな信号と付加雑音信号を用いて観測方程式を構成している。以下、雑音抑圧処理部１４０で行われるカルマンフィルタに基づく雑音抑圧処理動作について詳細に説明する。
【００２５】
雑音抑圧処理部１４０に入力される観測信号ｒ（ｎ）は、信号源からのクリアな信号（例えば、音声信号など）ｄ（ｎ）以外に付加雑音信号ｖ（ｎ）を含んでおり、次の式（１）を満たす。
【数１】
2008236270.tif 000003
ここで、ｎとは、装置の時刻ｎである。時刻ｎは、サンプリング部１２０で生成された離散的な時間系列において、処理開始時刻を時刻０と仮定したときに、そこからｎ番目の時刻のことを意味する。
【００２６】
従来のＡＲ過程に基づくモデル化方法（従来手法）では、信号ｄ（ｎ）は、ＡＲ係数を用いて次の式（２）でモデル化されると仮定している。
【数２】
2008236270.tif 000004
ここで、αｋ（ｎ）は時刻ｎでのＡＲ係数、ＫはＡＲ係数の次数、ｅ（ｎ）は信号ｄ（ｎ）が、式（２）に示すＫ次のＡＲ過程でモデル化されるとした場合の予測誤差（モデルリング誤差）である。
【００２７】
公知のように、従来手法では、付加雑音信号ｖ（ｎ）は、零平均であり、白色雑音であることが前提条件である。言い換えると、従来手法では、信号ｄ（ｎ）と付加雑音信号ｖ（ｎ）は無相関であり、つまり、次の式（３）を満たす。
【数３】
2008236270.tif 000005
ここで、Ｅ［・］はアンサンブル平均を示す。
【００２８】
ＡＲ係数は、ＡＲ係数推定アルゴリズムにより推定される。従来手法において最も重要な点は、カルマンフィルタを用いた高性能の雑音抑圧を達成するために、ＡＲ係数の正確な推定を必要とすることである。このことからも、カルマンフィルタの雑音抑圧能力がＡＲ係数の推定精度に大きく依存しているため雑音抑圧能力が大きく劣化することは容易に想像可能である。
【００２９】
本手法では、ＡＲ係数を用いずにカルマンフィルタの状態空間モデル、つまり音源信号からのクリアな信号のみを用いて状態方程式、およびそのクリアな信号と付加雑音信号とを用いて観測方程式を構成している。
【００３０】
以下、本手法での状態空間モデルの構成法について説明する。表記を容易にするために、まずＫ×１次の信号ベクトルｘｐ（ｎ）を次の式（４）で定義する。添え字“ｐ”は本発明により考案された表現であることを示す。
【数４】
2008236270.tif 000006
【００３１】
次に、式（２）と同様にして信号ｄ（ｎ）のモデルを構成する。ここでは、ＡＲ係数を用いないとしても、従来手法と同様に、未知のＫ×Ｋ次の状態遷移行列Φｐ（ｎ＋１）と、Ｋ×１次の駆動雑音ベクトルδｐ（ｎ＋１）を導入する。これにより、本実施の形態の状態空間モデル（信号ベクトルにより記述される）の状態方程式の形として、次の式（５）が定まる。
【数５】
2008236270.tif 000007
【００３２】
式（５）が等式としてＡＲ係数を用いずに成立するためには、Ｋ×Ｋ次の状態遷移行列Φｐ（ｎ＋１）、Ｋ×１次の駆動雑音ベクトルδｐ（ｎ＋１）は、次の式（６）および式（７）を満たすことが求められる。
【数６】
2008236270.tif 000008
【数７】
2008236270.tif 000009
【００３３】
したがって、式（５）は次の式（８）に書き直せる。
【数８】
2008236270.tif 000010
【００３４】
式（６）で表されるＫ×Ｋ次の状態遷移行列Φｐ（ｎ＋１）の行列要素は０と１のみであり、特に第１行がすべて０である。このことは、本手法の状態空間モデルの特徴の一つである。また、上記から明らかなように、式（５）、式（６）、式（７）の導出において、駆動雑音ベクトルδｐ（ｎ＋１）に対して前提条件を与えていないことに注意すべきである。すなわち、駆動雑音ベクトルδｐ（ｎ＋１）は有色であってよい。この理由については後述する。
【００３５】
本手法のカルマンフィルタアルゴリズムでは、特異値分解を避けるためにＫ×１次の観測ベクトルｙｐ（ｎ＋１）を次の式（９）で定義する。
【数９】
2008236270.tif 000011
【００３６】
よって、式（１）および式（９）から、本手法の状態空間モデルの観測方程式として、次の式（１０）が導出される。
【数１０】
2008236270.tif 000012
ここで、ＭｐはＫ×Ｋ次の状態遷移行列、εｐはＫ×１次の付加雑音ベクトルであり、それぞれ次の式（１１）および式（１２）で定義される。
【数１１】
2008236270.tif 000013
【数１２】
2008236270.tif 000014
【００３７】
以上より、本手法のカルマンフィルタアルゴリズムで用いる状態空間モデルの状態方程式および観測方程式として、式（５）および式（１０）がそれぞれ求められた。
【００３８】
図２は、本実施の形態に係る雑音抑圧装置の状態空間モデルを説明するブロック線図である。
【００３９】
図２において、５００は時刻ｎにおける信号ベクトルｘｐ（ｎ）、５０１は時刻ｎ＋１における信号ベクトルｘｐ（ｎ＋１）、５０２は時刻ｎにおける観測ベクトルｙｐ（ｎ）、５０３は時刻ｎにおける付加雑音ベクトルεｐ（ｎ）、５０４は時刻ｎ＋１における駆動雑音ベクトルδｐ（ｎ＋１）、５０５は状態遷移行列Φｐ、５０６は状態遷移行列Ｍｐである。
【００４０】
図２から理解できるように、本手法のカルマンフィルタアルゴリズムは、駆動雑音が有色であるにもかかわらず、従来のカルマンフィルタアルゴリズムと同様の手順で実行できる（理由は後述する）。本手法のカルマンフィルタアルゴリズムに基づいた信号推定は、後述するフローチャートに従って行われる。
【００４１】
図３は、本実施の形態に係る雑音抑圧装置のカルマンフィルタアルゴリズムに基づく信号推定手順を説明するフローチャートである。
【００４２】
まず、推定信号ベクトルの初期値ｘｅ（０｜０）、推定誤差の共分散行列の初期値Ｐ（０｜０）およびスカラー量である付加雑音ベクトルの共分散Ｒεｐ（ｎ）［ｉ，ｊ］の値を、次の式（１３）に示すように設定する（ＳＴ１００）。
【数１３】
2008236270.tif 000015
【００４３】
ここで、Ｉは単位行列である。また、σｖ２は付加雑音信号ｖ（ｎ）の雑音分散であり、既知と仮定している。もし付加雑音信号ｖ（ｎ）が白色雑音であり零平均であれば、σｖ２は、次の式（１４）で与えられる。
【数１４】
2008236270.tif 000016
ここで、Ｎは所定のサンプル数である。
【００４４】
次に、カルマンフィルタアルゴリズムを実行する（ＳＴ１０１）。このカルマンフィルタアルゴリズムの処理手順は、後で詳細に説明する。
【００４５】
次に、時刻ｎが所定のサンプル数Ｎに達したか否かを判定し（ＳＴ１０２）、所定のサンプル数Ｎに達していない場合は（ＳＴ１０２：ＮＯ）、ステップＳＴ１０１に戻って、カルマンフィルタアルゴリズムを繰り返し、所定のサンプル数Ｎに達した場合は（ＳＴ１０２：ＹＥＳ）、上記一連の処理を終了する。
【００４６】
本手法では、ＡＲ係数を用いずに状態空間モデルを設定している。そのため、従来手法で必要であったＡＲ係数を推定するステップを削減することができ、１段階処理で信号推定が可能となる。このことは本発明の大きな特徴の一つである。
【００４７】
図４は、図３のカルマンフィルタアルゴリズムの処理内容を示すフローチャートである。
【００４８】
ここで、開始時の時刻はｎ＋１であり、以下のフローでは、時刻ｎ＋１での観測信号ｒ（ｎ＋１）および観測ベクトルｙｐ（ｎ＋１）と、時刻ｎでの推定誤差の共分散行列Ｐ（ｎ｜ｎ）とから、時刻ｎ＋１での推定信号ベクトルｘpe（ｎ＋１｜ｎ）を推定し、同時に推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ＋１）を更新する。
【００４９】
まず、次の式（１５）を用いて、駆動雑音ベクトルの共分散Ｒδｐ（ｎ）［ｉ，ｊ］の値を計算する（ＳＴ２００）。
【数１５】
2008236270.tif 000017
【００５０】
次に、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）を計算し、更新する（ＳＴ２０１）。ここで、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）の更新式は、次の式（１６）である。ここでは、以下、添え字“ｐｅ”および“ｅ”は、推定値であることを示す。
【数１６】
2008236270.tif 000018
【００５１】
この式（１６）を展開すると、次の式（１７）となる。
【数１７】
2008236270.tif 000019
【００５２】
ここで、式（１６）の第３項および第４項に注目する。第３項および第４項は、推定信号ベクトルｘｐｅ（ｎ｜ｎ）と駆動雑音ベクトルδｐ（ｎ｜ｎ）のアンサンブル平均を含んでいる、したがって、推定信号ベクトルｘｐｅ（ｎ｜ｎ）と駆動雑音ベクトルδｐ（ｎ｜ｎ）とが無相関でない、言い換えれば、駆動雑音ベクトルδｐ（ｎ｜ｎ）が有色のとき、第３項および第４項は０にならない。すなわち、本手法のカルマンフィルタアルゴリズムが、駆動雑音信号が有色信号の場合に、従来と同様なカルマンフィルタアルゴリズム手順で適用可能か否かは、第３項および第４項が、式（１５）で表される推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）の更新にどの程度影響するかによる。
【００５３】
そこで、式（１７）の第３項に着目する。第３項は、次の式（１８）のように、
【数１８】
2008236270.tif 000020
と書き表せる。
【００５４】
ここで、次の式（１９）のように、
【数１９】
2008236270.tif 000021
とすれば、次の式（２０）となる。
【数２０】
2008236270.tif 000022
【００５５】
式（５）および次の式（２１）
【数２１】
2008236270.tif 000023
を考慮すれば、式（１８）は、次の式（２２）のように書き直すことができる。
【数２２】
2008236270.tif 000024
ただし、次の式（２３）である。
【数２３】
2008236270.tif 000025
ここで、ｄｅ（ｎ）は時刻ｎでの推定信号である。
【００５６】
式（２３）の第１項は真値のアンサンブル平均、第２項は推定値のアンサンブル平均を示す。したがって、その差は小さく、ｅ（ｌ＋１）は無視できる値である。よって、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）の更新式（１５）は、次の式（２４）のように書き直される。このことから、駆動雑音が有色であっても、本手法の式（５）および式（９）から構成される状態空間モデル(状態方程式と観測方程式)であれば、従来手法と同様なカルマンフィルタアルゴリズム手順が適用可能となる。
【数２４】
2008236270.tif 000026
【００５７】
次に、カルマンゲインＫｐ（ｎ＋１）を計算する（ＳＴ２０２）。ここで、カルマンゲインＫｐ（ｎ＋１）は、次の式（２５）で与えられる。
【数２５】
2008236270.tif 000027
【００５８】
そして、これらの値から、まず推定信号ベクトルｘｐｅ（ｎ＋１｜ｎ）を計算する（ＳＴ２０３）。ここで、推定信号ベクトルｘｐｅ（ｎ＋１｜ｎ）は、次の式（２６）で与えられる。
【数２６】
2008236270.tif 000028
【００５９】
そして、推定信号ベクトルｘｐｅ（ｎ＋１｜ｎ＋１）を計算する（ＳＴ２０４）。ここで、推定信号ベクトルｘｐｅ（ｎ＋１｜ｎ＋１）は、次の式（２７）で与えられる。
【数２７】
2008236270.tif 000029
【００６０】
最後に、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ＋１）を計算、更新し（ＳＴ２０５）、フローを終了する。ここで、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ＋１）は、次の式（２８）で与えられる。
【数２８】
2008236270.tif 000030
【００６１】
雑音抑圧処理部１４０は、上述した手順、アルゴリズムによって推定された推定信号をを出力部１５０に出力する。出力部１５０は、例えば、スピーカやディスプレイ、通信手段、記憶装置などで構成されている。
【００６２】
ＡＲ係数の推定を必要とする従来のカルマンフィルタの問題点は、ＡＲ係数の推定精度にカルマンフィルタアルゴリズムの能力が依存していることである。これに対して、本手法は、カルマンフィルタアルゴリズムのみで実行しているため、ＡＲ係数の推定精度に依存しない。
【００６３】
なお、本手法のカルマンアルゴリズムは、以下に説明する更に改良された方法（以下「改良手法」という）によっても実行できる。以下、この改良手法について説明する。
【００６４】
改良手法のカルマンアルゴリズムの状態方程式は、上記式（５）と同一である。ただし、観測方程式を、次の式（２９）のように書き直す。改良手法の観測方程式（２９）は、本手法の観測方程式（１０）の一部の要素を取り出してきたものである。
【数２９】
2008236270.tif 000031
ここで、ｙ’ｐ（ｎ＋１）、ｍＴｐ、ε’ｐ（ｎ＋１）は、それぞれ改良手法における観測信号、雑音信号、状態遷移ベクトルであり、次の式（３０）を満たす。
【数３０】
2008236270.tif 000032
【００６５】
改良手法の大きな特徴の一つは、本手法の観測方程式（１０）と改良手法の観測方程式（３０）の比較からわかるように、改良前の方法（本手法）における、観測ベクトルｙｐ（ｎ＋１）、雑音ベクトルεｐ（ｎ＋１）、および状態遷移行列ｍｐが、それぞれスカラーとベクトルに変更されていることである。したがって、改良手法はより少ない演算量で実行できるであろうことは明らかである。また、上記から明らかなように、改良手法においても、本手法と同様に、駆動雑音ベクトルδｐ（ｎ＋１）に対して前提条件を与えていないことに注意すべきである。すなわち、改良手法においても駆動雑音ベクトルδｐ（ｎ＋１）は有色であってよい。
【００６６】
すなわち、式（１７）で与えられた推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）を再度、行列形式で書き下すと、次の式（３１）となる。
【数３１】
2008236270.tif 000033
ここで、また以下の式で、行列式中、斜線網掛け部分で示される部分は駆動雑音が有色であるときに影響を受ける部分である。
【００６７】
改良手法を用いたカルマンフィルタアルゴリズムは、図４に示す本手法のカルマンフィルタアルゴリズムとは、カルマンゲインＫｐ（ｎ＋１）の計算（ＳＴ２０２）および推定信号ベクトルｘｐｅ（ｎ＋１｜ｎ＋１）の計算（ＳＴ２０４）における計算式が異なる。
【００６８】
すなわち、カルマンゲインＫｐ（ｎ＋１）の計算式は、次の式（３２）、つまり、
【数３２】
2008236270.tif 000034
となる。ただし、雑音信号ε’ｐ（ｎ＋１）の共分散をσ２ｖとしたときに、次の式（３３）が成り立つ。
【数３３】
2008236270.tif 000035
【００６９】
また、推定信号ベクトルｘpｅ（ｎ＋１｜ｎ＋１）の計算式は、次の式（３４）となる。
【数３４】
2008236270.tif 000036
となる
【００７０】
ここで、式（３４）のｘpｅ（ｎ＋１｜ｎ＋１）の第１要素に着目すると、斜線網掛けではない（つまり、有色雑音の影響を受けない）ことから、駆動雑音が有色であってもクリアな信号の推定信号を得られる。したがって、改良手法も、本手法と同様に、従来のカルマンフィルタアルゴリズムで実行可能である。既に述べたように、改良手法では、本手法に比べて、演算量を低減できる。しかしながら、本手法が観測ベクトルｙｐ（ｎ＋１）を用いているのに対して、改良手法方法では、スカラー量である観測信号ｙ‘ｐ（ｎ＋１）を用いている。そのサイズの違いから、過去の観測量をより積極的に用いることができる本手法のほうが、より雑音抑圧能力が高いことは明らかである。
【００７１】
以上をまとめると、本発明における上記二つの手法（本手法および改良手法）は、カルマンフィルタに必要な状態空間モデルの観測方程式を変更することによって、演算量を大幅に減少させることが可能である。より具体的には、本発明における上記二つの手法では、ＡＲ係数の推定を必要としないため、その演算量は従来手法に比べて少ない。また、改良手法では、演算に用いる変数のサイズが本手法に比べて小さいため、その演算量は本手法に比べて少ない。すなわち、演算量に関して、
改良手法＜本手法＜従来手法
である。
【００７２】
また、本発明における上記二つの手法を、例えば、半導体集積回路や半導体回路などのハードウエアとして実施する場合、また、パーソナルコンピュータなどで実行可能なソフトウエアとして実施する場合のいずれにおいても、その構成は従来手法よりも単純化される。したがって、本発明における二つの手法を用いれば、回路規模やプログラム量が大幅に低減できるであろうことは明らかである。
【００７３】
また、本発明における上記二つの手法は、カルマンフィルタに必要な状態空間モデルの観測方程式を変更することによって、雑音抑圧能力を大幅に向上させることが可能である。より具体的には、本発明における上記二つの手法では、ＡＲ係数の推定を必要としないため、雑音抑圧能力がＡＲ係数の推定精度に依存する従来手法に比べて、雑音抑圧能力が高い。また、改良手法では、演算に用いる変数のサイズが本手法に比べて小さいため、過去の観測量を積極的に用いていない。したがって、その雑音抑圧能力は本手法に比べて低い。すなわち、雑音抑圧能力に関して、
従来手法＜改良手法＜本手法
である。
【００７４】
例えば、音響分野において、音質を若干落としても（ただし、従来手法よりは高い）、演算速度を速くしたい場合に、改良手法は有効である。
【００７５】
ＡＲ係数を用いない本発明における上記二つの手法は、状態方程式の駆動雑音が有色信号となっており、従来のカルマンフィルタの仮定、つまり、駆動雑音が白色信号であるという仮定に反している。しかしながら、本発明の状態空間モデルにおける状態方程式と観測方程式の構造的な性質により、駆動雑音が白色信号である従来のカルマンフィルタアルゴリズと全く同じアルゴリズムで実行可能であり、その理由は上記の通りである。また、本発明の改良手法においても、同様に実行可能であることは上記より明らかである。
【００７６】
本発明における上記二つの方法は、ＡＲ係数を使わないようにカルマンフィルタを構成する。言い換えると、カルマンフィルタのみで実行可能な雑音抑圧方法であるため、従来の雑音抑圧能力の問題を解決していることは明らかである。
【００７７】
（実施の形態２）
実施の形態２は、実施の形態１に示す雑音抑圧装置を音声に適用した場合である。
【００７８】
図５は、本発明の実施の形態２に係る雑音抑圧装置の構成を示すブロック図である。
【００７９】
図５に示す雑音抑圧装置２００は、本実施の形態の雑音抑圧処理を実行できるパーソナルコンピュータ２１０、マイクロホン２２０、サンプリング部１２０、およびＡ／Ｄ変換部１３０を有する。
【００８０】
パーソナルコンピュータ２１０は、操作装置２１１、ディスプレイ２１２、バスインタフェース２１３、記録装置２１４、主記憶メモリ２１５、および中央演算装置２１６を有する。操作装置２１１は、典型的にはキーボートやマウスなどであるが、音声認識装置などを用いてもよい。使用者は、操作装置２１１を用い、ディスプレイ２１２で確認をしながらコンピュータを操作できる。
【００８１】
パーソナルコンピュータ２１０において本実施の形態の雑音抑圧処理を実行させるプログラムソフトウエアは、記録装置２１４に格納されていてもよいし、バスインタフェース２１３を介して外部からダウンロードされてきてもよい。記録装置２１４は、典型的にはハードディスク装置であるが、ＣＤ−ＲＯＭ装置やＤＶＤ装置、フラッシュメモリなどの可搬性のあるものであってもよい。また、それらの組み合わせであってもよい。
【００８２】
サンプリング部１２０およびＡ／Ｄ変換部１３０は、パーソナルコンピュータ２１０内部に格納された内蔵カード（ボード）であってもよいし、バスインタフェース２１３を経由して接続された外部設置型機器であってもよい。
【００８３】
マイクロホン２２０からの観測音声信号は、サンプリング部１２０に入力される。サンプリング部１２０は、所定のサンプリング周波数（例えば、１６ｋＨｚ）で、入力されたアナログの観測音声信号をサンプリング処理し、Ａ／Ｄ変換部１３０に出力する。Ａ／Ｄ変換部１３０は、サンプリングされた観測音声信号の振幅値を所定の分解能（例えば、８ｂｉｔ）でＡ／Ｄ変換処理し、一時格納する。Ａ／Ｄ変換部１３０は、所定のサンプリング数Ｎの音声フレーム単位で、ディジタル化した観測音声信号をパーソナルコンピュータ２１０のバスインタフェース２１３に出力する。
【００８４】
パーソナルコンピュータ２１０はバスインタフェース２１３に出力された観測音声信号を一時、主記憶メモリ２１５に格納し、その後、所定の音声フレーム（サンプリング数）単位で、雑音抑圧処理を施した上で、再度主記憶メモリ２１５に格納する。雑音抑圧処理は、主記憶メモリ２１５や記録装置２１４に格納されたソフトウエアをバスインタフェース２１３経由で中央演算装置２１６に呼び出し、実行させることで行われる。
【００８５】
パーソナルコンピュータ２１０は、使用者の操作により、処理を実行したり、中断、終了させたりする。また使用者の操作により、処理した推定音声信号を外部に出力してもよい。
【００８６】
次に、本実施の形態における音声の雑音抑圧処理について、図面を参照しつつ説明する。
【００８７】
本実施の形態における音声の雑音抑圧の性能評価の目的で音声波形の数値シミュレーションを行った。
【００８８】
図６は、本実施の形態における音声波形シミュレーションの第１の例の結果を示す図であり、図７は、本実施の形態における音声波形シミュレーションの第２の例の結果を示す図である。また、図８は、本実施の形態における音声波形シミュレーションの第３の例の結果を示す図であり、図９は、本実施の形態における音声波形シミュレーションの第４の例の結果を示す図である。
【００８９】
日本人成人男性および女性のオリジナル音声信号（クリアな音声信号）は、無響室において、１６ｋＨｚのサンプリングレートでサンプリング、ディジタル化した。本例では、二つの音声信号サンプルを検討した。二つの音声信号サンプルは、（Ａ−１）図６（Ａ）に示す無声区間のないクリアな音声信号、（Ａ−２）図７（Ａ）に示す無声区間をもつクリアな音声信号である。それぞれの音声信号サンプルは、以後、音声（Ａ−１）および音声（Ａ−２）と参照する。
【００９０】
雑音信号は、人工的な付加雑音信号であり、本例では、二つの雑音信号サンプルを検討した。二つの雑音信号サンプルは、（Ｂ−１）図６（Ｂ）および図７（Ｂ）に示す付加白色ガウス雑音、（Ｂ−２）図８（Ｂ）および図９（Ｂ）に示す付加有色バブル雑音である。それぞれの雑音信号サンプルは、以後、雑音（Ｂ−１）および雑音（Ｂ−２）と参照される。
【００９１】
雑音信号の雑音分散σｖ２は既知として次の式（３５）で表されるものとする。すなわち、
【数３５】
2008236270.tif 000037
である。ここで、Ｌは全音声信号サンプル数である。
【００９２】
また、信号雑音比ＳＮＲｉｎを式（３６）で定義する。
【数３６】
2008236270.tif 000038
【００９３】
従来手法と本発明の方法（本手法）とによる雑音抑圧の結果を、音声（Ａ−１）および音声（Ａ−２）と雑音（Ｂ−１）の組み合わせ、ならびに音声（Ａ−１）および音声（Ａ−２）と雑音（Ｂ−２）の組み合わせからなる観測音声信号の信号波形から比較した。
【００９４】
図６（Ｃ）と図７（Ｃ）は、それぞれ音声（Ａ−１）と雑音（Ｂ−１）、音声（Ａ−２）と雑音（Ｂ−１）からなる観測音声信号波形である。
【００９５】
図６（Ｄ）と図７（Ｄ）は、それぞれ音声（Ａ−１）と雑音（Ｂ−１）、音声（Ａ−２）と雑音（Ｂ−１）の合成波形に対して従来手法で雑音抑圧を行った後の推定音声信号波形である。
【００９６】
一方、図６（Ｅ）と図７（Ｅ）は、それぞれ音声（Ａ−１）と雑音（Ｂ−１）、音声（Ａ−２）と雑音（Ｂ−１）の合成波形に対して従来手法で雑音抑圧を行った後の推定音声信号波形である。
【００９７】
図６（Ｄ）と図７（Ｄ）に対して図６（Ａ）と図７（Ａ）に示すオリジナルの音声信号とを参照すると、図６（Ｄ）と図７（Ｄ）に対して図６（Ａ）と図７（Ａ）に示すクリアな音声信号とを比較すると、従来手法による雑音抑圧では、雑音抑圧後に推定音声信号の振幅が小さくなっており、クリアな音声信号が抑圧されていることがわかる。加えて、従来手法による雑音抑圧では、サンプリング数の増加とともに、雑音抑圧後の推定音声信号の波形が、図６（Ａ）と図７（Ａ）に示すクリアな音声信号の波形から変形していく。
【００９８】
さらに、従来手法の雑音抑圧では、無声区間を有する音声（Ａ−２）に対して、推定音声信号が抑圧されるだけでなく、無声区間においてオリジナルの雑音信号と異なる雑音が観察されている。これは、従来手法では、無声区間では信号ｄ（ｎ）がゼロであるにもかかわらず、式（２）でＡＲ係数を求めようとするため、ＡＲ係数の値は発散し、不安定な状態を与えるからと推測される。
【００９９】
また、このことから、雑音信号が有色の場合、従来手法の適用は困難であろうことは容易に推測される。
【０１００】
従来手法とは対照的に、図６（Ｅ）および図７（Ｅ）に示す本発明の雑音抑圧では、雑音抑圧後の推定音声信号の波形が、図６（Ａ）と図７（Ａ）に示すクリアな音声信号の波形と非常に似通っている。
【０１０１】
次に、図８（Ｄ）と図９（Ｄ）に対して図８（Ａ）と図９（Ａ）に示すクリアな音声信号とを比較すると、従来手法による雑音抑圧では、雑音（Ｂ−２）を含む観測音声信号に対して、非常に劣った結果を与えていることがわかる。これは、従来手法では、有色雑音である雑音（Ｂ−２）を含んだ観測音声信号に対してＡＲ係数を正確に推定することが困難であるためである。
【０１０２】
一方、本発明の雑音抑圧方法では、雑音（Ｂ−２）の場合も雑音（Ｂ−１）の場合と同程度の雑音抑圧が達成されている。
【０１０３】
すなわち、上述のように、本発明の雑音抑制方法は、白色、有色雑音、無声区間の有無に関わらず有効である。これは、本発明の雑音抑制方法の大きな特徴の一つである。
【０１０４】
本発明の雑音抑制方法が有色の付加雑音にも適用できることは、既述のカルマンアルゴリズムにおける推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）の更新式（２３）の導出過程からも明らかであるが、上述の本実施の形態の数値シミュレーションは、このことを支持している。
【０１０５】
次に、雑音抑制性能を数値的に比較するため、雑音抑圧性能ＳＮＲｏｕｔを式（３７）で定義し、その数値シミュレーションを行った。
【数３７】
2008236270.tif 000039
ここで、ｄｅｉ（ｎ）は推定音声信号を表す。
【０１０６】
表１は、本実施の形態における雑音抑圧性能の数値シミュレーションの一例の結果を示す表であり、音声（Ａ−１）と雑音（Ｂ−１）の組み合わせにおける、幾つかのＳＮＲｉｎ、Ｋ（従来手法においては式（２）におけるＡＲ係数の次数、本発明では状態遷移行列Φｐの次数）の値における、従来手法と本発明の方法の雑音抑圧性能ＳＮＲｏｕｔを比較して示している。
【０１０７】
【表１】
2008236270.tif 000040
【０１０８】
表２は、本実施の形態における雑音抑圧性能の数値シミュレーションの他の例の結果を示す表であり、音声（Ａ−１）と雑音（Ｂ−１）の組み合わせにおける、幾つかのＳＮＲｉｎ、Ｋの値における、従来手法と本発明の方法の雑音抑圧性能ＳＮＲｏｕｔを比較して示している。
【０１０９】
【表２】
2008236270.tif 000041
【０１１０】
表１および表２を参照すると、本発明の方法は、すべてのＳＮＲｉｎ、Ｋの値において、従来の方法に比べて雑音抑圧能力を改善していることがわかる。
【０１１１】
特に、表２に示す有色雑音の場合には、従来手法は非常に劣った結果を与えているのにも関わらず、本発明の方法は表１に示す白色雑音の場合と同程度の結果を示している。すなわち、本発明の雑音抑圧方法は、白色雑音、有色雑音両者に効果的で雑音の性質に堅牢な雑音抑圧方法であるといえる。
【０１１２】
また、表１および表２に見られるように、本発明の方法では、Ｋの値に対して雑音抑圧性能ＳＮＲｏｕｔは安定であり、Ｋの値の増加に伴い増加する傾向にある。対照的に従来の方法では、表１に見られるように、Ｋの値に対して雑音抑圧性能ＳＮＲｏｕｔは不安定である。これは、従来手法では、最適なＫの値、つまりＡＲ係数の次数を決定することが困難であることを意味している。
【０１１３】
ＡＲ係数推定を必要とする従来手法において最も問題になることは、ＡＲ係数の次数の正確な推定は一般に困難であることである。なぜなら、ＡＲ係数の次数の正確な推定は，例えば雑音除去であれば、クリアな音声信号に依存しているからである。
【０１１４】
このことは、クリアな音声信号が既知でなければならないことを意味しているため、リアルタイム処理は困難となる。ＡＲ係数の次数が正確でない場合には、カルマンフィルタアルゴリズの性能が劣化することは容易に想像可能である。また、何らかの手法でリアルタイムに推定することが可能となったとしても、処理が増加することより演算量などの問題を避けることは不可能である。
【０１１５】
次に、本実施の形態に係る推定音声信号の音声品質を評価するためにリスニングテストによる主観的評価を行った。
【０１１６】
音声品質評価に用いたオリジナル音声信号とオリジナル雑音信号は第２の実施の形態と同一であり、その説明は省略する。雑音信号は、異なるＳＮＲｉｎ（０、５、および１０［ｄＢ］）で音声信号に加えた。
【０１１７】
音声品質評価は、ＡＣＲ（絶対範疇評価）に基づいた５段階ＭＯＳ（平均オピニオン値）を用いたリスニングテストにより行った。５０人の聴取者が雑音抑圧により得られた推定音声信号のうち幾つかを評価した。各々の聴取者は、ポイント１からポイント５を決定する。ポイント５が最良である。
【０１１８】
図１０は、本実施の形態における雑音抑圧後の音声品質の主観的評価結果の一つの例を示す図であり、音声（Ａ−１）と雑音（Ｂ−１）の組み合わせにおける、従来手法と本発明の方法とリスニングテストの結果を比較して示している。
【０１１９】
また、図１１は、本実施の形態における雑音抑圧後の音声品質の主観的評価結果の他の例を示す図であり、音声（Ａ−１）と雑音（Ｂ−２）の組み合わせにおける、従来手法と本発明の方法とリスニングテストの結果を比較して示している。
【０１２０】
図１０および図１１から、本発明の方法で推定した音声信号のスコアは、すべてのＳＮＲｉｎ値において従来手法のスコアより高いことがわかる。特にその差は、音声（Ａ−１）と雑音（Ｂ−２）の組み合わせに対して大きい。
【０１２１】
以上より、本発明の雑音抑圧方法は、音声信号の音声品質を犠牲にすることのない、白色雑音、有色雑音に効果的な優れた雑音抑圧方法であるといえる。
【０１２２】
（実施の形態３）
本発明に係る雑音抑圧装置は、雑音抑圧以外の用途にも応用可能である。以下、その一つの応用例について説明する。
【０１２３】
図１２は、本発明の雑音抑圧方法、つまりＡＲ係数を必要としないカルマンフィルタが適用されたマルチキャリア受信装置の構成を示すブロック図である。
【０１２４】
図１２に示すマルチキャリア受信装置３００は、主に、検波部３０１、ＧＩ（ガードインターバル）除去部３０２、チャネル等化部３０３、チャネル推定部３０４、ＦＦＴ部３０５、復調部３０６、および復号部３０７を有する。
【０１２５】
本実施の形態において、チャネル推定部３０４は、本発明の状態空間モデルを用いたカルマンフィルタに基づいたチャネル推定を実行できるように構成されている。より具体的には、チャネル推定部３０４は、サンプリング部１２０、Ａ／Ｄ変換部１３０、バッファ１４０、雑音抑圧処理部１５０、出力部１６０を有する構成をとる。
【０１２６】
検波部３０１は伝送路上で周波数選択性フェージング等の影響を受けた受信信号を、中間周波数で直交検波し、ＧＩ除去部４０２に出力する。ＧＩ除去部３０２は、受信信号のガードインターバルを除去し、シンボル単位に連なった信号をチャネル等化部３０３と、チャネル推定部３０４に出力する。
【０１２７】
チャネル推定部３０４に入力された信号は、入力部１１０に出力される。入力部は信号に所定の入力信号処理を施し、サンプリング部１２０に出力する。サンプリング部１２０は、所定のサンプリング周波数（例えば１６ｋＨｚ）で、入力されたアナログの受信信号をサンプリング処理し、Ａ／Ｄ変換部１３０に出力する。Ａ／Ｄ変換部１３０は、サンプリングされた受信信号の振幅値を所定の分解能（例えば８ｂｉｔ）でＡ／Ｄ変換処理し、バッファ１４０に送る。バッファ１４０は所定のサンプリング数Ｎの受信信号フレーム（ブロック）を雑音抑圧処理部１５０に出力する。
【０１２８】
雑音抑圧処理部１５０は、カルマンフィルタに基づいたチャネル推定をサブキャリア全体、または各サブキャリアついて行い、推定された（サブ）チャネルゲインを出力部１６０に出力する。出力部１６０は、入力された推定チャネルゲインをチャネル等化部３０２に出力する。雑音処理部１５０での処理の詳細については後述する。
【０１２９】
チャネル等化部３０３は、入力された推定チャネルゲインを用いて、ＧＩ除去部３０２から入力した信号の同期検波を行い、結果をＦＦＴ部３０５に出力する。ＦＦＴ部３０５は、フーリエ変換処理を行い、受信信号を各サブキャリア信号成分に分離して、復調部３０６に出力する。復調部３０６は信号の復調処理を行い、結果を復号部３０７に出力する。復号部３０７は、信号の復号を行い、デジタルデータを出力する。
【０１３０】
本実施の形態では、まず、フェージングチャンネルを既述の状態空間モデルで表現する。マルチキャリア通信において、第ｋサブキャリアのみに着目すると、その受信信号ｙｋ（ｔ）は、次の式（３８）
【数３８】
2008236270.tif 000042
となる。ここでｎは時刻、Ｎは総キャリア数、Ｓｋ（ｔ）、Ｈｋ（ｔ）は第ｋキャリアの送信信号、チャネルゲインである。
【０１３１】
表記を容易にするために、Ｎ×１次チャネルゲインベクトルｈp（ｎ）を次の式（３９）で定義する。
【数３９】
2008236270.tif 000043
ここで、Ｋは後述の状態遷移行列の次数である。
【０１３２】
本実施の形態では、本発明の状態空間モデルを用いてカルマンフィルタに基づき、チャネルゲインの推定値を適応的に求める。状態空間モデルの状態ベクトルｘｐ（ｎ＋１）を次の式（４０）で定義する。
【数４０】
2008236270.tif 000044
【０１３３】
伝送路特性の時間変動を表す状態方程式は、次の式（４１）で記述される。
【数４１】
2008236270.tif 000045
また、観測方程式は、次の式（４２）で記述される。
【数４２】
2008236270.tif 000046
ここで、δｐ（ｎ＋１）は駆動雑音ベクトル、εｐ（ｎ＋１）は雑音信号ベクトル、Ｃｐ、Ｄｐは状態遷移行列、Ｇは定数行列である。
【０１３４】
駆動雑音ベクトルδｐ（ｎ＋１）、雑音信号ベクトルεｐ（ｎ＋１）、状態遷移行列Ｃｐ、Ｄｐ、定数行列Ｇは以下の式（４３）から式（４８）で定義される。
【数４３】
2008236270.tif 000047
【数４４】
2008236270.tif 000048
【数４５】
2008236270.tif 000049
【数４６】
2008236270.tif 000050
【数４７】
2008236270.tif 000051
【数４８】
2008236270.tif 000052
ここで、ＩＮ×Ｎ、０Ｎ×ＮはＮ×Ｎの単位行列、零行列である。
【０１３５】
上記の定式化は、従来Vector Kermanと呼ばれているチャネルゲイン推定法に、本発明の状態空間モデルを適用したものである。フェージングチャネル間の相関だけでなく、サブチャネル間の相関も利用しているため良い推定精度を有する。しかしながら、例えば、式（３８）の状態ベクトルｘｐ（ｎ＋１）のサイズはＫ×Ｎであり、計算量が増大するといった問題がある。
【０１３６】
この問題を改良したものが、従来Per Subcarrier Kalmanと呼ばれている手法である。この方法は、上記方法をサブキャリア単位に分割して処理を行うものである。本発明の状態空間モデルを適用した場合、上記、式（４０）から式（４８）をサブキャリアｋ単位に分割することで所望の式が得られる。
【０１３７】
サブキャリア単位への分割は、状態遷移行列Ｃｐを例にとると、式（４２）において、Ｎ×Ｎの単位行列、零行列であるＩＮ×Ｎ、０Ｎ×Ｎを、スカラー量である１，０に置き換えることに相当する。したがって、行列のサイズは１／Ｎになり、計算量を低減することができる。無線伝送の場合には、音響分野とは異なり、推定される信号の品質によりも処理速度が要求される。そのため、本方法はより実用的である。
【０１３８】
次に、本実施の形態におけるチャネル推定精度の比較結果の一例について、図面を参照しつつ説明する。
【０１３９】
チャネル推定精度の比較検討は数値シミュレーションにより行った。数値シミュレーションの条件は次のように設定した。送信フレームは６４のシンボルからなり、そのうち４シンボルがパイロットシンボル、６０シンボルがデータシンボルである。総送信フレーム数は２００、総送信データシンボル数は１２０００である．評価量はＮＭＳＥ（Normarized Mean Square Error）を用いた。
【０１４０】
図１３は、本実施の形態におけるチャネル推定精度数値シミュレーションの一つの例の結果を示す図であり、ｆＤＴ＝０.０４５、（ＮＴｓ）＝０.０８μｓにおける、信号雑音比ＳＮＲに対するＮＭＳＥ特性を示している。ここでｆＤは最大ドップラー周波数、ＴはＯＦＤＭシンボル周期、Ｔｓはサンプリング間隔、τｍａｘは最大遅延スプレッドである。図１３から、本発明の方法の推定精度が従来手法に比べ向上していることがわか。これは、本発明の方法では、ＡＲ係数の推定誤差によるチャネルの推定精度の劣化が軽減できたためと考えられる。また、ｆＤが１８０Ｈｚより大きい場合、本発明の方法の性能向上がより大きなものとなることは明らかである。
【０１４１】
図１４は、本実施の形態におけるチャネル推定精度数値シミュレーションの他の例の結果を示す図であり、ＳＮＲ＝２０ｄＢ、τｍａｘ／（ＮＴｓ）＝０.０８μｓにおける、最大ドップラー周波数ｆＤに対するＮＭＳＥ特性を示している。この結果より、本発明の方法は、ＮＭＳＥが最大ドップラー周波数ｆＤに依存しないため、従来手法より良好な性能を有していることがわかる。すなわち、チャネル間干渉の影響が深刻となるフェージング変動の激しい環境においても、本発明の方法の有効性が確認できる。
【０１４２】
本発明は、上記各実施の形態に限定されるものではない。
【０１４３】
本発明の雑音抑圧装置は、ノイズが含まれた音声信号（得られた情報）からクリアな音声信号（必要な情報）を取り出すことが可能である。その一つの実施の形態として、カーナビゲーション装置などに必要不可欠な音声認識装置の前処理雑音除去装置が考えられる。
【０１４４】
また、画像分野においては、何らかの原因でぼけてしまったぼけ画像（得られた情報）からぼけのとれたクリアな画像（必要な情報）を取り出すことが可能であり、画像処理装置として活用可能である。
【０１４５】
さらに、従来、ＡＲ過程によるモデル化とカルマンフィルタアルゴリズムとを組み合わせを用いた通信・信号処理全般にわたり、本発明が適応可能であることはいうまでもない。
【０１４６】
また、医療分野では、従来、妊婦の胎児の状況を検査するには、個人が購入することができない高価な装置と高い専門知識とが必要があったが、本発明によれば、妊婦の体から母胎の心音や胎児の心音やその他の音（得られた情報）から不必要な音（情報）を取り除き胎児の心音（必要な情報）のみを取り出すことが可能になり、通院せずとも自宅で胎児の健康状態を、その心音から容易に確認することが可能となる。
【０１４７】
上記各実施の形態の説明に用いた各機能要素は、例えば、集積回路として実現される。これらは、個別に１チップ化されてもよいし、一部または全てを含むように１チップ化されてもよい。また、集積回路製造後にプログラムすることが可能なＦＰＧＡ（Ｆield Programmable Gate Array）や、回路を再構成可能なリコンフィギュラブル・プロセッサを利用してもよい。
【０１４８】
さらに、上記各実施の形態はハードウエアに限定されるものではなく、ソフトウエアによってもよい。その逆も真である。また、それらの組み合わせであってもよい。
【産業上の利用可能性】
【０１４９】
本発明に係る雑音抑圧装置および雑音抑圧方法は、カルマンフィルタを用いつつ、ＡＲ係数の推定を必要とすることなく、シンプルな構成で、雑音抑圧能力を向上することができる雑音抑圧装置および雑音抑圧方法として有用である。
【図面の簡単な説明】
【０１５０】
【図１】本発明の実施の形態１に係る雑音抑圧装置の構成を示すブロック図
【図２】本実施の形態に係る雑音抑圧装置の状態空間モデルを説明するブロック線図
【図３】本実施の形態に係る雑音抑圧装置のカルマンフィルタアルゴリズムに基づく信号推定手順を説明するフローチャート
【図４】図３のカルマンフィルタアルゴリズムの処理内容を示すフローチャート
【図５】本発明の実施の形態２に係る雑音抑圧装置の構成を示すブロック図
【図６】本実施の形態における音声波形数値シミュレーションの第１の例の結果を示す図
【図７】本実施の形態における音声波形数値シミュレーションの第２の例の結果を示す図
【図８】本実施の形態における音声波形数値シミュレーションの第３の例の結果を示す図
【図９】本実施の形態における音声波形数値シミュレーションの第４の例の結果を示す図
【図１０】本実施の形態における雑音抑圧後の音声品質の主観的評価結果の一例を示す図
【図１１】本実施の形態における雑音抑圧後の音声品質の主観的評価結果の他の例を示す図
【図１２】本発明の雑音抑圧方法が適用されたマルチキャリア受信装置の構成を示すブロック図
【図１３】本実施の形態におけるチャネル推定精度数値シミュレーションの一つの例の結果を示す図
【図１４】本実施の形態におけるチャネル推定精度数値シミュレーションの他の例の結果を示す図
【符号の説明】
【０１５１】
１００、２００    雑音抑圧装置
１１０    入力部
１２０    サンプリング部
１３０    Ａ／Ｄ変換部
１４０    バッファ
１５０    雑音抑圧処理部
１６０    出力部
２１０    パーソナルコンピュータ
２１１    操作装置
２１２    ディスプレイ
２１３    バスインタフェース
２１４    記録装置
２１５    主記憶メモリ
２１６    中央演算装置
２２０    マイクロホン
３００    マルチキャリア受信装置
３０１    検波部
３０２    ＧＩ除去部
３０３    チャネル等化部
３０４    チャネル推定部
３０５    ＦＦＴ部
３０６    復調部
３０７    復号部
５００、５０１    信号ベクトル
５０２    観測ベクトル
５０３    付加雑音信号ベクトル
５０４    駆動雑音ベクトル
５０５    状態遷移行列
５０６    状態遷移行列
【図１】
2008236270.tif 000072
【図２】
2008236270.tif 000073
【図３】
2008236270.tif 000074
【図４】
2008236270.tif 000075
【図５】
2008236270.tif 000076
【図６】
2008236270.tif 000077
【図７】
2008236270.tif 000078
【図８】
2008236270.tif 000079
【図９】
2008236270.tif 000080
【図１０】
2008236270.tif 000081
【図１１】
2008236270.tif 000082
【図１２】
2008236270.tif 000083
【図１３】
2008236270.tif 000084
【図１４】
2008236270.tif 000085