音響信号フィルタとそのフィルタリング方法と、そのプログラムと記録媒体
【課題】分析窓長よりも長い有限長インパルス応答フィルタを時間領域において音響信号に畳み込む演算と等価な演算を、周波数領域において実現する。
【解決手段】この発明の音響信号フィルタの巡回時間シフト演算部は、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数フィルタのフレーム列に、巡回時間シフト演算子を乗算して周波数フィルタを生成する。畳み込み演算部は、離散値化された時間領域音響信号のフレーム長Ｍ＋Ｎ−１の各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する。窓関数適用部は、全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する。
【解決手段】この発明の音響信号フィルタの巡回時間シフト演算部は、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数フィルタのフレーム列に、巡回時間シフト演算子を乗算して周波数フィルタを生成する。畳み込み演算部は、離散値化された時間領域音響信号のフレーム長Ｍ＋Ｎ−１の各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する。窓関数適用部は、全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、例えば収音信号に有限長インパルス応答フィルタを畳み込んで雑音を抑圧するのに用いられる音響信号フィルタとそのフィルタリング方法と、そのプログラムと記録媒体に関する。
【背景技術】
【０００２】
ディジタル信号処理における音響信号のフィルタは、例えば短時間フーリエ変換領域で音響信号に有限長インパルス応答フィルタを畳み込んで実現される。一般に有限長インパルス応答フィルタのフィルタ長が短時間フーリエ変換の分析窓長より十分に短い場合、周波数領域におけるフィルタの畳み込み計算は、周波数信号と周波数フィルタの積として計算することができる。
離散時間のインデックスをｔ、時間領域の音響信号をｘｔ、時間領域の有限長インパルス応答フィルタをｆｔと表わすものとする。時間シフトしながら短時間分析窓を適用することで切り出される時間フレームのインデックスをｎ、周波数信号の周波数ビンに関するインデックスをｍ、Ｘｎ，ｍを音響信号ｘｔの短時間フーリエ変換のフレームｎにおける周波数ビンｍの値、Ｆｍを周波数ビンｍにおける有限長インパルス応答フィルタの短時間フーリエ変換の値とする。音響信号ｘｔに有限長インパルス応答フィルタｆｔを畳み込んで得られるフィルタリングされた信号ｙｔの短時間フーリエ変換値Ｙｎ，ｍは式（１）で近似的に計算することができる。
【０００３】
【数１】
この短時間フーリエ変換と同様の手法である高速フーリエ変換を用いたフィルタ１０００の機能構成を、図１０に示して簡単に説明する（非特許文献１参照）。従来の高速フーリエ変換を用いたフィルタ１０００は、高速フーリエ変換部１０１と、高速フーリエ変換部１０２と、畳み込み演算部１０３と、高速フーリエ逆変換部１０４を備える。高速フーリエ変換部１０１は、音響信号ｘｔを高速フーリエ変換して上記した周波数ビンＸｎ，ｍを生成する。高速フーリエ変換部１０２は、有限長インパルス応答フィルタｆｔを高速フーリエ変換して上記した有限長インパルス応答フィルタの短時間フーリエ変換の値Ｆｍを生成する。畳み込み演算部１０３は、上記した式（１）を演算してフィルタリングされた信号ｙｔの短時間フーリエ変換値Ｙｎ，ｍを生成する。高速フーリエ逆変換部１０４は、短時間フーリエ変換値Ｙｎ，ｍを上記した時間領域の信号ｙｔに変換する。
【非特許文献１】三谷政昭著「ディジタルフィルタデザイン」株式会社昭晃堂、平成６年４月２０日初版７刷発行、ｐ97
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかし、従来の方法では、短時間フーリエ変換領域でフィルタの畳み込みを計算するためには、フィルタ長よりも十分に長い分析窓を使う必要がある。フィルタ長より分析窓長が短い場合は、長い有限長インパルス応答フィルタの畳み込みを周波数領域に閉じて精度良く計算することができなかった。単純に分析窓長を大きくしても、分析の時間分解能が低くなるので音声信号や音楽信号のように短時間で変化する信号を精度良くフィルタリングすることができない。また、いったん時間領域の信号に変換し直して畳み込み計算してから再び周波数領域の信号に戻す方法も考えられるが、計算量が増えてしまう。また、二乗誤差最小基準のような何らかのコスト関数に基づく最適化処理を行なう場合は、異なる領域における演算が必要であると共にコスト関数の最適化の計算が複雑になる課題がある。
【０００５】
この発明は、このような点に鑑みてなされたものであり、分析窓長よりも長い有限長インパルス応答フィルタを時間領域において音響信号に畳み込む演算と等価な演算を、周波数領域において実現する音響信号フィルタとそのフィルタリング方法と、プログラムと記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００６】
この発明の音響信号フィルタは、巡回時間シフト演算部と、畳み込み演算部と、窓関数適用部を具備する。巡回時間シフト演算部は、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数フィルタのフレーム列に、巡回時間シフト演算子を乗算して周波数フィルタを生成する。畳み込み演算部は、離散値化された時間領域音響信号のフレーム長Ｍ＋Ｎ−１の各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する。窓関数適用部は、全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する。そして、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２である。
【発明の効果】
【０００７】
この発明の音響信号フィルタは、分析窓長よりも長い有限長インパルス応答フィルタを時間領域において音響信号に畳み込む演算と等価な演算を、周波数領域において可能にする。このフィルタリング方法は、周波数領域においても畳み込みとして計算できるので、効率的な計算アルゴリズムにすることができる。
【発明を実施するための最良の形態】
【０００８】
以下、この発明の実施の形態を図面を参照して説明する。複数の図面中同一のものには同じ参照符号を付し、説明は繰り返さない。
〔発明の基本的な考え〕
実施例の説明の前に、この発明の音響信号フィルタリング方法の基本的な考えを説明する。まず、時間領域信号に対する畳み込み演算を周波数領域で計算するための関係式を導出する。
ディジタル化した離散時間信号の標本インデックスをｔとし、時間領域の音響信号ｘｔにタップ長Ｋの有限長インパルス応答フィルタｃｔ（ｔ＝１〜Ｋ）を時間領域で畳み込んで得られる出力信号をｙｔとする。時刻ｔ０で始まる出力信号ｙｔの短時間セグメントを時間窓で切り出した信号は、ｚ変換領域において式（２）で表現できる。
【０００９】
【数２】
ここで、ｙ（ｚ）＝ｃ（ｚ）ｘ（ｚ）とし、ＷＮ（・）は時間領域における長さＮの窓関数に相当する関数とした。ＷＮ（ｃ（ｚ））は、ｃ（ｚ）中のｚの−Ｎ+１次から０次の項を取り出し、窓の形に比例して各係数を変更して窓の外の項を除去する。ｚｔ0は、時刻ｔ０で始まるフレームを窓関数の中に移動する時間シフト演算子である。さらに、フィルタｃｔのｔ０タップ目から長さＭのフレームを切り出すことを式（３）で表現する。
【００１０】
【数３】
ここでＷＭＲ（・）は、長さＭの方形窓を表わすものとする。するとｃ（Ｚ）は式（４）で表わせ、式（２）は以下のように書き換えることができる。
【００１１】
【数４】
ここで、Ｋｃはタップ長Ｋを方形窓長Ｍで除したＫｃ＝「Ｋ/Ｍ」とした。「Ｋ/Ｍ」は、
Ｋ/Ｍ以上の最小の整数を表わすものとする。また、Ｍはフレームシフトに対応する。式（５）中で窓関数の引数に含まれる項の内、窓の外の項を除去することで式（６）が導かれる。式（６）の窓関数中の引数の内、式（７）で示す引数は各フレーム同士の時間領域での畳み込みに相当する。
【００１２】
【数５】
つまり、フィルタｃｔのτＭタップ目から長さＭのフレームを切り出し、音響信号ｘｔの時刻ｔ０−Ｍ＋１−τＭから長さＭ＋Ｎ−１のフレームを切り出し、各フレームをｚ領域で乗算したものである。よく知られているように、ｚ領域における乗算は、時間領域における畳み込み演算に相当する。さらに、時間領域における音響信号のフレームとフィルタのフレームの畳み込み演算は、短時間フーリエ変換のポイントをＮｐ、フィルタｃｔのフレーム長をＭ、音響信号ｘｔのフレーム長をＭ＋Ｎ−１としたときに、Ｎｐ≧２Ｍ＋Ｎ−２であるとき、短時間フーリエ変換領域での積で正確に表現することができる。次に、式（６）の窓関数中の引数のうち、ｚＭ−１はｚ領域におけるＭ−１次の時間シフト演算子であり、時間領域においてサンプルの時間インデックスをＭ−１だけ小さくすることに相当する。これは短時間フーリエ変換領域では、Ｍ−１次の巡回時間シフト演算子との積として表わすことができる。本明細書では、これをＧ（Ｍ−１）と書くことにする。さらに、式（６）中のＷＮは、時間領域における窓関数に相当するので、短時間フーリエ変換領域においては、窓関数の短時間フーリエ変換との巡回畳み込みで表現できることが知られている。この巡回畳み込みは、後述するように巡回畳み込み行列との積として表現できる。最後に、式（６）中の総和演算Σは、短時間フーリエ変換領域においても総和演算として表現することができる。また、この総和演算において、τが１増加するに応じて、式（７）の二つのフレームはインデックスがＭずつ増加する。これは、Ｍが上記短時間フーリエ変換の計算におけるフレームシフトに対応していることを意味している。
【００１３】
この発明の音響信号フィルタリング方法は、上記した関係を用いることで時間領域での分析窓長よりも長い有限長インパルス応答フィルタの畳み込み計算を、周波数領域において可能にしたものである。
なお、この発明の説明における短時間フーリエ変換は式（８）と式（９）で定義されるものとする。
【００１４】
【数６】
ここで、Ｔは行列またはベクトルの転置を表わし、Ｘτ、ｍは、時刻ｔτで始まる時間フレームτに関する音響信号ｘｔの短時間フーリエ変換の周波数ビンｍにおける値を表わす。Ｘτは、フレームτにおける全周波数ビンの値を並べたベクトルとする。
【実施例１】
【００１５】
図１にこの発明の音響信号フィルタリング方法を用いた音響信号フィルタ１００の機能構成例を実施例１として示す。その動作フローを図２に示す。音響信号フィルタ１００は、巡回時間シフト演算部１０と、畳み込み演算部１２と、窓関数適用部１４を備える。この例の音響信号フィルタ１００は、例えばＲＯＭ、ＲＡＭ、ＣＰＵ等で構成されるコンピュータに所定のプログラムが読み込まれて、ＣＰＵがそのプログラムを実行することで実現されるものである。以降で説明する音響信号フィルタも同様である。
巡回時間シフト演算部１０は、時間領域の有限長インパルス応答フィルタのフレーム長Ｍ、フレームシフトＭで切り出される各フレームを、それぞれポイント数Ｎｐにより短時間フーリエ変換した周波数フィルタのフレーム列Ｃｎ（ｎ＝１〜Ｋｃ）を入力として、式（１０）に示す演算を行い周波数フィルタＣｎ〜（ｎ＝１〜Ｋｃ）を生成する（ステップＳ１０）。各ｎに対するＣｎは、式（９）と同様、各周波数ビンを要素にもつベクトルである。
【００１６】
【数７】
式（１０）は各周波数ビン毎の周波数フィルタのフレーム列Ｃｎに、ｋ＝Ｍ−１として巡回シフト演算子Ｇ（ｋ）を乗算するものである。ここでｄｉａｇ（Ｘ）は、ベクトルＸを対角要素に持つ対角行列を表わす。巡回時間シフト演算子Ｇ（ｋ）は、周波数領域の表現で式（１１）で表わされるものである。
【００１７】
【数８】
畳み込み演算部１２は、時間領域の音響信号ｘｔのフレーム長Ｍ＋Ｎ−１、フレームシフトＭで切り出される各フレームを、ポイント数Ｎｐにより短時間フーリエ変換した周波数信号のフレーム列を入力として、周波数ビン毎に周波数フィルタＣｎ〜を畳み込み演算して全周波数ビンをまとめたベクトルＹｎを式（１２）の演算により生成する（ステップＳ１２）。
【００１８】
【数９】
窓関数適用部１４は、式（１３）に示す様に全周波数ビンをまとめたベクトルＹｎに窓関数の短時間フーリエ変換値Ｈを一列目に持つエルミートテプリッツ行列ＷＮを乗じて周波数領域出力信号ＷＮＹｎを生成する（ステップＳ１４）。
いま、Ｈを窓関数の短時間フーリエ変換の各周波数ビンを並べたベクトルとし、式（１３）に示すように定義する。
【００１９】
【数１０】
また、エルミートテプリッツ行列ＷＮを式（１４）に示すように定義する。
【００２０】
【数１１】
ここで＊は複素共役を表わす。
【００２１】
【数１２】
周波数領域でエルミートテプリッツ行列ＷＮを乗ずることは、時間領域で対応する時間窓を適用することと等価な演算である。
なお、窓関数は一般に低域通過フィルタの特性を示すので、窓関数の短時間フーリエ変換Ｈは、低域周波数にのみ大きなゲインを持つ。このため、周波数領域における窓関数の演算においては、高域周波数に関する計算を省略しても大きな誤差を生じない。これは、ｍ＝０またはｍ＝Ｎｐ-1近傍の（式（８）を参照）の周波数ビン以外の周波数ビンに対応するＨの値は０と近似しても良いことを意味する。このように近似することで、エルミートテプリッツ行列ＷＮの積はＨの低域周波数近傍の周波数ビンのみに関する畳み込み演算とすることができ、効率的に計算することができる。
【００２２】
以上述べたように、この実施例によれば、時間領域の畳み込み演算を周波数領域における周波数ビン毎の畳み込み演算で表現することができ、時間領域での分析窓長よりも長い有限長インパルス応答フィルタの畳み込み計算を周波数領域において可能にする。この実施例で説明したフィルタリング方法は、周波数領域でも畳み込みとして計算できるので効率的な計算アルゴリズムにすることができる。
なお、周波数領域出力信号ＷＮＹｎを時間領域の信号に変換する場合は、図１に破線で示すように周波数領域出力信号ＷＮＹｎを時間領域の信号に変換する逆短時間フーリエ変換部１９と、時間領域の短時間フレームの時系列をオーバラップ加算するオーバラップ加算部２０とをさらに備え、それぞれで逆短時間フーリエ変換過程（ステップＳ１９）とオーバラップ加算過程（ステップＳ２０）の処理を行なえばよい。これらは、従来技術で容易に実現できるので説明は省略する。
〔変形例１〕
図３に実施例１を変形した音響信号フィルタ３００の機能構成例を示す。その動作フローを図４に示す。音響信号フィルタ３００は、周波数音響信号のフレーム列に巡回時間シフト演算子を乗算して周波数信号を生成するようにし、有限長インパルス応答フィルタ係数の周波数ビン毎に、巡回時間シフト演算子を乗算した周波数信号を畳み込み演算するようにしたものである。
【００２３】
音響信号フィルタ３００は、巡回時間シフト演算部３０と、畳み込み演算部３２と、窓関数適用部１４を備える。窓関数適用部１４は実施例１と同じものである。巡回時間シフト演算部３０は、離散値化された時間領域の音響信号のフレーム長Ｍ、フレームシフトＭで切り出した各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数音響信号のフレーム列Ｘｎに、ｋ＝Ｍ−１として巡回時間シフト演算子Ｇ（ｋ）を乗算して周波数音響信号Ｘｎ〜を生成する（式（１６））（ステップＳ３０）。
【００２４】
【数１３】
畳み込み演算部３２は、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍ＋Ｎ−１、フレームシフトＭで切り出した各フレームに対して、上記ポイント数Ｎｐで短時間フーリエ変換を適用して得られる周波数フィルタＣτ（τ＝１〜Ｋｃ）を、各周波数ビン毎に、上記周波数信号に畳み込み演算して全周波数ビンをまとめたベクトル生成する（式（１７））（ステップＳ３２）。
【００２５】
【数１４】
窓関数適用部１４は、実施例１と同じ動作を行い周波数領域出力信号ＷＮＹｎを生成する。以上述べたように、周波数音響信号のフレーム列に巡回時間シフト演算子を乗算して周波数信号を生成するようにしても、実施例１と同様に時間領域の畳み込み演算を周波数領域における周波数ビン毎の畳み込み演算で表現することができる。
【実施例２】
【００２６】
音響信号の時間窓長（フレーム長）Ｎを固定としたとき、有限長インパルス応答フィルタのフレーム長Ｍ（フレームシフトと同じ値）を大きくするほどフィルタ長Ｋｃを小さくできる。したがって上記した式（１２）の総和における加算の回数を小さくでき、効率的なフィルタリング計算が行える。ただし、ＭがＮよりも大きくなると、全周波数ビンをまとめたベクトルＹｎのフレームシフトが時間窓長Ｎよりも大きくなることに相当し、飛び飛びの時間フレームに対応するＹｎしか求められなくなる。そこで、効率的にフィルタリング計算を行うためには、Ｎ＝Ｍとすることが一つの解決策となる。Ｎ＝Ｍとすることで、実施例１で示した式（６）は式（１８）に書き換えることができる。
【００２７】
【数１５】
ｆｔは、同じ長さＮの時間フレーム同士の畳み込みの和である。短時間フーリエ変換のポ
イント数Ｎｐが、Ｎｐ≧２Ｎ−１のとき、短時間フーリエ変換の積の和として表現できる。
また、ｚＮは、Ｎ次のオーダの時間シフトオペレータであり、Ｎｐ≧２Ｎの短時間フーリエ変換を用いる場合、時間シフト計算が適切に行えることを意味する。
時間領域の窓関数が、周波数領域では巡回畳み込みで表現されることを考慮すると、式
（１８）は式（２０）に書き換えることができる。
【００２８】
【数１６】
音響信号フィルタ５００は式（２０）に示す動作を行なうようにしたものである。音響
信号フィルタ５００の機能構成例を図５にその動作フローを図６に示す。音響信号フィル
タ５００は、巡回シフト演算部５０と、畳み込み演算部５２と、窓関数適用部１４を備え
る。窓関数適用部１４は実施例１と同じものである。
巡回時間シフト演算部５０は、離散値化された時間領域の音響信号のフレーム長Ｎ、フレームシフトＭで切り出した各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数信号Ｘｎ〜を生成する（式（２１））（ステップＳ５０）。
【００２９】
畳み込み演算部５２は、時間領域の有限長インパルス応答フィルタ係数からフレーム長Ｍ、フレームシフトＭで切り出した各フレームに対して、ポイント数Ｎｐで短時間フーリエ変換を適用して得られる周波数フィルタを、周波数ビン毎に、連結周波数信号に畳み込み演算して全周波数ビンをまとめたベクトルＹｎを生成する（式（２０））（ステップＳ５２）。
このように構成しても実施例１と同様に時間領域の畳み込み演算を周波数領域における周波数ビン毎の畳み込み演算で表現することができる。Ｎ＝Ｍとすることで、音響信号フィルタ５００の入力信号のフレーム長と出力信号である周波数出力信号のフレーム長が同一になるので計算上の都合が良い。
【００３０】
なお、実施例２ではフレームシフト幅Ｍと窓長Ｎとが一致するため、隣接する時間窓がオーバラップしない。オーバラップしない信号を、そのまま時間領域の信号に変換すると雑音が増加する場合がある。それを防止するためには、例えば、時間領域信号から周波数信号の時系列を生成する際に適用する短時間フーリエ変換の時間シフトを時間窓長の半分にする。そして、奇数番目及び偶数番目のフレームインデックスを持つ二種類のフレーム系列を取り出し、各系列に別々に周波数フィルタを畳み込んだ後に一つの系列に戻すことで、隣接する時間窓が２分の１ずつオーバラップするフレーム系列を作ることができる。
〔変形例２〕
図７に実施例２を変形した音響信号フィルタ７００の機能構成例を示す。その動作フローを図８に示す。音響信号フィルタ７００は、実施例１に対する変形例１と同じように、巡回時間シフト演算子を乗算する信号系列を変更したものである。音響信号フィルタ７００は、巡回シフト演算部７０と、畳み込み演算部７２と、窓関数適用部１４を備える。窓関数適用部１４は実施例１と同じものである。
【００３１】
巡回時間シフト演算部７０は、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｎ、フレームシフトＭで切り出した各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタの各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数フィルタＣｎ〜を生成する（式（２２））（ステップＳ７０）。
畳み込み演算部７２は、時間領域の音響信号のフレーム長Ｎの時系列からフレーム長Ｍ、フレームシフトＭで切り出した各フレームに対して、ポイント数Ｎｐで短時間フーリエ変換を適用して得られる周波数音響信号を、周波数ビン毎に、連結周波数フィルタＣｎ〜に畳み込み演算して全周波数ビンをまとめたベクトルＹｎを生成する（式（２３））（ステップＳ７２）。
【００３２】
【数１７】
このように構成しても実施例２と同様に、周波数領域における周波数ビン毎の畳み込み演算を効率良くおこなうことができる。
【００３３】
〔シミュレーション結果〕
この発明の音響信号フィルタの動作を確認するシミュレーションを行った。シミュレー
ションは、長さ２０４８サンプルの時系列データを入力とし、長さ１０２４サンプルの有
限長インパルス応答フィルタを−１〜１の範囲乱数で生成する条件で行った。その結果を
図９に示す。図９（ａ）は時間領域の畳み込み演算でフィルタリングした場合を示す。図
９（ｂ）は有限長インパルス応答フィルタの時間窓長Ｍ＝６４、音響信号の時間窓長Ｎ＝
２５６として実施例１の方法でフィルタリングした場合を示す。図９（ｃ）はＭ＝２５６、Ｎ＝２５６として実施例２の方法でフィルタリングした場合を示す。それぞれの横軸は、時系列データのサンプル番号、縦軸は振幅である。
【００３４】
図９（ｂ）と（ｃ）の振幅０に引かれた破線は、時間領域畳み込み演算のフィルタリン
グ結果との差分を示す。このように実施例１と２のどちらのフィルタリング方法でも時間
領域畳み込み演算でフィルタリングしたのと同じ結果が得られた。
このようにこの発明のフィルタリング方法によれば、分析窓長よりも長い有限長インパ
ルス応答フィルタを時間領域において音響信号に畳み込む演算と等価な演算を、周波数領
域において可能にする。
なお、この発明の音響信号フィルタリング方法及び音響信号フィルタは上述の実施形態に限定されるものではなく、この発明の趣旨を逸脱しない範囲で適宜変更が可能であり、この発明の技術思想を、音響信号以外の２系統のディジタル信号系列に対して適用することも可能である。また、上記方法及び装置において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【００３５】
また、上記装置における処理手段をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、各装置における処理手段がコンピュータ上で実現される。
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）/ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【００３６】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記録装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
また、各手段は、コンピュータ上で所定のプログラムを実行させることにより構成することにしてもよいし、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図面の簡単な説明】
【００３７】
【図１】この発明の音響信号フィルタ１００の機能構成例を示す図。
【図２】音響信号フィルタ１００の動作フローを示す図。
【図３】この発明の音響信号フィルタ３００の機能構成例を示す図。
【図４】音響信号フィルタ３００の動作フローを示す図。
【図５】この発明の音響信号フィルタ５００の機能構成例を示す図。
【図６】音響信号フィルタ５００の動作フローを示す図。
【図７】この発明の音響信号フィルタ７００の機能構成例を示す図。
【図８】音響信号フィルタ７００の動作フローを示す図。
【図９】シミュレーション結果を示す図であり、（ａ）は時間領域の畳み込み演算でフィルタリングした場合、（ｂ）は有限長インパルス応答フィルタの時間窓長Ｍ＝６４、音響信号の時間窓長Ｎ＝２５６として実施例１の方法でフィルタリングした場合、（ｃ）はＭ＝２５６、Ｎ＝２５６として実施例２の方法でフィルタリングした場合を示す。
【図１０】非特許文献１に開示された従来の高速フーリエ変換を用いたフィルタの機能構成例を示す図。
【技術分野】
【０００１】
この発明は、例えば収音信号に有限長インパルス応答フィルタを畳み込んで雑音を抑圧するのに用いられる音響信号フィルタとそのフィルタリング方法と、そのプログラムと記録媒体に関する。
【背景技術】
【０００２】
ディジタル信号処理における音響信号のフィルタは、例えば短時間フーリエ変換領域で音響信号に有限長インパルス応答フィルタを畳み込んで実現される。一般に有限長インパルス応答フィルタのフィルタ長が短時間フーリエ変換の分析窓長より十分に短い場合、周波数領域におけるフィルタの畳み込み計算は、周波数信号と周波数フィルタの積として計算することができる。
離散時間のインデックスをｔ、時間領域の音響信号をｘｔ、時間領域の有限長インパルス応答フィルタをｆｔと表わすものとする。時間シフトしながら短時間分析窓を適用することで切り出される時間フレームのインデックスをｎ、周波数信号の周波数ビンに関するインデックスをｍ、Ｘｎ，ｍを音響信号ｘｔの短時間フーリエ変換のフレームｎにおける周波数ビンｍの値、Ｆｍを周波数ビンｍにおける有限長インパルス応答フィルタの短時間フーリエ変換の値とする。音響信号ｘｔに有限長インパルス応答フィルタｆｔを畳み込んで得られるフィルタリングされた信号ｙｔの短時間フーリエ変換値Ｙｎ，ｍは式（１）で近似的に計算することができる。
【０００３】
【数１】
この短時間フーリエ変換と同様の手法である高速フーリエ変換を用いたフィルタ１０００の機能構成を、図１０に示して簡単に説明する（非特許文献１参照）。従来の高速フーリエ変換を用いたフィルタ１０００は、高速フーリエ変換部１０１と、高速フーリエ変換部１０２と、畳み込み演算部１０３と、高速フーリエ逆変換部１０４を備える。高速フーリエ変換部１０１は、音響信号ｘｔを高速フーリエ変換して上記した周波数ビンＸｎ，ｍを生成する。高速フーリエ変換部１０２は、有限長インパルス応答フィルタｆｔを高速フーリエ変換して上記した有限長インパルス応答フィルタの短時間フーリエ変換の値Ｆｍを生成する。畳み込み演算部１０３は、上記した式（１）を演算してフィルタリングされた信号ｙｔの短時間フーリエ変換値Ｙｎ，ｍを生成する。高速フーリエ逆変換部１０４は、短時間フーリエ変換値Ｙｎ，ｍを上記した時間領域の信号ｙｔに変換する。
【非特許文献１】三谷政昭著「ディジタルフィルタデザイン」株式会社昭晃堂、平成６年４月２０日初版７刷発行、ｐ97
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかし、従来の方法では、短時間フーリエ変換領域でフィルタの畳み込みを計算するためには、フィルタ長よりも十分に長い分析窓を使う必要がある。フィルタ長より分析窓長が短い場合は、長い有限長インパルス応答フィルタの畳み込みを周波数領域に閉じて精度良く計算することができなかった。単純に分析窓長を大きくしても、分析の時間分解能が低くなるので音声信号や音楽信号のように短時間で変化する信号を精度良くフィルタリングすることができない。また、いったん時間領域の信号に変換し直して畳み込み計算してから再び周波数領域の信号に戻す方法も考えられるが、計算量が増えてしまう。また、二乗誤差最小基準のような何らかのコスト関数に基づく最適化処理を行なう場合は、異なる領域における演算が必要であると共にコスト関数の最適化の計算が複雑になる課題がある。
【０００５】
この発明は、このような点に鑑みてなされたものであり、分析窓長よりも長い有限長インパルス応答フィルタを時間領域において音響信号に畳み込む演算と等価な演算を、周波数領域において実現する音響信号フィルタとそのフィルタリング方法と、プログラムと記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００６】
この発明の音響信号フィルタは、巡回時間シフト演算部と、畳み込み演算部と、窓関数適用部を具備する。巡回時間シフト演算部は、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数フィルタのフレーム列に、巡回時間シフト演算子を乗算して周波数フィルタを生成する。畳み込み演算部は、離散値化された時間領域音響信号のフレーム長Ｍ＋Ｎ−１の各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する。窓関数適用部は、全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する。そして、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２である。
【発明の効果】
【０００７】
この発明の音響信号フィルタは、分析窓長よりも長い有限長インパルス応答フィルタを時間領域において音響信号に畳み込む演算と等価な演算を、周波数領域において可能にする。このフィルタリング方法は、周波数領域においても畳み込みとして計算できるので、効率的な計算アルゴリズムにすることができる。
【発明を実施するための最良の形態】
【０００８】
以下、この発明の実施の形態を図面を参照して説明する。複数の図面中同一のものには同じ参照符号を付し、説明は繰り返さない。
〔発明の基本的な考え〕
実施例の説明の前に、この発明の音響信号フィルタリング方法の基本的な考えを説明する。まず、時間領域信号に対する畳み込み演算を周波数領域で計算するための関係式を導出する。
ディジタル化した離散時間信号の標本インデックスをｔとし、時間領域の音響信号ｘｔにタップ長Ｋの有限長インパルス応答フィルタｃｔ（ｔ＝１〜Ｋ）を時間領域で畳み込んで得られる出力信号をｙｔとする。時刻ｔ０で始まる出力信号ｙｔの短時間セグメントを時間窓で切り出した信号は、ｚ変換領域において式（２）で表現できる。
【０００９】
【数２】
ここで、ｙ（ｚ）＝ｃ（ｚ）ｘ（ｚ）とし、ＷＮ（・）は時間領域における長さＮの窓関数に相当する関数とした。ＷＮ（ｃ（ｚ））は、ｃ（ｚ）中のｚの−Ｎ+１次から０次の項を取り出し、窓の形に比例して各係数を変更して窓の外の項を除去する。ｚｔ0は、時刻ｔ０で始まるフレームを窓関数の中に移動する時間シフト演算子である。さらに、フィルタｃｔのｔ０タップ目から長さＭのフレームを切り出すことを式（３）で表現する。
【００１０】
【数３】
ここでＷＭＲ（・）は、長さＭの方形窓を表わすものとする。するとｃ（Ｚ）は式（４）で表わせ、式（２）は以下のように書き換えることができる。
【００１１】
【数４】
ここで、Ｋｃはタップ長Ｋを方形窓長Ｍで除したＫｃ＝「Ｋ/Ｍ」とした。「Ｋ/Ｍ」は、
Ｋ/Ｍ以上の最小の整数を表わすものとする。また、Ｍはフレームシフトに対応する。式（５）中で窓関数の引数に含まれる項の内、窓の外の項を除去することで式（６）が導かれる。式（６）の窓関数中の引数の内、式（７）で示す引数は各フレーム同士の時間領域での畳み込みに相当する。
【００１２】
【数５】
つまり、フィルタｃｔのτＭタップ目から長さＭのフレームを切り出し、音響信号ｘｔの時刻ｔ０−Ｍ＋１−τＭから長さＭ＋Ｎ−１のフレームを切り出し、各フレームをｚ領域で乗算したものである。よく知られているように、ｚ領域における乗算は、時間領域における畳み込み演算に相当する。さらに、時間領域における音響信号のフレームとフィルタのフレームの畳み込み演算は、短時間フーリエ変換のポイントをＮｐ、フィルタｃｔのフレーム長をＭ、音響信号ｘｔのフレーム長をＭ＋Ｎ−１としたときに、Ｎｐ≧２Ｍ＋Ｎ−２であるとき、短時間フーリエ変換領域での積で正確に表現することができる。次に、式（６）の窓関数中の引数のうち、ｚＭ−１はｚ領域におけるＭ−１次の時間シフト演算子であり、時間領域においてサンプルの時間インデックスをＭ−１だけ小さくすることに相当する。これは短時間フーリエ変換領域では、Ｍ−１次の巡回時間シフト演算子との積として表わすことができる。本明細書では、これをＧ（Ｍ−１）と書くことにする。さらに、式（６）中のＷＮは、時間領域における窓関数に相当するので、短時間フーリエ変換領域においては、窓関数の短時間フーリエ変換との巡回畳み込みで表現できることが知られている。この巡回畳み込みは、後述するように巡回畳み込み行列との積として表現できる。最後に、式（６）中の総和演算Σは、短時間フーリエ変換領域においても総和演算として表現することができる。また、この総和演算において、τが１増加するに応じて、式（７）の二つのフレームはインデックスがＭずつ増加する。これは、Ｍが上記短時間フーリエ変換の計算におけるフレームシフトに対応していることを意味している。
【００１３】
この発明の音響信号フィルタリング方法は、上記した関係を用いることで時間領域での分析窓長よりも長い有限長インパルス応答フィルタの畳み込み計算を、周波数領域において可能にしたものである。
なお、この発明の説明における短時間フーリエ変換は式（８）と式（９）で定義されるものとする。
【００１４】
【数６】
ここで、Ｔは行列またはベクトルの転置を表わし、Ｘτ、ｍは、時刻ｔτで始まる時間フレームτに関する音響信号ｘｔの短時間フーリエ変換の周波数ビンｍにおける値を表わす。Ｘτは、フレームτにおける全周波数ビンの値を並べたベクトルとする。
【実施例１】
【００１５】
図１にこの発明の音響信号フィルタリング方法を用いた音響信号フィルタ１００の機能構成例を実施例１として示す。その動作フローを図２に示す。音響信号フィルタ１００は、巡回時間シフト演算部１０と、畳み込み演算部１２と、窓関数適用部１４を備える。この例の音響信号フィルタ１００は、例えばＲＯＭ、ＲＡＭ、ＣＰＵ等で構成されるコンピュータに所定のプログラムが読み込まれて、ＣＰＵがそのプログラムを実行することで実現されるものである。以降で説明する音響信号フィルタも同様である。
巡回時間シフト演算部１０は、時間領域の有限長インパルス応答フィルタのフレーム長Ｍ、フレームシフトＭで切り出される各フレームを、それぞれポイント数Ｎｐにより短時間フーリエ変換した周波数フィルタのフレーム列Ｃｎ（ｎ＝１〜Ｋｃ）を入力として、式（１０）に示す演算を行い周波数フィルタＣｎ〜（ｎ＝１〜Ｋｃ）を生成する（ステップＳ１０）。各ｎに対するＣｎは、式（９）と同様、各周波数ビンを要素にもつベクトルである。
【００１６】
【数７】
式（１０）は各周波数ビン毎の周波数フィルタのフレーム列Ｃｎに、ｋ＝Ｍ−１として巡回シフト演算子Ｇ（ｋ）を乗算するものである。ここでｄｉａｇ（Ｘ）は、ベクトルＸを対角要素に持つ対角行列を表わす。巡回時間シフト演算子Ｇ（ｋ）は、周波数領域の表現で式（１１）で表わされるものである。
【００１７】
【数８】
畳み込み演算部１２は、時間領域の音響信号ｘｔのフレーム長Ｍ＋Ｎ−１、フレームシフトＭで切り出される各フレームを、ポイント数Ｎｐにより短時間フーリエ変換した周波数信号のフレーム列を入力として、周波数ビン毎に周波数フィルタＣｎ〜を畳み込み演算して全周波数ビンをまとめたベクトルＹｎを式（１２）の演算により生成する（ステップＳ１２）。
【００１８】
【数９】
窓関数適用部１４は、式（１３）に示す様に全周波数ビンをまとめたベクトルＹｎに窓関数の短時間フーリエ変換値Ｈを一列目に持つエルミートテプリッツ行列ＷＮを乗じて周波数領域出力信号ＷＮＹｎを生成する（ステップＳ１４）。
いま、Ｈを窓関数の短時間フーリエ変換の各周波数ビンを並べたベクトルとし、式（１３）に示すように定義する。
【００１９】
【数１０】
また、エルミートテプリッツ行列ＷＮを式（１４）に示すように定義する。
【００２０】
【数１１】
ここで＊は複素共役を表わす。
【００２１】
【数１２】
周波数領域でエルミートテプリッツ行列ＷＮを乗ずることは、時間領域で対応する時間窓を適用することと等価な演算である。
なお、窓関数は一般に低域通過フィルタの特性を示すので、窓関数の短時間フーリエ変換Ｈは、低域周波数にのみ大きなゲインを持つ。このため、周波数領域における窓関数の演算においては、高域周波数に関する計算を省略しても大きな誤差を生じない。これは、ｍ＝０またはｍ＝Ｎｐ-1近傍の（式（８）を参照）の周波数ビン以外の周波数ビンに対応するＨの値は０と近似しても良いことを意味する。このように近似することで、エルミートテプリッツ行列ＷＮの積はＨの低域周波数近傍の周波数ビンのみに関する畳み込み演算とすることができ、効率的に計算することができる。
【００２２】
以上述べたように、この実施例によれば、時間領域の畳み込み演算を周波数領域における周波数ビン毎の畳み込み演算で表現することができ、時間領域での分析窓長よりも長い有限長インパルス応答フィルタの畳み込み計算を周波数領域において可能にする。この実施例で説明したフィルタリング方法は、周波数領域でも畳み込みとして計算できるので効率的な計算アルゴリズムにすることができる。
なお、周波数領域出力信号ＷＮＹｎを時間領域の信号に変換する場合は、図１に破線で示すように周波数領域出力信号ＷＮＹｎを時間領域の信号に変換する逆短時間フーリエ変換部１９と、時間領域の短時間フレームの時系列をオーバラップ加算するオーバラップ加算部２０とをさらに備え、それぞれで逆短時間フーリエ変換過程（ステップＳ１９）とオーバラップ加算過程（ステップＳ２０）の処理を行なえばよい。これらは、従来技術で容易に実現できるので説明は省略する。
〔変形例１〕
図３に実施例１を変形した音響信号フィルタ３００の機能構成例を示す。その動作フローを図４に示す。音響信号フィルタ３００は、周波数音響信号のフレーム列に巡回時間シフト演算子を乗算して周波数信号を生成するようにし、有限長インパルス応答フィルタ係数の周波数ビン毎に、巡回時間シフト演算子を乗算した周波数信号を畳み込み演算するようにしたものである。
【００２３】
音響信号フィルタ３００は、巡回時間シフト演算部３０と、畳み込み演算部３２と、窓関数適用部１４を備える。窓関数適用部１４は実施例１と同じものである。巡回時間シフト演算部３０は、離散値化された時間領域の音響信号のフレーム長Ｍ、フレームシフトＭで切り出した各フレームをポイント数Ｎｐで短時間フーリエ変換した周波数音響信号のフレーム列Ｘｎに、ｋ＝Ｍ−１として巡回時間シフト演算子Ｇ（ｋ）を乗算して周波数音響信号Ｘｎ〜を生成する（式（１６））（ステップＳ３０）。
【００２４】
【数１３】
畳み込み演算部３２は、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍ＋Ｎ−１、フレームシフトＭで切り出した各フレームに対して、上記ポイント数Ｎｐで短時間フーリエ変換を適用して得られる周波数フィルタＣτ（τ＝１〜Ｋｃ）を、各周波数ビン毎に、上記周波数信号に畳み込み演算して全周波数ビンをまとめたベクトル生成する（式（１７））（ステップＳ３２）。
【００２５】
【数１４】
窓関数適用部１４は、実施例１と同じ動作を行い周波数領域出力信号ＷＮＹｎを生成する。以上述べたように、周波数音響信号のフレーム列に巡回時間シフト演算子を乗算して周波数信号を生成するようにしても、実施例１と同様に時間領域の畳み込み演算を周波数領域における周波数ビン毎の畳み込み演算で表現することができる。
【実施例２】
【００２６】
音響信号の時間窓長（フレーム長）Ｎを固定としたとき、有限長インパルス応答フィルタのフレーム長Ｍ（フレームシフトと同じ値）を大きくするほどフィルタ長Ｋｃを小さくできる。したがって上記した式（１２）の総和における加算の回数を小さくでき、効率的なフィルタリング計算が行える。ただし、ＭがＮよりも大きくなると、全周波数ビンをまとめたベクトルＹｎのフレームシフトが時間窓長Ｎよりも大きくなることに相当し、飛び飛びの時間フレームに対応するＹｎしか求められなくなる。そこで、効率的にフィルタリング計算を行うためには、Ｎ＝Ｍとすることが一つの解決策となる。Ｎ＝Ｍとすることで、実施例１で示した式（６）は式（１８）に書き換えることができる。
【００２７】
【数１５】
ｆｔは、同じ長さＮの時間フレーム同士の畳み込みの和である。短時間フーリエ変換のポ
イント数Ｎｐが、Ｎｐ≧２Ｎ−１のとき、短時間フーリエ変換の積の和として表現できる。
また、ｚＮは、Ｎ次のオーダの時間シフトオペレータであり、Ｎｐ≧２Ｎの短時間フーリエ変換を用いる場合、時間シフト計算が適切に行えることを意味する。
時間領域の窓関数が、周波数領域では巡回畳み込みで表現されることを考慮すると、式
（１８）は式（２０）に書き換えることができる。
【００２８】
【数１６】
音響信号フィルタ５００は式（２０）に示す動作を行なうようにしたものである。音響
信号フィルタ５００の機能構成例を図５にその動作フローを図６に示す。音響信号フィル
タ５００は、巡回シフト演算部５０と、畳み込み演算部５２と、窓関数適用部１４を備え
る。窓関数適用部１４は実施例１と同じものである。
巡回時間シフト演算部５０は、離散値化された時間領域の音響信号のフレーム長Ｎ、フレームシフトＭで切り出した各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数信号Ｘｎ〜を生成する（式（２１））（ステップＳ５０）。
【００２９】
畳み込み演算部５２は、時間領域の有限長インパルス応答フィルタ係数からフレーム長Ｍ、フレームシフトＭで切り出した各フレームに対して、ポイント数Ｎｐで短時間フーリエ変換を適用して得られる周波数フィルタを、周波数ビン毎に、連結周波数信号に畳み込み演算して全周波数ビンをまとめたベクトルＹｎを生成する（式（２０））（ステップＳ５２）。
このように構成しても実施例１と同様に時間領域の畳み込み演算を周波数領域における周波数ビン毎の畳み込み演算で表現することができる。Ｎ＝Ｍとすることで、音響信号フィルタ５００の入力信号のフレーム長と出力信号である周波数出力信号のフレーム長が同一になるので計算上の都合が良い。
【００３０】
なお、実施例２ではフレームシフト幅Ｍと窓長Ｎとが一致するため、隣接する時間窓がオーバラップしない。オーバラップしない信号を、そのまま時間領域の信号に変換すると雑音が増加する場合がある。それを防止するためには、例えば、時間領域信号から周波数信号の時系列を生成する際に適用する短時間フーリエ変換の時間シフトを時間窓長の半分にする。そして、奇数番目及び偶数番目のフレームインデックスを持つ二種類のフレーム系列を取り出し、各系列に別々に周波数フィルタを畳み込んだ後に一つの系列に戻すことで、隣接する時間窓が２分の１ずつオーバラップするフレーム系列を作ることができる。
〔変形例２〕
図７に実施例２を変形した音響信号フィルタ７００の機能構成例を示す。その動作フローを図８に示す。音響信号フィルタ７００は、実施例１に対する変形例１と同じように、巡回時間シフト演算子を乗算する信号系列を変更したものである。音響信号フィルタ７００は、巡回シフト演算部７０と、畳み込み演算部７２と、窓関数適用部１４を備える。窓関数適用部１４は実施例１と同じものである。
【００３１】
巡回時間シフト演算部７０は、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｎ、フレームシフトＭで切り出した各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタの各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数フィルタＣｎ〜を生成する（式（２２））（ステップＳ７０）。
畳み込み演算部７２は、時間領域の音響信号のフレーム長Ｎの時系列からフレーム長Ｍ、フレームシフトＭで切り出した各フレームに対して、ポイント数Ｎｐで短時間フーリエ変換を適用して得られる周波数音響信号を、周波数ビン毎に、連結周波数フィルタＣｎ〜に畳み込み演算して全周波数ビンをまとめたベクトルＹｎを生成する（式（２３））（ステップＳ７２）。
【００３２】
【数１７】
このように構成しても実施例２と同様に、周波数領域における周波数ビン毎の畳み込み演算を効率良くおこなうことができる。
【００３３】
〔シミュレーション結果〕
この発明の音響信号フィルタの動作を確認するシミュレーションを行った。シミュレー
ションは、長さ２０４８サンプルの時系列データを入力とし、長さ１０２４サンプルの有
限長インパルス応答フィルタを−１〜１の範囲乱数で生成する条件で行った。その結果を
図９に示す。図９（ａ）は時間領域の畳み込み演算でフィルタリングした場合を示す。図
９（ｂ）は有限長インパルス応答フィルタの時間窓長Ｍ＝６４、音響信号の時間窓長Ｎ＝
２５６として実施例１の方法でフィルタリングした場合を示す。図９（ｃ）はＭ＝２５６、Ｎ＝２５６として実施例２の方法でフィルタリングした場合を示す。それぞれの横軸は、時系列データのサンプル番号、縦軸は振幅である。
【００３４】
図９（ｂ）と（ｃ）の振幅０に引かれた破線は、時間領域畳み込み演算のフィルタリン
グ結果との差分を示す。このように実施例１と２のどちらのフィルタリング方法でも時間
領域畳み込み演算でフィルタリングしたのと同じ結果が得られた。
このようにこの発明のフィルタリング方法によれば、分析窓長よりも長い有限長インパ
ルス応答フィルタを時間領域において音響信号に畳み込む演算と等価な演算を、周波数領
域において可能にする。
なお、この発明の音響信号フィルタリング方法及び音響信号フィルタは上述の実施形態に限定されるものではなく、この発明の趣旨を逸脱しない範囲で適宜変更が可能であり、この発明の技術思想を、音響信号以外の２系統のディジタル信号系列に対して適用することも可能である。また、上記方法及び装置において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【００３５】
また、上記装置における処理手段をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、各装置における処理手段がコンピュータ上で実現される。
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）/ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【００３６】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記録装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
また、各手段は、コンピュータ上で所定のプログラムを実行させることにより構成することにしてもよいし、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図面の簡単な説明】
【００３７】
【図１】この発明の音響信号フィルタ１００の機能構成例を示す図。
【図２】音響信号フィルタ１００の動作フローを示す図。
【図３】この発明の音響信号フィルタ３００の機能構成例を示す図。
【図４】音響信号フィルタ３００の動作フローを示す図。
【図５】この発明の音響信号フィルタ５００の機能構成例を示す図。
【図６】音響信号フィルタ５００の動作フローを示す図。
【図７】この発明の音響信号フィルタ７００の機能構成例を示す図。
【図８】音響信号フィルタ７００の動作フローを示す図。
【図９】シミュレーション結果を示す図であり、（ａ）は時間領域の畳み込み演算でフィルタリングした場合、（ｂ）は有限長インパルス応答フィルタの時間窓長Ｍ＝６４、音響信号の時間窓長Ｎ＝２５６として実施例１の方法でフィルタリングした場合、（ｃ）はＭ＝２５６、Ｎ＝２５６として実施例２の方法でフィルタリングした場合を示す。
【図１０】非特許文献１に開示された従来の高速フーリエ変換を用いたフィルタの機能構成例を示す図。
【特許請求の範囲】
【請求項１】
時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームを、ポイント数Ｎｐで短時間フーリエ変換したフレーム列に、巡回時間シフト演算子を乗算して周波数フィルタを生成する巡回時間シフト演算部と、
離散値化された時間領域の音響信号のフレーム長Ｍ＋Ｎ−１の各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、上記周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算部と、
上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用部と、
を備え、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２であることを特徴とする音響信号フィルタ。
【請求項２】
離散値化された時間領域の音響信号のフレーム長Ｍの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号のフレーム列に、巡回時間シフト演算子を乗算して周波数信号を生成する巡回時間シフト演算部と、
時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍ＋Ｎ−１の各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数ビン毎に、上記周波数信号を畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算部と、
上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用部と、
を備え、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２であることを特徴とする音響信号フィルタ。
【請求項３】
離散値化された時間領域の音響信号のフレーム長Ｎの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数信号を生成する巡回時間シフト演算部と、
時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタの周波数ビン毎に、上記連結周波数信号を畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算部と、
上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用部と、
を具備し、上記フレーム長Ｎは上記フレーム長Ｍと等しく、且つ上記ポイント数ＮｐがＮｐ≧２Ｎであることを特徴とする音響信号フィルタ。
【請求項４】
時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｎの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタの各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数フィルタを生成する巡回時間シフト演算部と、
離散値化された時間領域の音響信号の上記フレーム長Ｍの各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、上記連結周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算部と、
上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用部と、
を具備し、上記フレーム長Ｎは上記フレーム長Ｍと等しく、且つ上記ポイント数ＮｐがＮｐ≧２Ｎであることを特徴とする音響信号フィルタ。
【請求項５】
巡回時間シフト演算部が、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームを、ポイント数Ｎｐで短時間フーリエ変換したフレーム列に、巡回時間シフト演算子を乗算して周波数フィルタを生成する巡回時間シフト演算過程と、
畳み込み演算部が、離散値化された時間領域の音響信号のフレーム長Ｍ＋Ｎ−１の各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、上記周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算過程と、
窓関数適用部が、上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用過程と、
を備え、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２であることを特徴とする音響信号フィルタリング方法。
【請求項６】
巡回時間シフト演算部が、離散値化された時間領域の音響信号のフレーム長Ｍの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号のフレーム列に、巡回時間シフト演算子を乗算して周波数信号を生成する巡回時間シフト演算過程と、
畳み込み演算部が、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍ＋Ｎ−１の各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数ビン毎に、上記周波数信号を畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算過程と、
窓関数適用部が、上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用過程と、
を備え、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２であることを特徴とする音響信号フィルタリング方法。
【請求項７】
巡回時間シフト演算部が、離散値化された時間領域の音響信号のフレーム長Ｎの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号のフレーム列の各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数信号を生成する巡回時間シフト演算過程と、
窓関数適用部が、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタのフレーム列の周波数ビン毎に、上記連結周波数信号を畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算過程と、
窓関数適用部が、上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用過程と、を備え、上記フレーム長Ｎは上記フレーム長Ｍと等しく、且つ上記ポイント数ＮｐがＮｐ≧２Ｎであることを特徴とする音響信号フィルタリング方法。
【請求項８】
巡回時間シフト演算部が、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｎの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタの各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数フィルタを生成する巡回時間シフト演算過程と、
畳み込み演算部が、離散値化された時間領域の音響信号の上記フレーム長Ｎの各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、上記連結周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算過程と、
窓関数適用部が、上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用過程と、
を具備し、上記フレーム長Ｎは上記フレーム長Ｍと等しく、且つ上記ポイント数ＮｐがＮｐ≧２Ｎであることを特徴とする音響信号フィルタリング方法。
【請求項９】
請求項１乃至４の何れかに記載した音響信号フィルタとしてコンピュータを機能させるための装置プログラム。
【請求項１０】
請求項９に記載した何れかの装置プログラムを記録したコンピュータで読み取り可能な記録媒体。
【請求項１】
時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームを、ポイント数Ｎｐで短時間フーリエ変換したフレーム列に、巡回時間シフト演算子を乗算して周波数フィルタを生成する巡回時間シフト演算部と、
離散値化された時間領域の音響信号のフレーム長Ｍ＋Ｎ−１の各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、上記周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算部と、
上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用部と、
を備え、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２であることを特徴とする音響信号フィルタ。
【請求項２】
離散値化された時間領域の音響信号のフレーム長Ｍの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号のフレーム列に、巡回時間シフト演算子を乗算して周波数信号を生成する巡回時間シフト演算部と、
時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍ＋Ｎ−１の各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数ビン毎に、上記周波数信号を畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算部と、
上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用部と、
を備え、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２であることを特徴とする音響信号フィルタ。
【請求項３】
離散値化された時間領域の音響信号のフレーム長Ｎの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数信号を生成する巡回時間シフト演算部と、
時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタの周波数ビン毎に、上記連結周波数信号を畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算部と、
上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用部と、
を具備し、上記フレーム長Ｎは上記フレーム長Ｍと等しく、且つ上記ポイント数ＮｐがＮｐ≧２Ｎであることを特徴とする音響信号フィルタ。
【請求項４】
時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｎの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタの各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数フィルタを生成する巡回時間シフト演算部と、
離散値化された時間領域の音響信号の上記フレーム長Ｍの各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、上記連結周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算部と、
上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用部と、
を具備し、上記フレーム長Ｎは上記フレーム長Ｍと等しく、且つ上記ポイント数ＮｐがＮｐ≧２Ｎであることを特徴とする音響信号フィルタ。
【請求項５】
巡回時間シフト演算部が、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームを、ポイント数Ｎｐで短時間フーリエ変換したフレーム列に、巡回時間シフト演算子を乗算して周波数フィルタを生成する巡回時間シフト演算過程と、
畳み込み演算部が、離散値化された時間領域の音響信号のフレーム長Ｍ＋Ｎ−１の各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、上記周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算過程と、
窓関数適用部が、上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用過程と、
を備え、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２であることを特徴とする音響信号フィルタリング方法。
【請求項６】
巡回時間シフト演算部が、離散値化された時間領域の音響信号のフレーム長Ｍの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号のフレーム列に、巡回時間シフト演算子を乗算して周波数信号を生成する巡回時間シフト演算過程と、
畳み込み演算部が、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍ＋Ｎ−１の各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数ビン毎に、上記周波数信号を畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算過程と、
窓関数適用部が、上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用過程と、
を備え、上記ポイント数ＮｐがＮｐ≧２Ｍ＋Ｎ−２であることを特徴とする音響信号フィルタリング方法。
【請求項７】
巡回時間シフト演算部が、離散値化された時間領域の音響信号のフレーム長Ｎの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号のフレーム列の各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数信号を生成する巡回時間シフト演算過程と、
窓関数適用部が、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｍの各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタのフレーム列の周波数ビン毎に、上記連結周波数信号を畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算過程と、
窓関数適用部が、上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用過程と、を備え、上記フレーム長Ｎは上記フレーム長Ｍと等しく、且つ上記ポイント数ＮｐがＮｐ≧２Ｎであることを特徴とする音響信号フィルタリング方法。
【請求項８】
巡回時間シフト演算部が、時間領域の有限長インパルス応答フィルタ係数のフレーム長Ｎの各フレームを、ポイント数Ｎｐで短時間フーリエ変換した周波数フィルタの各周波数ビンに、上記各周波数ビンの１個前の周波数ビンに巡回時間シフト演算子を乗算した値を加算して連結周波数フィルタを生成する巡回時間シフト演算過程と、
畳み込み演算部が、離散値化された時間領域の音響信号の上記フレーム長Ｎの各フレームを、上記ポイント数Ｎｐで短時間フーリエ変換した周波数音響信号の周波数ビン毎に、上記連結周波数フィルタを畳み込み演算して全周波数ビンをまとめたベクトルを生成する畳み込み演算過程と、
窓関数適用部が、上記全周波数ビンをまとめたベクトルに窓関数の短時間フーリエ変換値を一列目に持つ巡回畳み込み行列を乗じて周波数領域出力信号を生成する窓関数適用過程と、
を具備し、上記フレーム長Ｎは上記フレーム長Ｍと等しく、且つ上記ポイント数ＮｐがＮｐ≧２Ｎであることを特徴とする音響信号フィルタリング方法。
【請求項９】
請求項１乃至４の何れかに記載した音響信号フィルタとしてコンピュータを機能させるための装置プログラム。
【請求項１０】
請求項９に記載した何れかの装置プログラムを記録したコンピュータで読み取り可能な記録媒体。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【公開番号】特開２００９−２１０６４２（Ｐ２００９−２１０６４２Ａ）
【公開日】平成２１年９月１７日（２００９．９．１７）
【国際特許分類】
【出願番号】特願２００８−５１１０７（Ｐ２００８−５１１０７）
【出願日】平成２０年２月２９日（２００８．２．２９）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【公開日】平成２１年９月１７日（２００９．９．１７）
【国際特許分類】
【出願日】平成２０年２月２９日（２００８．２．２９）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
[ Back to top ]