JP6950490B2 - フィルタリング装置及びフィルタリング装置のテーブル作成方法 - Google Patents

Description

本発明は、複数のセンサからの出力に対し整相処理するフィルタリング装置及びフィルタリング装置のテーブル作成方法に関する。

従来、アレイを構成する複数のセンサに対して、レベルの補正（シェーディング）と遅延時間の付与とを行い、それらの結果を加算する整相処理が知られている。レベルの補正と遅延時間の付与とは、センサ出力とフィルタ係数との畳み込み演算で実現される。フィルタ係数とは、レベルの補正と付与する遅延時間とに対応する時間領域の係数である。アレイをＭ個のセンサで構成した場合の整相処理の実現例を図１２に示す。図１２の整相処理は（１）式で表される。

・・・（１）

（１）式において、ｘｍ（ｎ）はｍ番目のセンサ出力であり、ｗｍ（ｉ）はｍ番目のセンサに与えるｉ番目のフィルタ係数であり、それぞれ各センサでＮ個の系列となっている。ただし、ｗｍ^＊（ｉ）にあるような、上付きの＊は複素共役を意味している。ｙ（ｎ）は時刻ｎにおける整相出力であり、Ｔｓはサンプリング周期である。しかしながら、整相処理では、各センサで所望の特性（方位特性、周波数特性）を得るために、フィルタ係数の長さＮは、数十〜数百と長くなることが多くなる。このような場合、ＦＦＴ（Ｆａｓｔ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）を用いて、（1）式のような時間領域での畳み込み演算を、周波数領域での積の演算によって高速化を図ることが一般的である。（１）式の整相出力を周波数領域で実現するために(２)式のように書き換える。

・・・（２）

・・・（３）

・・・（４）

・・・（５）

時間領域での畳み込み演算を周波数領域での積の演算に変換することによって、処理の高速化を図ることができるものの、ＦＦＴでは入力されるセンサ出力が周期的であると仮定しているという問題がある。この仮定によって、ＦＦＴに入力されるセンサ出力は実際には有限長の非周期的な波形であるため、センサ出力の両端が不連続であることの影響が、整相処理出力にも現れてしまう。そこで、非周期的な信号にも対応できる実現方法として、図１３に示すＯｖｅｒｌａｐ−ｓａｖｅ法を用いる方法が知られている。

図１３は、Ｏｖｅｒｌａｐ−ｓａｖｅ法を用いた整相処理の処理ブロック図である。Ｏｖｅｒｌａｐ−ｓａｖｅ法では、ｍ番目のセンサのフィルタ係数ｗｍ（ｎ）の後半にゼロ詰めをしてＦＦＴを実施したＷｍ（ｋ）を、入力信号ｘｍ（ｎ）のＦＦＴ結果Ｘｍ（ｋ）と掛け合わせる。この処理を全センサに対して実施し、各センサの処理結果の加算を行うことで、ｋ番目の周波数ビンの整相出力Ｙ（ｋ）を得る。各周波数ビンに同様の処理を実施した後、ＩＦＦＴ（Ｉｎｖｅｒｓｅ ＦＦＴ）を実施し、波形に戻した後に後半部分（ゼロ詰めした部分）のみを抽出する。

また、図１３のフィルタ係数計算部に示すように、フィルタ係数を作る際は、整相器としての所望の方位特性を実現するように、各センサでフィルタ周波数特性を作成しておき、それをＩＦＦＴしてゼロ詰め前のフィルタ係数を作成する必要がある。ここで、フィルタ周波数特性とは、フィルタ係数の周波数領域での表現である。整相器の特性が時間変化しない場合においては、図１３のフィルタ周波数特性の作成処理は、整相処理を実施する前に、事前に実施すればよい。

John Shynk,’Frequency-Domain and Multirate Adaptive Filtering’, 1992

図１３に示した方法はＦＦＴを用いた高速化手法が実現できるが、図１３における構成ではフィルタ係数が時間変化するような場合、例えば従来整相処理において時々刻々音速が変化する場合、もしくは後述する適応整相処理の場合は、図１３のフィルタ周波数特性の作成処理をリアルタイムで計算する必要がある。特にフィルタ係数作成部のＩＦＦＴ演算及びフィルタ係数にゼロ詰めした後のＦＦＴ演算の処理は演算量が多い。非特許文献１の図５にある適応フィルタの例でも、毎周期変化するフィルタ周波数特性に対して、ＩＦＦＴ演算、ゼロ詰め及びＦＦＴを毎周期実施している。

さらに、整相処理では一般的に複数（数十〜数百）の方位を同時に整相することが多いため、ＦＦＴ／ＩＦＦＴ演算がセンサ数×整相方位数だけ毎周期発生し、処理量が膨大になる。このように、従来のＯｖｅｒｌａｐ−ｓａｖｅ法を用いた整相処理は、フィルタ係数が時間変化である場合にＦＦＴ/ゼロ詰め／ＩＦＦＴ処理が毎周期発生し、計算量が増大する問題がある。

本発明は、上記のような課題を背景としてなされたもので、フィルタ周波数特性の時間変化部分の一部を事前にテーブル化し、整相処理を実施する際にはそのテーブルを参照の参照と後述する位相回転を実施することにより、従来技術では毎周期必要であったフィルタ係数のＩＦＦＴ及びＦＦＴ処理を削減し高速化することができるフィルタリング装置及びフィルタリング装置のテーブル作成方法を提供するものである。

本発明に係るフィルタリング装置は、複数のセンサからの出力に対するフィルタリング装置であって、フィルタの位相特性に与える遅延時間として、遅延時間の範囲を複数の区間に分割し、分割された各区間において遅延時間に対応する位相特性を求め、所定の振幅を有するバンドパスフィルタを周波数領域上で構築する構築手段と、構築手段によって構築されたバンドパスフィルタの周波数特性に基づいて時間領域上のインパルス応答を求める変換手段と、変換手段によって求められた時間領域上のインパルス応答の主要部を抽出する抽出手段と、抽出手段によって抽出されたインパルス応答の主要部の両端にゼロデータを挿入する挿入手段と、挿入手段によってゼロデータが挿入されたインパルス応答に基づいて得られる周波数特性をテーブルとするテーブル化手段と、を有するテーブル作成部と、周波数特性のテーブルに基づいて、フィルタの位相特性に与える遅延時間のうちサンプリング周波数以下の詳細遅延時間に対応する周波数特性を取得し、サンプリング周波数の整数倍のサンプル遅延時間に位相回転を実施して最終フィルタ周波数特性を求める特性取得部と、を備えることを特徴とする。

本発明によれば、フィルタ周波数特性の時間変化部分の一部を事前にテーブル化し、整相処理を実施する際にはそのテーブルを参照して位相回転を実施することにより、従来技術では毎周期必要であったフィルタ係数のＩＦＦＴ及びＦＦＴ処理を削減し高速化することができる。

本実施の形態１に係るフィルタリング装置１を示すブロック図である。 座標系を説明する図である。 本実施の形態１に係るテーブル作成部２３を示すブロック図である。 本実施の形態１に係るテーブルを示す模式図である。 本実施の形態１に係るフィルタリング装置１の処理結果を示すグラフである。 本実施の形態２に係るフィルタリング装置１００を示すブロック図である。 本実施の形態２に係るテーブル作成部１２４を示すブロック図である。 本実施の形態２に係るテーブルを示す模式図である。 本実施の形態２に係る平均部１２２及び補正量計算部１２３の機能を説明する図である。 本実施の形態２に係る加算手段１４３の機能を説明する図である。 本実施の形態２に係るフィルタリング装置１００の処理結果を示すグラフである。 従来の整相処理を説明する図である。 従来のフィルタリング装置を示すブロック図である。

実施の形態１．
以下、本発明に係るフィルタリング装置及びフィルタリング装置のテーブル作成方法の実施の形態について、図面を参照しながら説明する。図１は、本実施の形態１に係るフィルタリング装置１を示すブロック図である。図１に示すように、フィルタリング装置１は、フィルタ処理部２と、フィルタ特性出力部３とを有している。フィルタ処理部２は、高速フーリエ変換を実施するＦＦＴ１１、逆高速フーリエ変換を実施するＩＦＦＴ１２及び切り出し１３の機能を有している。フィルタ特性出力部３は、遅延計算部２１と、分解部２２と、テーブル作成部２３と、特性取得部２４とを備えている。

（遅延計算部２１）
遅延計算部２１は、整相方位（θ,φ）、ｍ番目センサのセンサ位置ｒｍ、音速ｃを入力として、ｍ番目センサ出力に与える遅延時間τｍを計算し、後段の分解部２２に出力する。図１において、アレイを構成するセンサの位置及び音速は、別の手段によって計測されているとする。ある整相方位（θ,φ）に対するｍ番目のセンサ出力に与える遅延時間τｍを、次の式で計算する。なお、座標系を図２に示す。

・・・（６）

ただし、ｒｍはｍ番目のセンサの位置ベクトル、ｂ（θ,φ）は整相方位（θ,φ）における方向ベクトル、ｃは音速である。なお、ｒｍＴにあるような、上付きのＴはベクトル又は行列の転置を意味するものとする。

・・・（７）

・・・（８）

（分解部２２）
分解部２２は、ｍ番目センサ出力に与える遅延時間τｍとサンプリング周波数ｆｓとを入力とし、サンプル遅延数（サンプル遅延時間）ｑｍと詳細遅延（詳細遅延時間）δｔｍとを出力する。分解部２２は、(６)式で計算した、ｍ番目のセンサ出力に与える遅延時間τｍとサンプリング周波数ｆｓを使用し、サンプル遅延数ｑｍと詳細遅延δｔｍとに分解する。

・・・（９）

・・・（１０）

ここで、サンプル遅延数ｑｍはサンプリング周波数の逆数の整数倍で表現できる遅延であり、詳細遅延δｔｍはサンプリング周期以下の短い遅延時間である。ｒｏｕｎｄ（Ａ）はＡを小数点第一位で四捨五入する処理である。

（テーブル作成部２３）
テーブル作成部２３は、詳細遅延δｔｍを入力とし、この詳細遅延δｔｍに対応するフィルタ周波数特性ｔｂｌｍ（ｋ）を、位相回転手段４１において位相回転して出力する。ここでｋは周波数ビン番号であり、ｋ＝０，１，・・・，Ｎ−１の複数の周波数ビンに対応するフィルタ特性が出力される。フィルタ周波数特性テーブルは、詳細遅延に対応するフィルタ周波数特性を格納したテーブルであり、詳細遅延δｔｍは、式（９）及び（１０）の定義から、−１／２ｆｓ〜１／２ｆｓの範囲の値を取る。

ここで、図３を用いて、フィルタ周波数特性テーブルの作成方法を述べる。テーブル作成部２３は、構築手段３１と、変換手段３２と、抽出手段３３と、挿入手段３４と、ＦＦＴ３５と、テーブル化手段３６とを有している。構築手段３１は、フィルタ周波数特性の振幅特性と位相特性を設計する。整相処理では、音速ｃや位置ベクトルｒｍが時間変化する状況では、式（６）により、ｍ番目センサ出力に与える遅延時間τｍが変化する。この場合でも、振幅特性は時不変であり、センサ出力に与えるシェーディング係数に基づいて設計すればよい。一方、位相特性は時間変化するので、この時間変化部分に対応する必要がある。

本実施の形態１では、この遅延時間τｍのうちサンプリング周期以下の短い遅延時間、つまり詳細遅延δｔｍに対応するように位相特性を設計する。ここで、詳細遅延δｔｍは上述のように−１／２ｆｓ〜１／２ｆｓの範囲の任意の値を取り得るが、ここでは−１／２ｆｓ〜１／２ｆｓの範囲をＴ個に分割し（ｄＴ（ｔ），ｔ＝１，２，・・・Ｔ）、この遅延時間ｄＴ（ｔ）に対応する位相特性が設計される。このとき、分割数Ｔは、整相処理の方位精度に影響するので、アプリケーションで許容できる方位誤差以下となるように決定される。このようにして設計した振幅特性と位相特性との積を取ったものを、フィルタ周波数特性の設定値ＦｉｌｔｅｒＤ（ｄＴ（ｔ），ｋ）とする。

次に、変換手段３２は、作成したフィルタ周波数特性ＦｉｌｔｅｒＤ（ｄＴ（ｔ），ｋ）をＩＦＦＴすることで、作成したフィルタ周波数特性の時間領域表現であるフィルタ係数を作成する。ここで、フィルタ係数の適当な並び替えによって、このフィルタ係数の時間的な中心位置がＮ／２＋１となるようにする。

次に、抽出手段３３は、変換手段３２で作成したフィルタ係数から、Ｎ／２＋１の位置を中心として、（Ｎ‘−１）／２サンプル前後のフィルタ係数（合計Ｎ’サンプル）を抽出する。

そして、挿入手段３４は、抽出したフィルタ係数に対して、前後にゼロ詰めを行う。挿入手段３４は、インパルス応答の主要部の両端に、時間変化する位相応答の影響を抑える分のゼロデータを挿入する。フィルタ係数へのゼロ詰めの種類について、３つ例示する。

１つ目のゼロ詰めでは、抽出手段３３において抽出したフィルタ係数の前側にＮｆサンプルのゼロ詰めを行う。このＮｆサンプルは、位相回転手段４１において負のサンプル遅延を与えても、フィルタ係数の非ゼロ部分が循環しないようにするためのマージンである。Ｎｆの値は、整相処理で各センサ出力に与える負の遅延の最大値に基づいて決定される。

２つ目のゼロ詰めでは、抽出手段３３において抽出したフィルタ係数の後ろ側にＮｂサンプルのゼロ詰めを行う。このＮｂサンプルは、位相回転手段４１において正のサンプル遅延を与えても、フィルタ係数の非ゼロ部分が循環しないようにするためのマージンである。Ｎｂの値は、整相処理で各センサ出力に与える正の遅延の最大値に基づいて決定される。

３つ目のゼロ詰めは、Ｏｖｅｒｌａｐ−ｓａｖｅ法のためのゼロ詰めであり、１つ目のゼロ詰め及び２つ目のゼロ詰めの後に、Ｎ０サンプルをフィルタ係数にゼロ詰めする。このとき、Ｎ０サンプルをゼロ詰めする位置は、フィルタ係数の前でも後ろでもよいが、本提案書では、後ろにゼロ詰めしたとして説明する。なお、このＮ０サンプルのゼロ詰めの位置によって、図１の切り出し１３におけるデータの切り出し位置が変わる。

次に、ＦＦＴ３５において、このゼロ詰めしたフィルタ係数をＦＦＴすることで、Ｏｖｅｒｌａｐ−ｓａｖｅ法で用いるフィルタ周波数特性Ｆｉｌｔｅｒ（ｄＴ（ｔ），ｋ）を得る。このとき、得られたフィルタ周波数特性は、設定値ＦｉｌｔｅｒＤ（ｄＴ（ｔ），ｋ）と比べて誤差があるので、誤差が許容値以上である場合は構築手段３１に戻り、フィルタ周波数特性の振幅特性及び位相特性を設定し直す。誤差が許容値以内であれば、テーブル化手段３６は、フィルタ周波数特性Ｆｉｌｔｅｒ（ｄＴ（ｔ），ｋ）をテーブル化する。

図４に、フィルタ周波数特性テーブルに、フィルタ周波数特性Ｆｉｌｔｅｒ（ｄＴ（ｔ），ｋ）が格納されているイメージを示す。フィルタ周波数特性テーブルには、詳細遅延の範囲−１／２ｆｓ〜１／２ｆｓをＴ分割したときの、各詳細遅延に対応するフィルタ周波数特性Ｆｉｌｔｅｒ（ｄＴ（ｔ），ｋ）が並んでいる。

このようにして生成されたフィルタ周波数特性テーブルＦｉｌｔｅｒ（ｄＴ（ｔ），ｋ）から、（１０）式で得られた詳細遅延δｔｍに最も近い遅延のテーブルを選択し、ｔｂｌｍ（ｋ）を得る。

・・・（１１）

ただし、ｋ＝０，１，・・・，Ｎ−１である。ここでＴａｂｌｅ（Ａ，ｋ）は、Ａに最も近い値に対応する、Ｆｉｌｔｅｒ（ｄＴ（ｔ），ｋ）を選択する関数である。

（特性取得部２４，位相回転手段４１）
特性取得部２４は、位相回転手段４１を有する。位相回転手段４１は、フィルタ周波数特性ｔｂｌｍ（ｋ）とサンプル遅延数ｑｍを入力として、フィルタ周波数特性ｔｂｌｍ（ｋ）に対してサンプル遅延数ｑｍで決まる位相回転を実施し、最終的なフィルタ周波数特性（最終フィルタ周波数特性）Ｗｍ（ｋ）を求めて、フィルタ処理部２に出力する。位相回転手段４１は、フィルタ周波数特性テーブルから取得したフィルタ周波数特性ｔｂｌｍ（ｋ）に対し、サンプル遅延数ｑｍで決まる位相回転を実施し、ｍ番目のセンサ出力に与えるフィルタ周波数特性Ｗｍ（ｋ）を各周波数ビンｋ（ｋ＝０，１，・・・Ｎ−１）で計算する。

・・・（１２）

この一連の計算が全センサの出力に対して実施され、Ｍ個のセンサ出力に対して与えるフィルタ周波数特性が得られる。

実施の形態１では、遅延計算部２１及び分解部２２において、各センサ出力に与える遅延時間をサンプル遅延数と詳細遅延（詳細遅延時間）に分離される。フィルタ周波数特性テーブルでは、この詳細遅延に基づいて、あらかじめ計算されているフィルタ周波数特性をテーブルから取り出すことによって、フィルタ周波数特性を毎周期計算する必要がなくなる。

加えて、サンプル遅延に応じた位相回転において、フィルタ周波数特性に対しサンプル遅延数に基づいた位相回転を与えることで、このフィルタ周波数特性にＩＦＦＴ１２を適用することなく、このフィルタ周波数特性の時間領域表現であるインパルス応答に対して時間軸方向のサンプル遅延を与えることができる。これは、挿入手段３４においてフィルタ係数に与えた、ＮｆサンプルとＮｂサンプルによる効果である。

図５に、実施の形態１による計算量の削減について、例を示す。この例では、アレイを構成するセンサ数及びビーム数を６４とし、ＦＦＴ１１のブロック長Ｎと整相処理に要する複素乗算の回数の関係を示している。なお、ＦＦＴ１回あたりの複素乗算の回数を２Ｎｌｏｇ_２Ｎとして試算している。

実施の形態１の乗算回数は、図２のフィルタ係数計算部のＩＦＦＴ１２と、ゼロ詰め後のＦＦＴ１１とが発生しないために処理量が少ない。ＦＦＴ次数が１０２４の場合、従来技術に対して、実施の形態１における乗算回数は約１／１５に削減できている。ＦＦＴ次数が大きくなるほど処理量削減の効果は大きい、ＦＦＴ次数が６４〜４０９６の範囲では、乗算回数は従来技術に対して、１／１０〜１／１８の乗算回数となる。

したがって、実施の形態１を導入することにより、大幅な高速化を図ることができる。

実施の形態２．
実施の形態１においては、遅延時間が時々刻々変化するが、フィルタ周波数特性の振幅特性は時不変であり、遅延時間の変化が処理する帯域に渡って同一な（直線位相特性を満たす）整相処理を対象とした。このような通常の整相処理に対して、周囲環境に応じてフィルタ周波数特性を動的に変化させることで、よりＳ／Ｎを改善することができる整相処理も広く使われている。この整相処理を適応整相処理と呼ぶ。

適応整相処理では入力データに対して最適なフィルタ周波数特性を求めるため、振幅特性・位相特性共に時間変化である。実施の形態１は、フィルタ周波数特性の振幅特性が時不変で、センサに与える遅延時間の変化が全帯域で同じ場合に対応する構成であるため、フィルタ周波数特性の振幅特性及び位相特性が周波数ごとに時間変化する適応整相処理には対応できない。実施の形態２は、適応整相処理に対して、波形の連続性を維持しつつ高速化を実現する手法である。

実施の形態２の構成を図６に示す。図６において太い枠で囲んだ部分が発明部分である。適応重みの計算は、ＭＶＤＲ（Ｍｉｎｉｍｕｍ Ｖａｒｉａｎｃｅ Ｄｉｓｔｏｒｔｉｏｎｌｅｓｓ Ｒｅｓｐｏｎｓｅ）法、ＬＭＳ（Ｌｅａｓｔ Ｍｅａｎ Ｓｑｕａｒｅ）法又はＤＭＩ（Ｄｉｒｅｃｔ Ｍａｔｒｉｘ Ｉｎｖｅｒｓｉｏｎ）法等といった任意のアルゴリズムで実現することができる。本実施の形態２では、適応重みと発明部分以外の部分とは実施の形態１と同じのため、異なる部分についてのみ説明する。

フィルタ特性出力部１０３は、重み計算部１２１と、平均部１２２と、補正量計算部１２３と、テーブル作成部１２４と、特性取得部１２５とを有している。

（重み計算部１２１）
重み計算部１２１は、フィルタ処理部２のＦＦＴ３５からの出力に基づいて適応重みＷｏｐｔ（ｋ）を計算する。

（平均部１２２）
平均部１２２は、ｍ番目センサに与える適応重みＷｏｐｔ（ｋ）に対して、周波数方向に平均化した適応重みＷｏｐｔｍｅａｎｍ（ｌ）を後段に出力する。ここで、ｌはフィルタバンク番号である。各フィルタバンクの振幅特性と位相特性とを補正するために、適応重みＷｏｐｔ（ｋ）の周波数方向での平均処理を行い、各フィルタバンクの中心周波数の振幅と位相とを得る。フィルタバンク数はＬである。ここで、平均化した適応重みＷｏｐｔｍｅａｎｍ（ｌ）は以下の式で計算する（ただしｌ＝１，２，・・・Ｌ）。

・・・(１３)

Ｋ‘は平均処理に使用するビン数であり、ｋ＝０，２，・・・Ｎ−１である。（１３）式においてｋｌは１番目のフィルタバンクの中心周波数のビン番号である。ここで、フィルタバンク数Ｌ、平均処理に使用するビン数Ｋ’は、振幅特性・位相特性のフィルタバンク内の変化量が、アプリケーションで求められる許容誤差になるように選ぶ。

（補正量計算部１２３）
補正量計算部１２３では、平均化した適応重みＷｏｐｔｍｅａｎｍ（ｌ）を入力として、振幅補正量Ａｍ（ｌ）、サンプル遅延数ｑｍ（ｌ）、詳細遅延δｔｍ（ｌ）が出力される。この処理では、平均部１２２で計算したＷｏｐｔｍｅａｎｍ（ｌ）を用いて、各フィルタバンクの振幅及び位相補正量を計算する。各フィルタバンクに対する位相補償量に基づいて詳細遅延δｔｍ（ｌ）とサンプル遅延数ｑｍ（ｌ）とを計算するために、ｌ番目のフィルタバンクに対応した平均化適応重みＷｏｐｔｍｅａｎｍ（ｌ）の位相ｐｈａｓｅｏｐｔｍ（ｌ）を計算する。

・・・（１４）

（１４）式におけるａｎｇｌｅ（Ａ）はＡの位相角を計算する処理であり、ｕｎｗｒａｐ（Ａ）はＡの位相アンラップ処理である。位相アンラップとは、本来は−π〜πの範囲で定義される位相の値を、位相の連続性を考慮して、−π〜πより外側の領域へ拡張する処理である。ここで位相ｐｈａｓｅｏｐｔｍ（ｌ）をｌ番目のフィルタバンクに対する遅延時間τｍ（ｌ）に変換する。

・・・（１５）

ただし、ｆ（ｋｌ）はｋｌ番目の周波数ビンに相当する周波数である。このτｍ（ｌ）を遅延分解することでサンプル遅延数ｑｍ（ｌ）、詳細遅延δｔｍ（ｌ）に分解する。

・・・（１６）

・・・（１７）

一方、振幅補正量はＷｏｐｔｍｅａｎｍ（ｌ）の絶対値を用いる。

・・・（１８）

ここで、（１８）のａｂｓ（Ａ）はＡの絶対値を計算する処理である。

（テーブル作成部１２４）
テーブル作成部１２４は、詳細遅延δｔｍ（ｌ）に応じたフィルタ周波数特性ｔｂｌｍ（ｋ，ｌ）を位相回転手段４１に送る。詳細遅延δｔｍ（ｌ）に相当するフィルタ周波数特性をテーブルから読み出す。実施の形態１と比較すると、テーブルをフィルタバンクごとに読み出す点が実施の形態１と異なる。ここで、詳細遅延δｔｍ（ｌ）のフィルタ周波数特性は１番目の周波数を中心としたバンドパスフィルタであり、Ｏｖｅｒｌａｐ−ｓａｖｅ法のために時間領域でゼロ詰めされたフィルタ係数をＦＦＴしたものである。フィルタテーブル作成方法を図７に示す。変換手段３２と、抽出手段３３と、挿入手段３４と、ＦＦＴ３５と、テーブル化手段３６とは、それぞれ、実施の形態１の変換手段３２と、抽出手段３３と、挿入手段３４と、ＦＦＴ３５と、テーブル化手段３６と同様の処理であるので説明を省略する。

構築手段１３１は、フィルタ周波数特性の振幅特性と位相特性を設計する。位相特性の設定方法は、実施の形態１の図３と同様に、−１／２ｆｓ〜１／２ｆｓの範囲をＴ個に分割し（ｄＴ（ｔ），ｔ＝１，２，・・・Ｔ）、この遅延時間ｄＴ（ｔ）に対応する位相特性を設計する。一方、振幅特性は、ｌ番目の中心周波数ｋｌを中心とした振幅１のバンドパスフィルタを設定する。図８に、フィルタテーブルに、フィルタ周波数特性Ｆｉｌｔｅｒ（ｄＴ（ｔ），ｋ，ｌ）が格納されているイメージを示す。フィルタテーブルには、詳細遅延の範囲−１／２ｆｓ〜１／２ｆｓをＴ分割し、分割した中でさらに、フィルタバンク数Ｌ個だけ各詳細遅延に対応するフィルタ周波数特性Ｆｉｌｔｅｒ（ｄＴ（ｔ），ｋ，ｌ）が並んでいる。

このようにして生成されたフィルタテーブルから、（１７）式で得られた詳細遅延δｔｍ（ｌ）に最も近い遅延のテーブルを選択し、ｔｂｌｍ（ｋ，ｌ）を得る。

・・・（１９）

ここでＴａｂｌｅ２（Ａ，ｋ，ｌ）は、Ａに最も近い値に対応する、Ｆｉｌｔｅｒ（ｄＴ（ｔ），ｋ，ｌ）を選択する関数である。

（特性取得部１２５，位相回転手段１４１）
特性取得部１２５は、位相回転手段１４１と、振幅補正手段１４２と、加算手段１４３とを有する。位相回転手段４１は、フィルタ周波数特性ｔｂｌｍ（ｋ，ｌ）に対して、サンプル遅延数ｑｍ（ｌ）で決まる、位相回転を行ったフィルタ周波数特性Ｐｍ（ｋ，ｌ）を振幅補正手段１４２に送る。詳細遅延δｔｍ（ｌ）に対応するフィルタ周波数特性ｔｂｌｍ（ｋ，ｌ）に対してサンプル遅延数ｑｍ（ｌ）に相当する位相回転を行い、位相回転したフィルタ周波数特性Ｐｍ（ｋ，ｌ）を得る。

・・・（２０）

（振幅補正手段１４２）
振幅補正手段１４２は、位相回転を行ったフィルタ周波数特性Ｐｍ（ｋ，ｌ）に対して、振幅補正量Ａｍ（ｌ）を掛け合わせた、ｌ番目フィルタバンクで計算したフィルタ周波数特性Ｗｏｐｔｂａｎｄｍ（ｋ，ｌ）を加算手段１４３に送る。位相回転したフィルタ周波数特性Ｐｍ（ｋ，ｌ）に対して、（１８）式で求めた振幅補正量Ａｍ（ｌ）を掛け合わせて、ｌ番目フィルタバンクのフィルタ周波数特性Ｗｏｐｔｂａｎｄｍ（ｋ，ｌ）を得る。

・・・（２１）

（加算手段１４３）
加算手段１４３は、各フィルタバンクで計算したフィルタ周波数特性Ｗｏｐｔｂａｎｄｍ（ｋ，ｌ）を加算することで、最終的にｍ番目センサ出力に与えるフィルタ周波数特性Ｗｏｐｔｏｌｓｍ（ｋ）をフィルタ処理部２に送る。フィルタバンクのフィルタ周波数特性は、各センサ出力に対してＬ個計算される。このＬ個のフィルタバンクのフィルタ周波数特性を加算することで、最終的にｍ番目センサ出力に与えるフィルタ周波数特性を得る。

・・・（２２）

この一連の計算を全センサに対して実施し、Ｍ個のセンサ出力に対して与えるフィルタ周波数特性を得る。

実施の形態２．
次に、実施の形態２の動作について説明する。ここでは、例として、一つの整相方位を持った整相器における場合を述べる。

実施の形態２は、フィルタ周波数特性の振幅特性及び位相特性が周波数ごとに時間変化するような適応整相処理に対応するものである。そこで、処理帯域内を複数のフィルタバンクに分割し、各フィルタバンクでＯｖｅｒｌａｐ−ｓａｖｅ法を適用する。

実施の形態２における、フィルタバンクによる考え方を、図９、図１０に示す。詳細については、各動作の説明で述べることとし、ここでは実施の形態２の概要を示す。図９において示している、９−１で示す適応重みの振幅特性は時間変化であり、位相特性についても時間変化であり、さらに周波数ごとに遅延時間が異なる。まず、狭い帯域幅の区間では振幅特性・位相特性ともに大きく変化しないという前提の下に、９−２にあるように振幅特性と周波数特性とをＬ個の帯域に分割する。

次に、９−３にあるように各帯域で平均処理を実施するが、これは平均部１２２で実施される処理である。平均化された適応重みの振幅と位相について、各フィルタバンクでＯｖｅｒｌａｐ−ｓａｖｅ法を実現することを考える。９−４において、ｌ番目フィルタバンクの平均化された適応重みの位相の値を、遅延時間に変換する。遅延時間を、サンプリング周期以下の短い遅延時間である詳細遅延と、サンプリング周波数の逆数の整数倍で表現できるサンプル遅延に分割する。詳細遅延に対しては、その値に応じて、９−５にあるように事前にテーブル化した、ｌ番目フィルタバンクのフィルタ周波数特性を読み出す。テーブルから参照したフィルタ周波数特性に対して、サンプル遅延を基にした、９−６の位相回転を実施したあと、９−７にあるように、ｌ番目フィルタバンクの平均化された適応重みの振幅の値を掛け合わせる。以上の動作により、適応重みの振幅・位相を再現した、ｌ番目フィルタバンクのフィルタ周波数特性を得る。

図１０は、図９に示した手法で作成した、各フィルタバンクのフィルタ周波数特性を加算するイメージを記す。１０−１は各フィルタバンクで作成したフィルタ周波数特性であり、これらを加算することにより、１０−２に示すように、最終的にセンサ出力に与えるフィルタ周波数特性を得る。なお、図１０は説明のため振幅特性のみ図示しているが、フィルタバンクのフィルタ周波数特性の加算は複素加算で実施する。これは加算手段１４３に対応する。

実施の形態１では時間変化し、かつ、周波数ごとに振幅特性／位相特性が変わる適応整相のフィルタ係数に対応できなかったが、実施の形態２でフィルタバンクごとに補正を行うように修正したことで、適応整相の高速化が可能になる。

補正量計算部１２３において、位相回転したフィルタ周波数特性Ｐｍ（ｋ，ｌ）に対して、定数（Ａｍ（ｌ））を掛け合わせているだけなので、Ｗｏｐｔｂａｎｄｍ（ｋ，ｌ）の逆ＦＦＴ結果であるフィルタ係数も定数倍されるだけである。そのため、Ｐｍ（ｋ，ｌ）同様、波形の連続を維持し、かつ適応重みの特性が付与されたフィルタ周波数特性になっている。また、加算手段１４３において、最終的にｍ番目センサ出力に与えるフィルタ周波数特性をＷｏｐｔｏｌｓｍ（ｋ）も、波形の連続性を維持したフィルタバンクの加算結果であるため、波形の連続性を維持している。

図１１に、実施の形態２による計算量の削減効果について、例を示す。実施の形態１と同じ評価を実施する。適応整相は、ビーム数をセンサ数よりも多くする必要があるので、アレイを構成するセンサ数を６４、ビーム数を３２０(センサ数の５倍)とした。

実施の形態２の乗算回数は、従来技術の乗算回数に対し、１／１１〜１／２０の乗算回数で処理を実現できる。また、従来技術の通常の整相処理と比べると、乗算回数が１／２〜１／４程度になり、従来技術の通常の整相処理よりも乗算回数を減らすことができる。

実施の形態１において、フィルタ周波数特性テーブルは振幅特性と位相特性を分割してテーブル化することができる。実施の形態１は、振幅特性(パターン数)×位相特性（遅延分割数）のテーブルを持つ必要があるが、分割すると振幅特性(パターン数)＋位相特性（遅延分割数）だけのテーブルを持てばよい。ただし、この場合振幅特性と位相特性の周波数領域上での掛け算になり、時間領域では振幅特性と位相特性の逆ＦＦＴ結果であるフィルタ係数同士の畳み込み演算となる。そのため、畳み込み演算をした結果が図５のＮ０に漏れこまないように、振幅特性を計算するときのＮｆ、Ｎｂ、そして位相特性を計算するときのＮｆ、Ｎｂを設計する必要がある。

また、実施の形態１・実施の形態２ともにサンプル遅延量による位相回転量を毎回計算する構成になっている。この位相回転量は、取り得る値を前もって見積もることができるため、事前に計算した位相回転量をテーブル化し、位相回転の際にはテーブルから読み出すことも可能である。

実施の形態１・実施の形態２ともに、Ｏｖｅｒｌａｐ−ｓａｖｅ法をベースにした高速化手法である。ＦＦＴの高速性を利用した、波形の連続性を保つ別手法として、Ｏｖｅｒｌａｐ−ａｄｄ法がある。本手法はＯｖｅｒｌａｐ−ａｄｄ法でも、提案部分であるフィルタの作成方法を変更することなく実現可能である。Ｏｖｅｒｌａｐ−ａｄｄ法を用いる場合、入力データの波形にゼロ詰めを行ったのちＦＦＴを実施し、ＩＦＦＴ後に切り出しではなく、出力データを前回出力データと所定の部分をオーバーラップして加算することで実現可能である。

１ フィルタリング装置、２ フィルタ処理部、３ フィルタ特性出力部、１１ ＦＦＴ、１２ ＩＦＦＴ、１３ 切り出し、２１ 遅延計算部、２２ 分解部、２３ テーブル作成部、２４ 特性取得部、３１ 構築手段、３２ 変換手段、３３ 抽出手段、３４ 挿入手段、３５ ＦＦＴ、３６ テーブル化手段、４１ 位相回転手段、１００ フィルタリング装置、１０３ フィルタ特性出力部、１２１ 重み計算部、１２２ 平均部、１２３ 補正量計算部、１２４ テーブル作成部、１２５ 特性取得部、１３１ 構築手段、１４１ 位相回転手段、１４２ 振幅補正手段、１４３ 加算手段。

Claims (7)

1. 複数のセンサからの出力に対するフィルタリング装置であって、
   フィルタの位相特性に与える遅延時間として、遅延時間の範囲を複数の区間に分割し、分割された各区間において遅延時間に対応する位相特性を求め、所定の振幅を有するバンドパスフィルタを周波数領域上で構築する構築手段と、
   前記構築手段によって構築されたバンドパスフィルタの周波数特性に基づいて時間領域上のインパルス応答を求める変換手段と、
   前記変換手段によって求められた時間領域上のインパルス応答の主要部を抽出する抽出手段と、
   前記抽出手段によって抽出されたインパルス応答の主要部の両端にゼロデータを挿入する挿入手段と、
   前記挿入手段によってゼロデータが挿入されたインパルス応答に基づいて得られる周波数特性をテーブルとするテーブル化手段と、を有するテーブル作成部と、
   前記周波数特性のテーブルに基づいて、フィルタの位相特性に与える遅延時間のうちサンプリング周波数以下の詳細遅延時間に対応する周波数特性を取得し、前記サンプリング周波数の整数倍のサンプル遅延時間に位相回転を実施して最終フィルタ周波数特性を求める特性取得部と、を備える
   ことを特徴とするフィルタリング装置。
2. 前記フィルタの位相特性に与える遅延時間を計算する遅延計算部と、
   前記遅延計算部によって計算された遅延時間を、前記サンプル遅延時間と、前記詳細遅延時間とに分解する分解部と、を更に備え、
   前記特性取得部は、
   前記分解部によって分解された詳細遅延時間に対応する周波数特性を取得し、前記分解部によって分解されたサンプル遅延時間に位相回転を実施して最終フィルタ周波数特性を求める位相回転手段を有する
   ことを特徴とする請求項１記載のフィルタリング装置。
3. 複数のセンサからの出力に対するフィルタリング装置であって、
   前記構築手段は、
   周波数帯域内で複数に分割されたフィルタバンク毎に、フィルタの位相特性に与える遅延時間として、遅延時間の範囲を複数の区間に分割するものである
   ことを特徴とする請求項１又は２記載のフィルタリング装置。
4. 周波数領域で求められた重み係数を前記フィルタバンクに分割して、前記フィルタバンク内で重み係数の平均値を算出する平均部と、
   前記平均部によって算出された重み係数の平均値に基づいて、振幅補正量と、前記サンプル遅延時間と、前記詳細遅延時間とを求める補正量計算部と、を更に備え、
   前記特性取得部は、
   前記補正量計算部によって求められた詳細遅延時間に対応する周波数特性を取得し、前記補正量計算部によって求められたサンプル遅延時間に位相回転を実施する位相回転手段と、
   前記位相回転手段によって位相回転されたサンプル遅延時間と、前記補正量計算部によって求められた振幅補正量とに基づいて、前記フィルタバンクの重み係数を求める振幅補正手段と、
   前記振幅補正手段によって求められた重み係数により重み付けされたフィルタバンクに基づいて最終フィルタ周波数特性を求める加算手段と、を有する
   ことを特徴とする請求項３記載のフィルタリング装置。
5. 前記挿入手段は、
   前記インパルス応答の主要部の両端に、時間変化する位相応答の影響を抑える分のゼロデータを挿入し、
   前記ゼロデータが挿入されたインパルス応答に対し、Ｏｖｅｒｌａｐ−ｓａｖｅ法又はＯｖｅｌａｐ−ａｄｄ法に必要なゼロデータを挿入する
   ことを特徴とする請求項１〜４のいずれか１項に記載のフィルタリング装置。
6. 複数のセンサからの出力に対するフィルタリング装置のテーブル作成方法であって、
   フィルタの位相特性に与える遅延時間として、遅延時間の範囲を複数の区間に分割し、分割された各区間において遅延時間に対応する位相特性を求め、所定の振幅を有するバンドパスフィルタを周波数領域上で構築するステップと、
   構築されたバンドパスフィルタの周波数特性に基づいて時間領域上のインパルス応答を求めるステップと、
   求められた時間領域上のインパルス応答の主要部を抽出するステップと、
   抽出されたインパルス応答の主要部の両端にゼロデータを挿入するステップと、
   ゼロデータが挿入されたインパルス応答に基づいて得られる周波数特性をテーブルとするステップと、を有する
   ことを特徴とするフィルタリング装置のテーブル作成方法。
7. 複数のセンサからの出力に対するフィルタリング装置のテーブル作成方法であって、
   周波数帯域内で複数に分割されたフィルタバンク毎に、フィルタの位相特性に与える遅延時間として、遅延時間の範囲を複数の区間に分割するステップを更に有する
   ことを特徴とする請求項６記載のフィルタリング装置のテーブル作成方法。

Images