JP5115952B2 - 雑音抑圧装置および雑音抑圧方法 - Google Patents

Description

本発明は、カルマンフィルタに基づく雑音抑圧装置および雑音抑圧方法に関する。

所望の情報に不必要な情報が混在した受信情報（付加雑音などにより破損した情報）から不必要な情報を取り除き、所望情報のみを抽出することは、音声や無線通信、画像、姿勢制御などの分野における重要な技術であり、近年盛んに研究開発が行われている。

例えば、音声分野における公知の雑音抑圧方法としては、単一のマイクロホンを用いた方法や、複数のマイクロホンから構成されるマイクロホンアレイを用いた方法が提案されている。

しかしながら、マイクロホンアレイを用いた方法では、雑音信号の数が増大すると、マイクロホンの数も比例して増加することが避けらず、コストが増大する。そのため、現在は、単一のマイクロホンを用いた雑音抑圧方法の開発が主流となっている。

単一のマイクロホンしか用いない従来の雑音抑圧方法のアルゴリズムとしては、以下のようなものが知られている。

非特許文献１記載のＡＮＣ（適応ノイズキャンセラ）アルゴリズムは、音声信号の周期性を利用してノイズ信号を低減する。

非特許文献２には、線形予測に基づいた雑音抑圧アルゴリズムが記載されている。このアルゴリズムは、ピッチ推定や、雑音パワースペクトラム、雑音の平均方向に関する事前知識を必要としない。

また、上記アルゴリズムとは別に、カルマンフィルタに基づいた雑音抑圧アルゴリズムが、非特許文献３に提案されている。このアルゴリズムは、音声信号を自己回帰（ＡＲ：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

しかしながら、非特許文献１に記載された公知のアルゴリズムは、音声信号のピッチ周期の正確な推定を必要とする。そのため、このアルゴリズムは、その雑音抑圧能力が付加雑音によって劣化してしまうという問題点を有している。

この点、非特許文献２記載のアルゴリズムは、音声信号のピッチ周期の正確な推定を必要とせずに、雑音抑圧が可能である。さらに、このアルゴリズムは、その原理が単純であり、演算量を少なくすることができるといった長所を有している。しかし、このアルゴリズムは、その雑音抑圧能力が入力音声信号の周期性や線形性などの特性に依存している。言い換えると、このアルゴリズムは、雑音抑圧能力の中に音声信号に依存するパラメータが存在しているため、その実用には一定の限界がある。

非特許文献３記載のアルゴリズムは、強力な雑音抑圧能力を有し、特に高い音質を得たい音響分野への応用に適した手法である。

しかしながら、一方で、このアルゴリズムは、ＡＲ係数を必要とするため、ＡＲ係数の推定精度に雑音抑圧性能（つまり、当該カルマンフィルタアルゴリズムの性能）が大きく依存してしまうという問題点を有している。すなわち、ＡＲ係数が正確に推定されない場合、雑音を除去し切れないのみならず、場合によっては雑音に加えて音声信号まで除去してしまう可能性がある。これらは、雑音除去された音声信号の音質の劣化を引き起こす要因となりうる。

この点、一般には、ＡＲ係数の正確な推定は困難である。ＡＲ係数の正確な推定は、例えば、雑音除去であれば、クリアな音声信号に依存しているからである。このことは、音声信号が既知でなければならないことを意味しているため、リアルタイム処理は困難となる。また、仮に何らかの手法でリアルタイムにＡＲ係数を正確に推定することが可能となったとしても、処理が増加するため演算量の問題は避けられない。

本発明は、かかる点に鑑みてなされたものであり、ＡＲ係数の推定を必要としないシンプルで雑音抑圧能力が高い、カルマンフィルタに基づく雑音抑圧装置および雑音抑圧方法を提供することを目的とする。

本発明の雑音抑圧装置は、所望の情報に不必要な情報が混在した情報を取得する取得手段と、カルマンフィルタのみを用いて、取得された情報から前記不必要な情報を除去して前記所望情報を抽出する抽出手段と、を有し、前記カルマンフィルタは、状態空間モデルの観測方程式において自己回帰モデルの係数を使用しないように構成されている、構成を採る。

本発明の雑音抑圧方法は、所望の情報に不必要な情報が混在した情報を取得する取得ステップと、カルマンフィルタのみを用いて、取得した情報から前記不必要な情報を除去して前記所望情報を抽出する抽出ステップと、を有し、前記カルマンフィルタは、状態空間モデルの観測方程式において自己回帰モデルの係数を使用しないように構成されている、ようにした。

本発明によれば、カルマンフィルタを用いつつ、ＡＲ係数の推定を必要とすることなく、シンプルな構成で、雑音抑圧能力を向上することができる。

以下、本発明の実施の形態について、図面を参照して詳細に説明する。

（実施の形態１）
図１は、本発明の実施の形態１に係る雑音抑圧装置の構成を示すブロック図である。

図１に示す雑音抑圧装置１００は、入力部１１０、サンプリング部１２０、Ａ／Ｄ変換部１３０、バッファ１４０、雑音抑圧処理部１５０、および出力部１６０を有する。

入力部１１０は、観測信号を入力する。観測信号は、信号源からのクリアな信号と、付加雑音信号とが合わさった信号である。入力部１１０は、例えば、入力したアナログの観測信号を入力処理して、サンプリング部１２０に出力する。入力処理は、例えば、帯域制限処理や自動利得制御処理などである。

サンプリング部１２０は、所定のサンプリング周波数（例えば、１６ｋＨｚ）で、入力されたアナログの観測信号をサンプリング処理し、Ａ／Ｄ変換部１３０に出力する。Ａ／Ｄ変換部１３０は、サンプリングされた観測信号の振幅値を所定の分解能（例えば、８ｂｉｔ）でＡ／Ｄ変換処理し、バッファ１４０に送る。バッファ１４０は、所定のサンプリング数Ｎの信号フレーム（ブロック）を雑音抑圧処理部１５０に出力する。

雑音抑圧処理部１４０は、本発明の特徴的な構成要素であり、ＡＲ係数を用いないカルマンフィルタを内蔵している。すなわち、カルマンフィルタに基づく従来の雑音抑圧方法（非特許文献３参照）（以下「従来手法」という）では、線形予測を用いてＡＲ係数を推定した後、その結果を用いてカルマンフィルタを実行することで雑音抑圧を実現しているのに対し、本発明の雑音抑圧方法（以下「本手法」という）では、カルマンフィルタのみを用いて雑音抑圧を実現している。そのため、本手法では、信号源からのクリアな信号のみを用いて状態方程式を構成し、そのクリアな信号と付加雑音信号を用いて観測方程式を構成している。以下、雑音抑圧処理部１４０で行われるカルマンフィルタに基づく雑音抑圧処理動作について詳細に説明する。

雑音抑圧処理部１４０に入力される観測信号ｒ（ｎ）は、信号源からのクリアな信号（例えば、音声信号など）ｄ（ｎ）以外に付加雑音信号ｖ（ｎ）を含んでおり、次の式（１）を満たす。

ここで、ｎとは、装置の時刻ｎである。時刻ｎは、サンプリング部１２０で生成された離散的な時間系列において、処理開始時刻を時刻０と仮定したときに、そこからｎ番目の時刻のことを意味する。

従来のＡＲ過程に基づくモデル化方法（従来手法）では、信号ｄ（ｎ）は、ＡＲ係数を用いて次の式（２）でモデル化されると仮定している。

ここで、α_ｋ（ｎ）は時刻ｎでのＡＲ係数、ＫはＡＲ係数の次数、ｅ（ｎ）は信号ｄ（ｎ）が、式（２）に示すＫ次のＡＲ過程でモデル化されるとした場合の予測誤差（モデルリング誤差）である。

公知のように、従来手法では、付加雑音信号ｖ（ｎ）は、零平均であり、白色雑音であることが前提条件である。言い換えると、従来手法では、信号ｄ（ｎ）と付加雑音信号ｖ（ｎ）は無相関であり、つまり、次の式（３）を満たす。

ここで、Ｅ［・］はアンサンブル平均を示す。

ＡＲ係数は、ＡＲ係数推定アルゴリズムにより推定される。従来手法において最も重要な点は、カルマンフィルタを用いた高性能の雑音抑圧を達成するために、ＡＲ係数の正確な推定を必要とすることである。このことからも、カルマンフィルタの雑音抑圧能力がＡＲ係数の推定精度に大きく依存しているため雑音抑圧能力が大きく劣化することは容易に想像可能である。

本手法では、ＡＲ係数を用いずにカルマンフィルタの状態空間モデル、つまり音源信号からのクリアな信号のみを用いて状態方程式、およびそのクリアな信号と付加雑音信号とを用いて観測方程式を構成している。

以下、本手法での状態空間モデルの構成法について説明する。表記を容易にするために、まずＫ×１次の信号ベクトルｘ_ｐ（ｎ）を次の式（４）で定義する。添え字“ｐ”は本発明により考案された表現であることを示す。

次に、式（２）と同様にして信号ｄ（ｎ）のモデルを構成する。ここでは、ＡＲ係数を用いないとしても、従来手法と同様に、未知のＫ×Ｋ次の状態遷移行列Φ_ｐ（ｎ＋１）と、Ｋ×１次の駆動雑音ベクトルδ_ｐ（ｎ＋１）を導入する。これにより、本実施の形態の状態空間モデル（信号ベクトルにより記述される）の状態方程式の形として、次の式（５）が定まる。

式（５）が等式としてＡＲ係数を用いずに成立するためには、Ｋ×Ｋ次の状態遷移行列Φ_ｐ（ｎ＋１）、Ｋ×１次の駆動雑音ベクトルδ_ｐ（ｎ＋１）は、次の式（６）および式（７）を満たすことが求められる。

したがって、式（５）は次の式（８）に書き直せる。

式（６）で表されるＫ×Ｋ次の状態遷移行列Φ_ｐ（ｎ＋１）の行列要素は０と１のみであり、特に第１行がすべて０である。このことは、本手法の状態空間モデルの特徴の一つである。また、上記から明らかなように、式（５）、式（６）、式（７）の導出において、駆動雑音ベクトルδ_ｐ（ｎ＋１）に対して前提条件を与えていないことに注意すべきである。すなわち、駆動雑音ベクトルδ_ｐ（ｎ＋１）は有色であってよい。この理由については後述する。

本手法のカルマンフィルタアルゴリズムでは、特異値分解を避けるためにＫ×１次の観測ベクトルｙ_ｐ（ｎ＋１）を次の式（９）で定義する。

よって、式（１）および式（９）から、本手法の状態空間モデルの観測方程式として、次の式（１０）が導出される。

ここで、Ｍ_ｐはＫ×Ｋ次の状態遷移行列、ε_ｐはＫ×１次の付加雑音ベクトルであり、それぞれ次の式（１１）および式（１２）で定義される。

以上より、本手法のカルマンフィルタアルゴリズムで用いる状態空間モデルの状態方程式および観測方程式として、式（５）および式（１０）がそれぞれ求められた。

図２は、本実施の形態に係る雑音抑圧装置の状態空間モデルを説明するブロック線図である。

図２において、５００は時刻ｎにおける信号ベクトルｘ_ｐ（ｎ）、５０１は時刻ｎ＋１における信号ベクトルｘ_ｐ（ｎ＋１）、５０２は時刻ｎにおける観測ベクトルｙ_ｐ（ｎ）、５０３は時刻ｎにおける付加雑音ベクトルε_ｐ（ｎ）、５０４は時刻ｎ＋１における駆動雑音ベクトルδ_ｐ（ｎ＋１）、５０５は状態遷移行列Φ_ｐ、５０６は状態遷移行列Ｍ_ｐである。

図２から理解できるように、本手法のカルマンフィルタアルゴリズムは、駆動雑音が有色であるにもかかわらず、従来のカルマンフィルタアルゴリズムと同様の手順で実行できる（理由は後述する）。本手法のカルマンフィルタアルゴリズムに基づいた信号推定は、後述するフローチャートに従って行われる。

図３は、本実施の形態に係る雑音抑圧装置のカルマンフィルタアルゴリズムに基づく信号推定手順を説明するフローチャートである。

まず、推定信号ベクトルの初期値ｘ_ｅ（０｜０）、推定誤差の共分散行列の初期値Ｐ（０｜０）およびスカラー量である付加雑音ベクトルの共分散Ｒε_ｐ（ｎ）［ｉ，ｊ］の値を、次の式（１３）に示すように設定する（ＳＴ１００）。

ここで、Ｉは単位行列である。また、σ_ｖ ^２は付加雑音信号ｖ（ｎ）の雑音分散であり、既知と仮定している。もし付加雑音信号ｖ（ｎ）が白色雑音であり零平均であれば、σ_ｖ ^２は、次の式（１４）で与えられる。

ここで、Ｎは所定のサンプル数である。

次に、カルマンフィルタアルゴリズムを実行する（ＳＴ１０１）。このカルマンフィルタアルゴリズムの処理手順は、後で詳細に説明する。

次に、時刻ｎが所定のサンプル数Ｎに達したか否かを判定し（ＳＴ１０２）、所定のサンプル数Ｎに達していない場合は（ＳＴ１０２：ＮＯ）、ステップＳＴ１０１に戻って、カルマンフィルタアルゴリズムを繰り返し、所定のサンプル数Ｎに達した場合は（ＳＴ１０２：ＹＥＳ）、上記一連の処理を終了する。

本手法では、ＡＲ係数を用いずに状態空間モデルを設定している。そのため、従来手法で必要であったＡＲ係数を推定するステップを削減することができ、１段階処理で信号推定が可能となる。このことは本発明の大きな特徴の一つである。

図４は、図３のカルマンフィルタアルゴリズムの処理内容を示すフローチャートである。

ここで、開始時の時刻はｎ＋１であり、以下のフローでは、時刻ｎ＋１での観測信号ｒ（ｎ＋１）および観測ベクトルｙ_ｐ（ｎ＋１）と、時刻ｎでの推定誤差の共分散行列Ｐ（ｎ｜ｎ）とから、時刻ｎ＋１での推定信号ベクトルｘ_pe（ｎ＋１｜ｎ）を推定し、同時に推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ＋１）を更新する。

まず、次の式（１５）を用いて、駆動雑音ベクトルの共分散Ｒδ_ｐ（ｎ）［ｉ，ｊ］の値を計算する（ＳＴ２００）。

次に、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）を計算し、更新する（ＳＴ２０１）。ここで、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）の更新式は、次の式（１６）である。ここでは、以下、添え字“ｐｅ”および“ｅ”は、推定値であることを示す。

この式（１６）を展開すると、次の式（１７）となる。

ここで、式（１６）の第３項および第４項に注目する。第３項および第４項は、推定信号ベクトルｘ_ｐｅ（ｎ｜ｎ）と駆動雑音ベクトルδ_ｐ（ｎ｜ｎ）のアンサンブル平均を含んでいる、したがって、推定信号ベクトルｘ_ｐｅ（ｎ｜ｎ）と駆動雑音ベクトルδ_ｐ（ｎ｜ｎ）とが無相関でない、言い換えれば、駆動雑音ベクトルδ_ｐ（ｎ｜ｎ）が有色のとき、第３項および第４項は０にならない。すなわち、本手法のカルマンフィルタアルゴリズムが、駆動雑音信号が有色信号の場合に、従来と同様なカルマンフィルタアルゴリズム手順で適用可能か否かは、第３項および第４項が、式（１５）で表される推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）の更新にどの程度影響するかによる。

そこで、式（１７）の第３項に着目する。第３項は、次の式（１８）のように、

と書き表せる。

ここで、次の式（１９）のように、

とすれば、次の式（２０）となる。

式（５）および次の式（２１）

を考慮すれば、式（１８）は、次の式（２２）のように書き直すことができる。

ただし、次の式（２３）である。

ここで、ｄ_ｅ（ｎ）は時刻ｎでの推定信号である。

式（２３）の第１項は真値のアンサンブル平均、第２項は推定値のアンサンブル平均を示す。したがって、その差は小さく、ｅ（ｌ＋１）は無視できる値である。よって、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）の更新式（１５）は、次の式（２４）のように書き直される。このことから、駆動雑音が有色であっても、本手法の式（５）および式（９）から構成される状態空間モデル(状態方程式と観測方程式)であれば、従来手法と同様なカルマンフィルタアルゴリズム手順が適用可能となる。

次に、カルマンゲインＫ_ｐ（ｎ＋１）を計算する（ＳＴ２０２）。ここで、カルマンゲインＫ_ｐ（ｎ＋１）は、次の式（２５）で与えられる。

そして、これらの値から、まず推定信号ベクトルｘ_ｐｅ（ｎ＋１｜ｎ）を計算する（ＳＴ２０３）。ここで、推定信号ベクトルｘ_ｐｅ（ｎ＋１｜ｎ）は、次の式（２６）で与えられる。

そして、推定信号ベクトルｘ_ｐｅ（ｎ＋１｜ｎ＋１）を計算する（ＳＴ２０４）。ここで、推定信号ベクトルｘ_ｐｅ（ｎ＋１｜ｎ＋１）は、次の式（２７）で与えられる。

最後に、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ＋１）を計算、更新し（ＳＴ２０５）、フローを終了する。ここで、推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ＋１）は、次の式（２８）で与えられる。

雑音抑圧処理部１４０は、上述した手順、アルゴリズムによって推定された推定信号をを出力部１５０に出力する。出力部１５０は、例えば、スピーカやディスプレイ、通信手段、記憶装置などで構成されている。

ＡＲ係数の推定を必要とする従来のカルマンフィルタの問題点は、ＡＲ係数の推定精度にカルマンフィルタアルゴリズムの能力が依存していることである。これに対して、本手法は、カルマンフィルタアルゴリズムのみで実行しているため、ＡＲ係数の推定精度に依存しない。

なお、本手法のカルマンアルゴリズムは、以下に説明する更に改良された方法（以下「改良手法」という）によっても実行できる。以下、この改良手法について説明する。

改良手法のカルマンアルゴリズムの状態方程式は、上記式（５）と同一である。ただし、観測方程式を、次の式（２９）のように書き直す。改良手法の観測方程式（２９）は、本手法の観測方程式（１０）の一部の要素を取り出してきたものである。

ここで、ｙ’_ｐ（ｎ＋１）、ｍ^Ｔ _ｐ、ε’_ｐ（ｎ＋１）は、それぞれ改良手法における観測信号、雑音信号、状態遷移ベクトルであり、次の式（３０）を満たす。

改良手法の大きな特徴の一つは、本手法の観測方程式（１０）と改良手法の観測方程式（３０）の比較からわかるように、改良前の方法（本手法）における、観測ベクトルｙ_ｐ（ｎ＋１）、雑音ベクトルε_ｐ（ｎ＋１）、および状態遷移行列ｍ_ｐが、それぞれスカラーとベクトルに変更されていることである。したがって、改良手法はより少ない演算量で実行できるであろうことは明らかである。また、上記から明らかなように、改良手法においても、本手法と同様に、駆動雑音ベクトルδ_ｐ（ｎ＋１）に対して前提条件を与えていないことに注意すべきである。すなわち、改良手法においても駆動雑音ベクトルδ_ｐ（ｎ＋１）は有色であってよい。

すなわち、式（１７）で与えられた推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）を再度、行列形式で書き下すと、次の式（３１）となる。

ここで、また以下の式で、行列式中、斜線網掛け部分で示される部分は駆動雑音が有色であるときに影響を受ける部分である。

改良手法を用いたカルマンフィルタアルゴリズムは、図４に示す本手法のカルマンフィルタアルゴリズムとは、カルマンゲインＫ_ｐ（ｎ＋１）の計算（ＳＴ２０２）および推定信号ベクトルｘ_ｐｅ（ｎ＋１｜ｎ＋１）の計算（ＳＴ２０４）における計算式が異なる。

すなわち、カルマンゲインＫ_ｐ（ｎ＋１）の計算式は、次の式（３２）、つまり、

となる。ただし、雑音信号ε’_ｐ（ｎ＋１）の共分散をσ^２ _ｖとしたときに、次の式（３３）が成り立つ。

また、推定信号ベクトルｘ_pｅ（ｎ＋１｜ｎ＋１）の計算式は、次の式（３４）となる。

となる

ここで、式（３４）のｘ_pｅ（ｎ＋１｜ｎ＋１）の第１要素に着目すると、斜線網掛けではない（つまり、有色雑音の影響を受けない）ことから、駆動雑音が有色であってもクリアな信号の推定信号を得られる。したがって、改良手法も、本手法と同様に、従来のカルマンフィルタアルゴリズムで実行可能である。既に述べたように、改良手法では、本手法に比べて、演算量を低減できる。しかしながら、本手法が観測ベクトルｙ_ｐ（ｎ＋１）を用いているのに対して、改良手法方法では、スカラー量である観測信号ｙ‘_ｐ（ｎ＋１）を用いている。そのサイズの違いから、過去の観測量をより積極的に用いることができる本手法のほうが、より雑音抑圧能力が高いことは明らかである。

以上をまとめると、本発明における上記二つの手法（本手法および改良手法）は、カルマンフィルタに必要な状態空間モデルの観測方程式を変更することによって、演算量を大幅に減少させることが可能である。より具体的には、本発明における上記二つの手法では、ＡＲ係数の推定を必要としないため、その演算量は従来手法に比べて少ない。また、改良手法では、演算に用いる変数のサイズが本手法に比べて小さいため、その演算量は本手法に比べて少ない。すなわち、演算量に関して、
改良手法＜本手法＜従来手法
である。

また、本発明における上記二つの手法を、例えば、半導体集積回路や半導体回路などのハードウエアとして実施する場合、また、パーソナルコンピュータなどで実行可能なソフトウエアとして実施する場合のいずれにおいても、その構成は従来手法よりも単純化される。したがって、本発明における二つの手法を用いれば、回路規模やプログラム量が大幅に低減できるであろうことは明らかである。

また、本発明における上記二つの手法は、カルマンフィルタに必要な状態空間モデルの観測方程式を変更することによって、雑音抑圧能力を大幅に向上させることが可能である。より具体的には、本発明における上記二つの手法では、ＡＲ係数の推定を必要としないため、雑音抑圧能力がＡＲ係数の推定精度に依存する従来手法に比べて、雑音抑圧能力が高い。また、改良手法では、演算に用いる変数のサイズが本手法に比べて小さいため、過去の観測量を積極的に用いていない。したがって、その雑音抑圧能力は本手法に比べて低い。すなわち、雑音抑圧能力に関して、
従来手法＜改良手法＜本手法
である。

例えば、音響分野において、音質を若干落としても（ただし、従来手法よりは高い）、演算速度を速くしたい場合に、改良手法は有効である。

ＡＲ係数を用いない本発明における上記二つの手法は、状態方程式の駆動雑音が有色信号となっており、従来のカルマンフィルタの仮定、つまり、駆動雑音が白色信号であるという仮定に反している。しかしながら、本発明の状態空間モデルにおける状態方程式と観測方程式の構造的な性質により、駆動雑音が白色信号である従来のカルマンフィルタアルゴリズと全く同じアルゴリズムで実行可能であり、その理由は上記の通りである。また、本発明の改良手法においても、同様に実行可能であることは上記より明らかである。

本発明における上記二つの方法は、ＡＲ係数を使わないようにカルマンフィルタを構成する。言い換えると、カルマンフィルタのみで実行可能な雑音抑圧方法であるため、従来の雑音抑圧能力の問題を解決していることは明らかである。

（実施の形態２）
実施の形態２は、実施の形態１に示す雑音抑圧装置を音声に適用した場合である。

図５は、本発明の実施の形態２に係る雑音抑圧装置の構成を示すブロック図である。

図５に示す雑音抑圧装置２００は、本実施の形態の雑音抑圧処理を実行できるパーソナルコンピュータ２１０、マイクロホン２２０、サンプリング部１２０、およびＡ／Ｄ変換部１３０を有する。

パーソナルコンピュータ２１０は、操作装置２１１、ディスプレイ２１２、バスインタフェース２１３、記録装置２１４、主記憶メモリ２１５、および中央演算装置２１６を有する。操作装置２１１は、典型的にはキーボートやマウスなどであるが、音声認識装置などを用いてもよい。使用者は、操作装置２１１を用い、ディスプレイ２１２で確認をしながらコンピュータを操作できる。

パーソナルコンピュータ２１０において本実施の形態の雑音抑圧処理を実行させるプログラムソフトウエアは、記録装置２１４に格納されていてもよいし、バスインタフェース２１３を介して外部からダウンロードされてきてもよい。記録装置２１４は、典型的にはハードディスク装置であるが、ＣＤ−ＲＯＭ装置やＤＶＤ装置、フラッシュメモリなどの可搬性のあるものであってもよい。また、それらの組み合わせであってもよい。

サンプリング部１２０およびＡ／Ｄ変換部１３０は、パーソナルコンピュータ２１０内部に格納された内蔵カード（ボード）であってもよいし、バスインタフェース２１３を経由して接続された外部設置型機器であってもよい。

マイクロホン２２０からの観測音声信号は、サンプリング部１２０に入力される。サンプリング部１２０は、所定のサンプリング周波数（例えば、１６ｋＨｚ）で、入力されたアナログの観測音声信号をサンプリング処理し、Ａ／Ｄ変換部１３０に出力する。Ａ／Ｄ変換部１３０は、サンプリングされた観測音声信号の振幅値を所定の分解能（例えば、８ｂｉｔ）でＡ／Ｄ変換処理し、一時格納する。Ａ／Ｄ変換部１３０は、所定のサンプリング数Ｎの音声フレーム単位で、ディジタル化した観測音声信号をパーソナルコンピュータ２１０のバスインタフェース２１３に出力する。

パーソナルコンピュータ２１０はバスインタフェース２１３に出力された観測音声信号を一時、主記憶メモリ２１５に格納し、その後、所定の音声フレーム（サンプリング数）単位で、雑音抑圧処理を施した上で、再度主記憶メモリ２１５に格納する。雑音抑圧処理は、主記憶メモリ２１５や記録装置２１４に格納されたソフトウエアをバスインタフェース２１３経由で中央演算装置２１６に呼び出し、実行させることで行われる。

パーソナルコンピュータ２１０は、使用者の操作により、処理を実行したり、中断、終了させたりする。また使用者の操作により、処理した推定音声信号を外部に出力してもよい。

次に、本実施の形態における音声の雑音抑圧処理について、図面を参照しつつ説明する。

本実施の形態における音声の雑音抑圧の性能評価の目的で音声波形の数値シミュレーションを行った。

図６は、本実施の形態における音声波形シミュレーションの第１の例の結果を示す図であり、図７は、本実施の形態における音声波形シミュレーションの第２の例の結果を示す図である。また、図８は、本実施の形態における音声波形シミュレーションの第３の例の結果を示す図であり、図９は、本実施の形態における音声波形シミュレーションの第４の例の結果を示す図である。

日本人成人男性および女性のオリジナル音声信号（クリアな音声信号）は、無響室において、１６ｋＨｚのサンプリングレートでサンプリング、ディジタル化した。本例では、二つの音声信号サンプルを検討した。二つの音声信号サンプルは、（Ａ−１）図６（Ａ）に示す無声区間のないクリアな音声信号、（Ａ−２）図７（Ａ）に示す無声区間をもつクリアな音声信号である。それぞれの音声信号サンプルは、以後、音声（Ａ−１）および音声（Ａ−２）と参照する。

雑音信号は、人工的な付加雑音信号であり、本例では、二つの雑音信号サンプルを検討した。二つの雑音信号サンプルは、（Ｂ−１）図６（Ｂ）および図７（Ｂ）に示す付加白色ガウス雑音、（Ｂ−２）図８（Ｂ）および図９（Ｂ）に示す付加有色バブル雑音である。それぞれの雑音信号サンプルは、以後、雑音（Ｂ−１）および雑音（Ｂ−２）と参照される。

雑音信号の雑音分散σ_ｖ ^２は既知として次の式（３５）で表されるものとする。すなわち、

である。ここで、Ｌは全音声信号サンプル数である。

また、信号雑音比ＳＮＲ_ｉｎを式（３６）で定義する。

従来手法と本発明の方法（本手法）とによる雑音抑圧の結果を、音声（Ａ−１）および音声（Ａ−２）と雑音（Ｂ−１）の組み合わせ、ならびに音声（Ａ−１）および音声（Ａ−２）と雑音（Ｂ−２）の組み合わせからなる観測音声信号の信号波形から比較した。

図６（Ｃ）と図７（Ｃ）は、それぞれ音声（Ａ−１）と雑音（Ｂ−１）、音声（Ａ−２）と雑音（Ｂ−１）からなる観測音声信号波形である。

図６（Ｄ）と図７（Ｄ）は、それぞれ音声（Ａ−１）と雑音（Ｂ−１）、音声（Ａ−２）と雑音（Ｂ−１）の合成波形に対して従来手法で雑音抑圧を行った後の推定音声信号波形である。

一方、図６（Ｅ）と図７（Ｅ）は、それぞれ音声（Ａ−１）と雑音（Ｂ−１）、音声（Ａ−２）と雑音（Ｂ−１）の合成波形に対して従来手法で雑音抑圧を行った後の推定音声信号波形である。

図６（Ｄ）と図７（Ｄ）に対して図６（Ａ）と図７（Ａ）に示すオリジナルの音声信号とを参照すると、図６（Ｄ）と図７（Ｄ）に対して図６（Ａ）と図７（Ａ）に示すクリアな音声信号とを比較すると、従来手法による雑音抑圧では、雑音抑圧後に推定音声信号の振幅が小さくなっており、クリアな音声信号が抑圧されていることがわかる。加えて、従来手法による雑音抑圧では、サンプリング数の増加とともに、雑音抑圧後の推定音声信号の波形が、図６（Ａ）と図７（Ａ）に示すクリアな音声信号の波形から変形していく。

さらに、従来手法の雑音抑圧では、無声区間を有する音声（Ａ−２）に対して、推定音声信号が抑圧されるだけでなく、無声区間においてオリジナルの雑音信号と異なる雑音が観察されている。これは、従来手法では、無声区間では信号ｄ（ｎ）がゼロであるにもかかわらず、式（２）でＡＲ係数を求めようとするため、ＡＲ係数の値は発散し、不安定な状態を与えるからと推測される。

また、このことから、雑音信号が有色の場合、従来手法の適用は困難であろうことは容易に推測される。

従来手法とは対照的に、図６（Ｅ）および図７（Ｅ）に示す本発明の雑音抑圧では、雑音抑圧後の推定音声信号の波形が、図６（Ａ）と図７（Ａ）に示すクリアな音声信号の波形と非常に似通っている。

次に、図８（Ｄ）と図９（Ｄ）に対して図８（Ａ）と図９（Ａ）に示すクリアな音声信号とを比較すると、従来手法による雑音抑圧では、雑音（Ｂ−２）を含む観測音声信号に対して、非常に劣った結果を与えていることがわかる。これは、従来手法では、有色雑音である雑音（Ｂ−２）を含んだ観測音声信号に対してＡＲ係数を正確に推定することが困難であるためである。

一方、本発明の雑音抑圧方法では、雑音（Ｂ−２）の場合も雑音（Ｂ−１）の場合と同程度の雑音抑圧が達成されている。

すなわち、上述のように、本発明の雑音抑制方法は、白色、有色雑音、無声区間の有無に関わらず有効である。これは、本発明の雑音抑制方法の大きな特徴の一つである。

本発明の雑音抑制方法が有色の付加雑音にも適用できることは、既述のカルマンアルゴリズムにおける推定誤差の共分散行列Ｐ（ｎ＋１｜ｎ）の更新式（２３）の導出過程からも明らかであるが、上述の本実施の形態の数値シミュレーションは、このことを支持している。

次に、雑音抑制性能を数値的に比較するため、雑音抑圧性能ＳＮＲ_ｏｕｔを式（３７）で定義し、その数値シミュレーションを行った。

ここで、ｄ_ｅｉ（ｎ）は推定音声信号を表す。

表１は、本実施の形態における雑音抑圧性能の数値シミュレーションの一例の結果を示す表であり、音声（Ａ−１）と雑音（Ｂ−１）の組み合わせにおける、幾つかのＳＮＲ_ｉｎ、Ｋ（従来手法においては式（２）におけるＡＲ係数の次数、本発明では状態遷移行列Φ_ｐの次数）の値における、従来手法と本発明の方法の雑音抑圧性能ＳＮＲ_ｏｕｔを比較して示している。

表２は、本実施の形態における雑音抑圧性能の数値シミュレーションの他の例の結果を示す表であり、音声（Ａ−１）と雑音（Ｂ−１）の組み合わせにおける、幾つかのＳＮＲ_ｉｎ、Ｋの値における、従来手法と本発明の方法の雑音抑圧性能ＳＮＲ_ｏｕｔを比較して示している。

表１および表２を参照すると、本発明の方法は、すべてのＳＮＲ_ｉｎ、Ｋの値において、従来の方法に比べて雑音抑圧能力を改善していることがわかる。

特に、表２に示す有色雑音の場合には、従来手法は非常に劣った結果を与えているのにも関わらず、本発明の方法は表１に示す白色雑音の場合と同程度の結果を示している。すなわち、本発明の雑音抑圧方法は、白色雑音、有色雑音両者に効果的で雑音の性質に堅牢な雑音抑圧方法であるといえる。

また、表１および表２に見られるように、本発明の方法では、Ｋの値に対して雑音抑圧性能ＳＮＲ_ｏｕｔは安定であり、Ｋの値の増加に伴い増加する傾向にある。対照的に従来の方法では、表１に見られるように、Ｋの値に対して雑音抑圧性能ＳＮＲ_ｏｕｔは不安定である。これは、従来手法では、最適なＫの値、つまりＡＲ係数の次数を決定することが困難であることを意味している。

ＡＲ係数推定を必要とする従来手法において最も問題になることは、ＡＲ係数の次数の正確な推定は一般に困難であることである。なぜなら、ＡＲ係数の次数の正確な推定は，例えば雑音除去であれば、クリアな音声信号に依存しているからである。

このことは、クリアな音声信号が既知でなければならないことを意味しているため、リアルタイム処理は困難となる。ＡＲ係数の次数が正確でない場合には、カルマンフィルタアルゴリズの性能が劣化することは容易に想像可能である。また、何らかの手法でリアルタイムに推定することが可能となったとしても、処理が増加することより演算量などの問題を避けることは不可能である。

次に、本実施の形態に係る推定音声信号の音声品質を評価するためにリスニングテストによる主観的評価を行った。

音声品質評価に用いたオリジナル音声信号とオリジナル雑音信号は第２の実施の形態と同一であり、その説明は省略する。雑音信号は、異なるＳＮＲ_ｉｎ（０、５、および１０［ｄＢ］）で音声信号に加えた。

音声品質評価は、ＡＣＲ（絶対範疇評価）に基づいた５段階ＭＯＳ（平均オピニオン値）を用いたリスニングテストにより行った。５０人の聴取者が雑音抑圧により得られた推定音声信号のうち幾つかを評価した。各々の聴取者は、ポイント１からポイント５を決定する。ポイント５が最良である。

図１０は、本実施の形態における雑音抑圧後の音声品質の主観的評価結果の一つの例を示す図であり、音声（Ａ−１）と雑音（Ｂ−１）の組み合わせにおける、従来手法と本発明の方法とリスニングテストの結果を比較して示している。

また、図１１は、本実施の形態における雑音抑圧後の音声品質の主観的評価結果の他の例を示す図であり、音声（Ａ−１）と雑音（Ｂ−２）の組み合わせにおける、従来手法と本発明の方法とリスニングテストの結果を比較して示している。

図１０および図１１から、本発明の方法で推定した音声信号のスコアは、すべてのＳＮＲ_ｉｎ値において従来手法のスコアより高いことがわかる。特にその差は、音声（Ａ−１）と雑音（Ｂ−２）の組み合わせに対して大きい。

以上より、本発明の雑音抑圧方法は、音声信号の音声品質を犠牲にすることのない、白色雑音、有色雑音に効果的な優れた雑音抑圧方法であるといえる。

（実施の形態３）
本発明に係る雑音抑圧装置は、雑音抑圧以外の用途にも応用可能である。以下、その一つの応用例について説明する。

図１２は、本発明の雑音抑圧方法、つまりＡＲ係数を必要としないカルマンフィルタが適用されたマルチキャリア受信装置の構成を示すブロック図である。

図１２に示すマルチキャリア受信装置３００は、主に、検波部３０１、ＧＩ（ガードインターバル）除去部３０２、チャネル等化部３０３、チャネル推定部３０４、ＦＦＴ部３０５、復調部３０６、および復号部３０７を有する。

本実施の形態において、チャネル推定部３０４は、本発明の状態空間モデルを用いたカルマンフィルタに基づいたチャネル推定を実行できるように構成されている。より具体的には、チャネル推定部３０４は、サンプリング部１２０、Ａ／Ｄ変換部１３０、バッファ１４０、雑音抑圧処理部１５０、出力部１６０を有する構成をとる。

検波部３０１は伝送路上で周波数選択性フェージング等の影響を受けた受信信号を、中間周波数で直交検波し、ＧＩ除去部４０２に出力する。ＧＩ除去部３０２は、受信信号のガードインターバルを除去し、シンボル単位に連なった信号をチャネル等化部３０３と、チャネル推定部３０４に出力する。

チャネル推定部３０４に入力された信号は、入力部１１０に出力される。入力部は信号に所定の入力信号処理を施し、サンプリング部１２０に出力する。サンプリング部１２０は、所定のサンプリング周波数（例えば１６ｋＨｚ）で、入力されたアナログの受信信号をサンプリング処理し、Ａ／Ｄ変換部１３０に出力する。Ａ／Ｄ変換部１３０は、サンプリングされた受信信号の振幅値を所定の分解能（例えば８ｂｉｔ）でＡ／Ｄ変換処理し、バッファ１４０に送る。バッファ１４０は所定のサンプリング数Ｎの受信信号フレーム（ブロック）を雑音抑圧処理部１５０に出力する。

雑音抑圧処理部１５０は、カルマンフィルタに基づいたチャネル推定をサブキャリア全体、または各サブキャリアついて行い、推定された（サブ）チャネルゲインを出力部１６０に出力する。出力部１６０は、入力された推定チャネルゲインをチャネル等化部３０２に出力する。雑音処理部１５０での処理の詳細については後述する。

チャネル等化部３０３は、入力された推定チャネルゲインを用いて、ＧＩ除去部３０２から入力した信号の同期検波を行い、結果をＦＦＴ部３０５に出力する。ＦＦＴ部３０５は、フーリエ変換処理を行い、受信信号を各サブキャリア信号成分に分離して、復調部３０６に出力する。復調部３０６は信号の復調処理を行い、結果を復号部３０７に出力する。復号部３０７は、信号の復号を行い、デジタルデータを出力する。

本実施の形態では、まず、フェージングチャンネルを既述の状態空間モデルで表現する。マルチキャリア通信において、第ｋサブキャリアのみに着目すると、その受信信号ｙ_ｋ（ｔ）は、次の式（３８）

となる。ここでｎは時刻、Ｎは総キャリア数、Ｓ_ｋ（ｔ）、Ｈ_ｋ（ｔ）は第ｋキャリアの送信信号、チャネルゲインである。

表記を容易にするために、Ｎ×１次チャネルゲインベクトルｈ_p（ｎ）を次の式（３９）で定義する。

ここで、Ｋは後述の状態遷移行列の次数である。

本実施の形態では、本発明の状態空間モデルを用いてカルマンフィルタに基づき、チャネルゲインの推定値を適応的に求める。状態空間モデルの状態ベクトルｘ_ｐ（ｎ＋１）を次の式（４０）で定義する。

伝送路特性の時間変動を表す状態方程式は、次の式（４１）で記述される。

また、観測方程式は、次の式（４２）で記述される。

ここで、δ_ｐ（ｎ＋１）は駆動雑音ベクトル、ε_ｐ（ｎ＋１）は雑音信号ベクトル、Ｃ_ｐ、Ｄ_ｐは状態遷移行列、Ｇは定数行列である。

駆動雑音ベクトルδ_ｐ（ｎ＋１）、雑音信号ベクトルε_ｐ（ｎ＋１）、状態遷移行列Ｃ_ｐ、Ｄ_ｐ、定数行列Ｇは以下の式（４３）から式（４８）で定義される。

ここで、Ｉ_Ｎ×Ｎ、０_Ｎ×ＮはＮ×Ｎの単位行列、零行列である。

上記の定式化は、従来Vector Kermanと呼ばれているチャネルゲイン推定法に、本発明の状態空間モデルを適用したものである。フェージングチャネル間の相関だけでなく、サブチャネル間の相関も利用しているため良い推定精度を有する。しかしながら、例えば、式（３８）の状態ベクトルｘ_ｐ（ｎ＋１）のサイズはＫ×Ｎであり、計算量が増大するといった問題がある。

この問題を改良したものが、従来Per Subcarrier Kalmanと呼ばれている手法である。この方法は、上記方法をサブキャリア単位に分割して処理を行うものである。本発明の状態空間モデルを適用した場合、上記、式（４０）から式（４８）をサブキャリアｋ単位に分割することで所望の式が得られる。

サブキャリア単位への分割は、状態遷移行列Ｃ_ｐを例にとると、式（４２）において、Ｎ×Ｎの単位行列、零行列であるＩ_Ｎ×Ｎ、０_Ｎ×Ｎを、スカラー量である１，０に置き換えることに相当する。したがって、行列のサイズは１／Ｎになり、計算量を低減することができる。無線伝送の場合には、音響分野とは異なり、推定される信号の品質によりも処理速度が要求される。そのため、本方法はより実用的である。

次に、本実施の形態におけるチャネル推定精度の比較結果の一例について、図面を参照しつつ説明する。

チャネル推定精度の比較検討は数値シミュレーションにより行った。数値シミュレーションの条件は次のように設定した。送信フレームは６４のシンボルからなり、そのうち４シンボルがパイロットシンボル、６０シンボルがデータシンボルである。総送信フレーム数は２００、総送信データシンボル数は１２０００である．評価量はＮＭＳＥ（Normarized Mean Square Error）を用いた。

図１３は、本実施の形態におけるチャネル推定精度数値シミュレーションの一つの例の結果を示す図であり、ｆ_ＤＴ＝０.０４５、（ＮＴ_ｓ）＝０.０８μｓにおける、信号雑音比ＳＮＲに対するＮＭＳＥ特性を示している。ここでｆ_Ｄは最大ドップラー周波数、ＴはＯＦＤＭシンボル周期、Ｔ_ｓはサンプリング間隔、τ_ｍａｘは最大遅延スプレッドである。図１３から、本発明の方法の推定精度が従来手法に比べ向上していることがわか。これは、本発明の方法では、ＡＲ係数の推定誤差によるチャネルの推定精度の劣化が軽減できたためと考えられる。また、ｆ_Ｄが１８０Ｈｚより大きい場合、本発明の方法の性能向上がより大きなものとなることは明らかである。

図１４は、本実施の形態におけるチャネル推定精度数値シミュレーションの他の例の結果を示す図であり、ＳＮＲ＝２０ｄＢ、τ_ｍａｘ／（ＮＴ_ｓ）＝０.０８μｓにおける、最大ドップラー周波数ｆ_Ｄに対するＮＭＳＥ特性を示している。この結果より、本発明の方法は、ＮＭＳＥが最大ドップラー周波数ｆ_Ｄに依存しないため、従来手法より良好な性能を有していることがわかる。すなわち、チャネル間干渉の影響が深刻となるフェージング変動の激しい環境においても、本発明の方法の有効性が確認できる。

本発明は、上記各実施の形態に限定されるものではない。

本発明の雑音抑圧装置は、ノイズが含まれた音声信号（得られた情報）からクリアな音声信号（必要な情報）を取り出すことが可能である。その一つの実施の形態として、カーナビゲーション装置などに必要不可欠な音声認識装置の前処理雑音除去装置が考えられる。

また、画像分野においては、何らかの原因でぼけてしまったぼけ画像（得られた情報）からぼけのとれたクリアな画像（必要な情報）を取り出すことが可能であり、画像処理装置として活用可能である。

さらに、従来、ＡＲ過程によるモデル化とカルマンフィルタアルゴリズムとを組み合わせを用いた通信・信号処理全般にわたり、本発明が適応可能であることはいうまでもない。

また、医療分野では、従来、妊婦の胎児の状況を検査するには、個人が購入することができない高価な装置と高い専門知識とが必要があったが、本発明によれば、妊婦の体から母胎の心音や胎児の心音やその他の音（得られた情報）から不必要な音（情報）を取り除き胎児の心音（必要な情報）のみを取り出すことが可能になり、通院せずとも自宅で胎児の健康状態を、その心音から容易に確認することが可能となる。

上記各実施の形態の説明に用いた各機能要素は、例えば、集積回路として実現される。これらは、個別に１チップ化されてもよいし、一部または全てを含むように１チップ化されてもよい。また、集積回路製造後にプログラムすることが可能なＦＰＧＡ（Ｆield Programmable Gate Array）や、回路を再構成可能なリコンフィギュラブル・プロセッサを利用してもよい。

さらに、上記各実施の形態はハードウエアに限定されるものではなく、ソフトウエアによってもよい。その逆も真である。また、それらの組み合わせであってもよい。

本発明に係る雑音抑圧装置および雑音抑圧方法は、カルマンフィルタを用いつつ、ＡＲ係数の推定を必要とすることなく、シンプルな構成で、雑音抑圧能力を向上することができる雑音抑圧装置および雑音抑圧方法として有用である。

本発明の実施の形態１に係る雑音抑圧装置の構成を示すブロック図 本実施の形態に係る雑音抑圧装置の状態空間モデルを説明するブロック線図 本実施の形態に係る雑音抑圧装置のカルマンフィルタアルゴリズムに基づく信号推定手順を説明するフローチャート 図３のカルマンフィルタアルゴリズムの処理内容を示すフローチャート 本発明の実施の形態２に係る雑音抑圧装置の構成を示すブロック図 本実施の形態における音声波形数値シミュレーションの第１の例の結果を示す図 本実施の形態における音声波形数値シミュレーションの第２の例の結果を示す図 本実施の形態における音声波形数値シミュレーションの第３の例の結果を示す図 本実施の形態における音声波形数値シミュレーションの第４の例の結果を示す図 本実施の形態における雑音抑圧後の音声品質の主観的評価結果の一例を示す図 本実施の形態における雑音抑圧後の音声品質の主観的評価結果の他の例を示す図 本発明の雑音抑圧方法が適用されたマルチキャリア受信装置の構成を示すブロック図 本実施の形態におけるチャネル推定精度数値シミュレーションの一つの例の結果を示す図 本実施の形態におけるチャネル推定精度数値シミュレーションの他の例の結果を示す図

符号の説明

１００、２００ 雑音抑圧装置
１１０ 入力部
１２０ サンプリング部
１３０ Ａ／Ｄ変換部
１４０ バッファ
１５０ 雑音抑圧処理部
１６０ 出力部
２１０ パーソナルコンピュータ
２１１ 操作装置
２１２ ディスプレイ
２１３ バスインタフェース
２１４ 記録装置
２１５ 主記憶メモリ
２１６ 中央演算装置
２２０ マイクロホン
３００ マルチキャリア受信装置
３０１ 検波部
３０２ ＧＩ除去部
３０３ チャネル等化部
３０４ チャネル推定部
３０５ ＦＦＴ部
３０６ 復調部
３０７ 復号部
５００、５０１ 信号ベクトル
５０２ 観測ベクトル
５０３ 付加雑音信号ベクトル
５０４ 駆動雑音ベクトル
５０５ 状態遷移行列
５０６ 状態遷移行列

Claims (9)

1. 所望の情報に雑音が混在した観測情報のみから前記所望情報を推定する雑音抑圧装置であって、
   前記観測情報を取得する取得手段と、
   有色雑音を駆動源として用いるカルマンフィルタを用いて、取得された前記観測情報から前記雑音を除去して前記所望情報を抽出する抽出手段と、を有し、
   前記カルマンフィルタは、
   状態空間モデルの状態方程式において自己回帰モデルの係数を使用しないように構成されている、
   雑音抑制装置。
2. 前記抽出手段は、
   時刻ｎのみの観測情報に対して、時刻ｎまでの情報により前記所望情報を含む時刻ｎ＋１のシステムの状態量を推定した場合の推定誤差の第１相関値行列を算出する第１の相関演算部と、
   時刻ｎのみの観測情報に対して、前記第１の相関演算部によって算出された前記推定誤差の第１相関値行列を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値ベクトルと、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値ベクトルと、前記観測情報を含む観測量の推定誤差ベクトルと、の関係を規定するための重み係数行列を算出する重み係数算出部と、
   時刻ｎのみの観測情報に対して、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の第１最適推定値ベクトルを算出する第１の最適推定値算出部と、
   時刻ｎのみの観測情報に対して、前記重み係数算出部によって算出された前記重み係数行列を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の第２最適推定値ベクトルを算出する第２の最適推定値算出部と、
   時刻ｎのみの観測情報に対して、時刻ｎ＋１までの情報により当該時刻の前記状態量を推定した場合の推定誤差の第２相関値行列を算出する第２の相関演算部と、
   を有する請求項１記載の雑音抑圧装置。
3. 前記第１の相関演算部は、
   所定の状態遷移行列、与えられた駆動源ベクトルの共分散の要素値、および与えられたまたは前回前記第２の相関演算部によって算出された前記推定誤差の第２相関値行列を用いて、前記推定誤差の第１相関値行列の算出を行い、
   前記重み係数算出部は、
   前記第１の相関演算部によって算出された前記推定誤差の第１相関値行列、与えられた観測遷移行列、および与えられた雑音ベクトルの共分散の要素値を用いて、前記重み係数行列の算出を行い、
   前記第１の最適推定値算出部は、
   前記状態遷移行列、および、与えられたまたは前回前記第２の最適推定値算出部によって算出された前記状態量の第２最適推定値ベクトルを用いて、前記状態量の第１最適推定値ベクトルの算出を行い、
   前記第２の最適推定値算出部は、
   前記第１の最適推定値算出部によって算出された前記状態量の第１最適推定値ベクトル、前記重み係数算出部によって算出された前記重み係数行列、前記観測遷移行列、および時刻ｎ＋１のみにおける観測量を用いて、前記状態量の第２最適推定値ベクトルの算出を行い、
   前記第２の相関演算部は、
   前記重み係数算出部によって算出された前記重み係数行列、前記観測遷移行列、および前記第１の相関演算部によって算出された前記推定誤差の第１相関値行列を用いて、前記推定誤差の第２相関値行列の算出を行う、
   請求項２記載の雑音抑圧装置。
4. 所望の情報に雑音が混在した観測情報のみから前記所望情報を推定する雑音抑圧方法であって、
   前記観測情報を取得する取得ステップと、
   有色雑音を駆動源として用いるカルマンフィルタを用いて、取得した前記観測情報から前記雑音を除去して前記所望情報を抽出する抽出ステップと、を有し、
   前記カルマンフィルタは、
   状態空間モデルの状態方程式において自己回帰モデルの係数を使用しないように構成されている、
   雑音抑圧方法。
5. 前記抽出ステップは、
   時刻ｎのみの観測情報に対して、時刻ｎまでの情報により前記所望情報を含む時刻ｎ＋１のシステムの状態量を推定した場合の推定誤差の第１相関値行列を算出する第１の相関演算工程と、
   時刻ｎのみの観測情報に対して、前記第１の相関演算工程で算出した前記推定誤差の第１相関値行列を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値ベクトルと、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値ベクトルと、前記観測情報を含む観測量の推定誤差ベクトルと、の関係を規定するための重み係数行列を算出する重み係数算出工程と、
   時刻ｎのみの観測情報に対して、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の第１最適推定値ベクトルを算出する第１の最適推定値算出工程と、
   時刻ｎのみの観測情報に対して、前記重み係数算出工程で算出した前記重み係数行列を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の第２最適推定値ベクトルを算出する第２の最適推定値算出工程と、
   時刻ｎのみの観測情報に対して、時刻ｎ＋１までの情報により当該時刻の前記状態量を推定した場合の推定誤差の第２相関値行列を算出する第２の相関演算工程と、
   を有する請求項４記載の雑音抑圧方法。
6. 前記第１の相関演算工程は、
   所定の状態遷移行列、与えられた駆動源ベクトルの共分散の要素値、および与えられたまたは前回前記第２の相関演算工程で算出した前記推定誤差の第２相関値行列を用いて、前記推定誤差の第１相関値行列の算出を行い、
   前記重み係数算出工程は、
   前記第１の相関演算工程で算出した前記推定誤差の第１相関値行列、与えられた観測遷移行列、および与えられた雑音ベクトルの共分散の要素値を用いて、前記重み係数行列の算出を行い、
   前記第１の最適推定値算出工程は、
   前記状態遷移行列、および、与えられたまたは前回前記第２の最適推定値算出工程で算出した前記状態量の第２最適推定値ベクトルを用いて、前記状態量の第１最適推定値ベクトルの算出を行い、
   前記第２の最適推定値算出工程は、
   前記第１の最適推定値算出工程で算出した前記状態量の第１最適推定値ベクトル、前記重み係数算出工程で算出した前記重み係数行列、前記観測遷移行列、および時刻ｎ＋１のみにおける観測量を用いて、前記状態量の第２最適推定値ベクトルの算出を行い、
   前記第２の相関演算工程は、
   前記重み係数算出工程で算出した前記重み係数行列、前記観測遷移行列、および前記第１の相関演算工程で算出した前記推定誤差の第１相関値行列を用いて、前記推定誤差の第２相関値行列の算出を行う、
   請求項５記載の雑音抑圧方法。
7. 所望の情報に雑音が混在した観測情報のみから前記所望情報を推定するための雑音抑圧プログラムであって、
   コンピュータに、
   前記観測情報を取得する取得ステップと、
   有色雑音を駆動源として用いるカルマンフィルタを用いて、取得した前記観測情報から前記雑音を除去して前記所望情報を抽出する抽出ステップと、前記カルマンフィルタは、状態空間モデルの状態方程式において自己回帰モデルの係数を使用しないように構成されている、
   を実行させるための雑音抑制プログラム。
8. 前記抽出ステップは、
   時刻ｎのみの観測情報に対して、時刻ｎまでの情報により前記所望情報を含む時刻ｎ＋１のシステムの状態量を推定した場合の推定誤差の第１相関値行列を算出する第１の相関演算工程と、
   時刻ｎのみの観測情報に対して、前記第１の相関演算工程で算出した前記推定誤差の第１相関値行列を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の最適推定値ベクトルと、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の最適推定値ベクトルと、前記観測情報を含む観測量の推定誤差ベクトルと、の関係を規定するための重み係数行列を算出する重み係数算出工程と、
   時刻ｎのみの観測情報に対して、時刻ｎまでの情報による時刻ｎ＋１での前記状態量の第１最適推定値ベクトルを算出する第１の最適推定値算出工程と、
   時刻ｎのみの観測情報に対して、前記重み係数算出工程で算出した前記重み係数行列を用いて、時刻ｎ＋１までの情報による当該時刻での前記状態量の第２最適推定値ベクトルを算出する第２の最適推定値算出工程と、
   時刻ｎのみの観測情報に対して、時刻ｎ＋１までの情報により当該時刻の前記状態量を推定した場合の推定誤差の第２相関値行列を算出する第２の相関演算工程と、
   を有する請求項７記載の雑音抑圧プログラム。
9. 前記第１の相関演算工程は、
   所定の状態遷移行列、与えられた駆動源ベクトルの共分散の要素値、および与えられたまたは前回前記第２の相関演算工程で算出した前記推定誤差の第２相関値行列を用いて、前記推定誤差の第１相関値行列の算出を行い、
   前記重み係数算出工程は、
   前記第１の相関演算工程で算出した前記推定誤差の第１相関値行列、与えられた観測遷移行列、および与えられた雑音ベクトルの共分散の要素値を用いて、前記重み係数行列の算出を行い、
   前記第１の最適推定値算出工程は、
   前記状態遷移行列、および、与えられたまたは前回前記第２の最適推定値算出工程で算出した前記状態量の第２最適推定値ベクトルを用いて、前記状態量の第１最適推定値ベクトルの算出を行い、
   前記第２の最適推定値算出工程は、
   前記第１の最適推定値算出工程で算出した前記状態量の第１最適推定値ベクトル、前記重み係数算出工程で算出した前記重み係数行列、前記観測遷移行列、および時刻ｎ＋１のみにおける観測量を用いて、前記状態量の第２最適推定値ベクトルの算出を行い、
   前記第２の相関演算工程は、
   前記重み係数算出工程で算出した前記重み係数行列、前記観測遷移行列、および前記第１の相関演算工程で算出した前記推定誤差の第１相関値行列を用いて、前記推定誤差の第２相関値行列の算出を行う、
   請求項８記載の雑音抑圧プログラム。

Images