JP5553334B2 - 正弦波パラメータ推定方法 - Google Patents

Description

本発明は、音響センシング、音響信号処理分野に関し、特に信号波形中に含まれる正弦波パラメータを推定する方法に関する。

さらに本発明は、音響や光など、広く波動場のセンシングや信号処理、通信やデータ記録における復調や復号化の技術においては、正弦波パラメータの推定は極めて重要な技術事項であるので、それらの分野にも強く関連する。

音声信号などの信号波形解析においては、その信号に含まれる正弦波のパラメータを推定することは基本的、且つ重要な事項である。

したがって、従来から信号中の正弦波パラメータ推定技術として種々のものが提案されている。広く利用されている技術は、ＦＦＴ（高速フーリエ変換）を利用したものであり、特にその得られたＦＦＴスペクトルを２次補間して正弦波パラメータを推定する方法は簡単で精度も高いので広く利用されている。

先行技術文献
正弦波パラメータ推定と関連性があると思われるいくつかの先行技術を挙げる。

下記特許文献１には、正弦波の振幅と周波数とを推定する方法等が開示されている。ここに開示されている技術は、サンプリング信号に基づき、これらの信号間を補間すること等によって正弦波の推定を行うことを特徴とする。

また、下記特許文献２には、未知の信号を推定するスペクトル分析装置が開示されている。ここに開示されているスペクトル分析装置によれば、入力信号の周波数や位相が推定できるとされている。

また、下記特許文献３には、入力音声の正弦波を抽出し、その正弦波を利用して音声合成する技術が開示されている。

また、下記特許文献４には、ノイズが音声信号に混入した場合、そのノイズの期間における正弦波を推定し、入力信号を補間することによってノイズキャンセルを行う装置が開示されている。

また、下記非特許文献１においては、ＦＦＴの２次補間による正弦波パラメータの推定の基準に関する研究が開示されている。特に、その設計パラメータ、窓関数や、窓の長さ、零詰めを行う際の零詰めの量等について研究がなされている。

特開平０７−０８３９６４号公報 特開２０００−２５８４７８号公報 特開２００４−１０９８０９号公報 特開２００６−１３５８０６号公報 安部, スミス, "ＦＦＴの２次補間に基づく正弦波パラメータ推定法の設計基準",信学技報, vol.104 No.304, pp.7-12, Nov. 2004.

このように、未知の信号に対して正弦波のパラメータを調べる従来の手法は、補間等のテクニックを用いた近似であり、精度の向上にも一定の限界があると考えられる。

そこで、近似ではなく、直接的に正弦波のパラメータ（振幅、位相）を推定できる手法が望ましいことは言うまでもない。

理論的には、ナイキストの定理から、サンプリング周期が十分に高ければそのサンプリングした信号の値から（量子化誤差を除けば）元の正弦波信号は完全に再現できるはずであり、そのような手法が望まれている。

本願発明は、このような背景に鑑みなされたものであり、補間等の近似ではなく、入力信号を計測した値から、直接的にその入力信号中に含まれる正弦波のパラメータを推定する技術を実現することを目的とする。

本願発明は、上記課題を解決するために、微分方程式の有限区間荷重積分に基づく正弦波パラメータの直接推定を行うことをその解決原理とするものである。

また、本特許においては、正弦波のパラメータとは、周波数と振幅と位相を意味する。なお、正弦波のパラメータという場合、周波数と振幅と位相とに加えて直流バイアスを意味する場合もある。

本特許においては、この直流バイアスを直接パラメータとして扱ってはいないが、周波数が直流という特別の場合の振幅として推定することができる。

具体的には、上記課題を解決するため以下のような手段を採用している。

（１）本発明は、上記課題を解決するために、解析の対象である信号中の正弦波のパラメータを推定する方法において、前記信号の微分を含む支配方程式を立てる工程と、前記微分を含む支配方程式に任意の荷重関数を導入し、有限区間の積分形式に置き換える工程と、前記積分形式に変換した方程式を解く工程と、
を含み、
前記積分形式の支配方程式を解く工程は、
前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程
を含むことを特徴とする正弦波パラメータ推定方法である。

（２）本発明は、解析の対象である信号中の正弦波のパラメータを推定する方法において、任意の荷重関数を導入した有限区間の積分形式で、前記信号の微分を含む支配方程式を構成する工程と、前記積分形式の支配方程式を解く工程と、
を含み、
前記積分形式の支配方程式を解く工程は、
前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程
を含むことを特徴とする正弦波パラメータ推定方法である。

（３）また、本発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記荷重関数は、前記有限期間の整数分の１の周期の三角関数であることを特徴とする正弦波パラメータ推定方法である。

整数分の１にとることによって、後述するように、未知数や計数の個数を近似無しで減じ、かつ、必要な個数は残すことが可能である。

（４）また、本発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記積分形式の支配方程式を解く工程は、異なる３個の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程と、を含むことを特徴とする正弦波パラメータ推定方法である。

（５）また、本発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記積分形式の支配方程式を解く工程は、異なる２個の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程と、を含むことを特徴とする正弦波パラメータ推定方法である。

（６）また、本発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記積分形式の支配方程式を解く工程は、３個以上の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程と、を含むことを特徴とする正弦波パラメータ推定方法である。

（７）また、他の関連発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記積分形式の支配方程式を解く工程は、１個の荷重周波数と窓関数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程と、を含むことを特徴とする正弦波パラメータ推定方法である。

（８）また、他の関連発明は、上記（４）〜（７）のいずれかに記載の正弦波パラメータ推定方法において、前記パラメータを推定する工程は、前記推定された周波数と、積分境界値項と、その推定された周波数の信号の振幅と位相とを正弦波のパラメータとして求める工程、を含むことを特徴とする正弦波パラメータ推定方法。ここで、積分境界値項とは積分区間の信号やその微分の両端の値の差をいう。

（９）また、他の関連発明は、上記（４）〜（７）のいずれかに記載の正弦波パラメータ推定方法において、前記パラメータを推定する工程は、前記推定された周波数に関して、位相と振幅とをパラメータとして含む正弦波モデルを構築し、この正弦波モデルを解析対象となる解析信号にフィッティングさせることによって、最もフィットする正弦波モデルの振幅と位相とを求める工程、を含むことを特徴とする正弦波パラメータ推定方法である。

（１０）また、他の関連発明は、上記（１）〜（９）のいずれかに記載の正弦波パラメータ推定方法において、前記積分形式に置き換える工程は、前記積分の積分境界を、前記信号を標本化するタイミングの中間に設定することを特徴とする正弦波パラメータ推定方法である。

参考発明の手段
以下、本件発明と関連する参考発明の手段を説明する。

上で述べた正弦波パラメータ推定方法は、以下のようなハードウェアを用いて実施することが好適である。

（１１）本参考発明は、上記（１）〜（１０）のいずれかに記載の正弦波パラメータ推定方法の前記荷重積分を支援する装置において、前記解析信号に前記荷重を乗算するＮ個の乗算手段と、前記Ｎ個の乗算手段毎にそれぞれ設けられ、前記乗算手段の各出力信号をＮ周期遅延させるＮ個の遅延手段と、前記Ｎ個の遅延手段毎にそれぞれ設けられ、遅延後の出力信号と、遅延前の信号との差を求め、その差を累積値に累積するＮ個の加算手段と、前記加算手段毎にそれぞれ設けられ、前記加算手段の出力信号を累積するＮ個の累積手段と、を含むことを特徴とする正弦波パラメータ推定方法支援装置である。ここで、前記Ｎは正の整数である。

（１２）さらに、上で述べた正弦波パラメータ推定方法を実施する為に、任意の荷重関数を導入した有限区間の積分形式で、前記信号の支配方程式を構成する手段と、前記方程式を解く手段と、を含むことを特徴とする正弦波パラメータ推定支援装置を用いることが好ましい。

以上述べたように、本発明によれば、解析信号中に含まれる正弦波のパラメータを直接推定することができる。

本実施例における反復法２の処理の流れのフローチャートである。 本実施例におけるシミュレーションした周波数の推定結果を表すグラフである。 本実施例におけるｓｉｎｃ関数の積分関数（正弦積分）と、複素正弦波荷重のｓｉｎｃ積分関数（ω＝２π／４）のグラフである。 本実施例における積分範囲を標本化メッシュに対して標本点の中間位置までシフトした場合の荷重分布と誤差分布のグラフである。 本実施例における移動有限フーリエ変換を実行する構成のブロック図である。 本実施例において、方程式の他の解法の窓関数と荷重関数の例を示すグラフである。

符号の説明

１０ 乗算手段
２０ 遅延部
２２ 単位遅延手段
３０ 累積部
３２ 加算手段

第１ 解析信号の周波数等の推定
ａ）支配方程式
解析の対象となる信号を、解析信号と呼ぶ。まず、支配方程式について説明する。

解析信号に対して下記式（１）

の一般解は、

となる。ここで、

であるから、その結果、下記式（２）

が得られる。すなわち、上で述べた式（１）（数１）が支配方程式である。複素正弦波の周波数は、この式中のａの虚部に相当し、その振幅変化率はａの実部に表現される。この支配方程式が、請求の範囲の「微分を含む支配方程式」の好適な一例に相当する。

ｂ）荷重積分観測方程式
観測期間を［−Ｔ／２，Ｔ／２］（すなわり、総時間はＴ）とし、その期間中、信号周波数は一定であると仮定する。すなわり、上記ａ）で述べた微分方程式が一様に成立していると仮定する。この条件は、任意の荷重関数

を導入すると、下記の同値関係式（３）によって、積分形式に置き換えることができる。

この支配方程式が、請求の範囲の「積分形式に変換した支配方程式」の好適な一例に相当する。この式（３）において、ｗ（ｔ）は、具体的には、区間［−Ｔ／２，Ｔ／２］で完備な関数系

を荷重関数として用いればよい。

この式（３）の下で微分を消去するために、式（３）の後段の積分を実行すると、

となる。ここで、

と選ぶと、これらは区間［−Ｔ／２，Ｔ／２］で完備な直交系をなし、

のように、荷重の微分による荷重積分が原関数の荷重積分に比例する。さらに、

のように、全ての

に関する積分境界値項が符号を除いて同一となる。ここで、積分境界値項とは積分区間の信号やその微分の両端の値の差を言う。以下、同様である。さてこの結果、

及び、荷重積分を、下記式（４）に示すように

と置くことによって、題意の荷重積分式は

を未知変数とする線形方程式（下記式（５））

に帰着させることができる。ａを任意の複素数としても未知数は４個であり、方程式は上記の実部と虚部の２個である。したがって、上記の線形方程式は、任意の２個の周波数で連立することによって解くことが可能になる。

ｃ）ＤＣ成分との連立法
さて、具体的にω＝０の式と上式とを連立させると、

より、下記の式（６）

と解かれる。

ここで注意すべきことは、ｆ（ｔ）は解析の対象である解析信号だということである。したがって、ｆ（ｔ）は虚部を有し、

にはこのｆ（ｔ）のヒルベルト変換のフーリエ変換が加わり、

も虚部を有する。解析信号化が上記の荷重積分で同時になされると考えることも可能ではあるが、それは何を解析信号であると認識するかの問題でもある。

なお、ヒルベルト変換のフーリエ変換が加わっても

は２倍されるだけであろうかという疑問が生じるかもしれないが、正しいヒルベルト変換ではないので２倍にはならない点にも注意すべきである。

ｄ）近接２周波数成分を用いた連立法
今、利用する２個の周波数をω_１、ω_２とおき、それらの周波数におけるｓ_ｔをｓ_１，ｓ_２とおくと、連立方程式は

となる。したがって、

となり、よって、下記（７）式、（８）式

を得る。

ｅ）振幅と位相の推定
さて、推定された周波数ａとして

とおくと（ここでＡは複素振幅を表す）、

より下記（９）式

が求められる。ここで、ａＴ＝２ｎπのときはｓｉｎ（ａＴ／２）＝０となって適用不可能となるが、この場合は、単一周波数スペクトルとなって、振幅も位相もこのスペクトル値から決定することができる。

第２ 実信号の周波数等の推定
ａ）支配方程式
解析の対象となる実信号に関して、まず上の記載と同様に、支配方程式についての説明を行う。実信号に対して下記式（１０）

の一般解は、周知のように、下記式（１１）

となる。ここで、ＡやＢは符号も含めて任意である。すなわち、上で述べた式（１０）（数２８）が支配方程式であり、上記式（１１）の一般解の形の意味でのａの推定が周波数を推定することに相当する。この支配方程式が、請求の範囲の「微分を含む支配方程式」の好適な一例に相当する。ここでａの符号は決定しない、又は意味を有しないと言い換えてもよい。

ｂ）荷重積分観測方程式
観測期間を、上の記載と同様に、［−Ｔ／２，Ｔ／２］（すなわち、総時間はＴ）とし、その期間中、信号周波数は一定、すなわち、上記ａ）で述べた微分方程式（支配方程式）が一様に成立していると仮定する。この条件は、上述したのと同様に、任意の荷重関数

を導入すると、下記の同値関係式（１２）によって、積分形式に置き換えることができる。

この支配方程式が、請求の範囲の「積分形式に変換した支配方程式」の好適な一例に相当する。この式（１２）（数３１）において、上述した荷重関数は、具体的には、区間［−Ｔ／２，Ｔ／２］で完備な関数系

を用いればよい。一般に、

であるので、前節で述べたように、

とおき、

とおくと、前節で述べたのと同様の論法によって、同値条件を構成する周波数ごとの方程式として、下記式（１３）

が得られる。

なお、支配方程式を荷重関数を用いて積分形式にすることは上で述べたような機械的な作業である。特に、荷重関数の周期を積分区間の整数分の一にとっていることも、単純な作業の実現に寄与している。そのため、コンピュータ等の上で実施することが好適である。このようにして構成した式に、実際に計測したデータを導入すれば、正弦波のパラメータを変数として含む方程式が完成する。

この処理に際して、まず支配方程式を決めて、それから荷重関数を導入して変形する処理をコンピュータの上で行わせるの好適であるし、また、積分形式の支配方程式の形を予め準備しておき、その係数を、荷重関数に応じて設定して方程式を構成することも好ましい。

請求の範囲において方程式を「立てる」としている動作は、上述のように方程式を構成する処理を意味し、コンピュータで実行することが迅速な処理の為には好ましい。

また、荷重関数の周期は、対象となる信号（正弦波）の周期と区別する為に、特に、荷重周期と呼ぶ場合がある。同様に、荷重関数の周波数は、対象となる信号の周波数と区別する為に、特に、荷重周波数と呼ぶ場合がある。

ｃ）ＤＣ成分を含む連立法
さて、ω＝０に対する下記の関係を用いると、

積分境界値が１個消去され、関係式は

となる。この方程式は複素変数の方程式であり、

の実部に注意し、上記方程式の実部のみを取り出すと、

が得られる。よって、以下の式（１４）のように、

と求められる。また、積分境界値は、

のように書くことができる（式（１５）及び式（１６））。

ｄ）ＤＣ成分を含む３周波数を用いた連立法
荷重積分観測方程式を解くもう一つの連立法は、上記ＤＣに加えてさらに２個の周波数ω_１、ω_２を連立する方法である。それぞれの周波数におけるｓ_ｔをｓ_１、ｓ_２とおくと、これらを連立することによって

であり、よってこの前段の第１式にｓ_２ω_２を乗じ、後段の第２式にｓ_１ω_１を乗じて差し引くことによって、

のように、残りの積分境界値も消去され、よって、下記式（１７）

が得られる。対称性が良く、従来の周波数成分の補間のほうほうと対応が取りやすいという利点を有する。

ｅ）近接３周波数成分を用いた連立法
３個の周波数をω_１、ω_２、ω_３とおき、そのときのｓ_ｔをそれぞれｓ_１、ｓ_２、ｓ_３とおくと、これらを連立することによって、

を得る。この方程式の解は下記式（１８）、式（１９）、式（２０）のようになる。

ただし、

である（式（２１））。なお、ここで、２×２と３×３の行列の逆行列の公式を用いている（下記数４９参照）。

また、念のため、ω_３＝０、ｓ_３＝１の場合について求めてみると

となり、前節の結果と一致することが理解されよう。

ｆ−１）振幅と位相の推定
さて、推定された周波数と、同時に推定された積分境界値項

から、振幅・位相を求めることが考えられる。すなわち、ａ＞０を推定された周波数として、

とおく。すると、

であるので、

すなわち、下記式（２２）

を求めることができ、この結果から、振幅と位相が求められる。なお、数５３の導出には、下記三角関数の和と積の公式が用いられている。

さて、ａＴ＝２ｎπのときには、ｓｉｎ（ａＴ／２）＝０となってこの方法を適用することが困難であるが、この場合は、ωＴ＝２ｎπの中の単一周波数のみにスペクトルが立つ（表れる）ので、周波数も位相もこのスペクトルから決定可能である。

したがって、アルゴリズムとしては、切り換えるための判断ステップ、又は、切り換えまで統合した推定アルゴリズムを用いることが好ましい。そのようなスペクトルのみが立つ場合を検出することはスペクトルを検査すれば比較的容易に判断できるので、その判断に基づきアルゴリズムを切り換えるのは当業者であれば容易に実施することができる。

ｆ−２ａ）振幅と位相の他の推定方式
また、他の方式として、推定に用いた周波数成分

へ正弦波モデル

をフィッティングする方法もある。この方法は、上述した単一の周波数のみにスペクトルが立つ場合でも、条件判断をしてアルゴリズムを切り換える必要がないので便利である。また、全ての周波数成分の影響を受ける積分境界値を用いないため、複数の周波数が混合する信号を取り扱う場合にも高い精度を維持できると考えられる。

すなわち、標記モデルの最小二乗評価関数は

を展開すると、

であり、Ａで微分して零とおくと、

よって、

となる。この式では、Ａは常に実数であるが、正の値だけでなく負にもなりえる点には注意すべきである。次に、φで微分して零とおくと、

すなわち、

すなわち、

を得る。この式では、位相が［−π／２，π／２］の範囲でしか求まらないが、上記Ａの符号を正に固定し、符号を位相φに含ませると、全範囲で位相を推定することが可能である。

ｆ−２ｂ）上記ｆ−２の別解
複素振幅をＡで表し、最小二乗評価関数

をＡ、Ａ^＊で微分して零とおくと、

よって、

となる。

ｇ）複素正弦波の同時推定法（複数推定法）
観測信号（「解析信号」、「実信号」とも呼び、解析の対象となる信号を意味する）に複数の正弦波が含まれる場合について説明する。この観測信号に含まれる正弦波の個数Ｎは既知であるとする。このＮ個は予め知ることができる場合も多い。例えば、離散スペクトル分布のピークの個数から予めＮを求めることができる。

このＮの算出もコンピュータ上で行うことが好ましい。スペクトル分布からピークを算出し、その個数をカウントすることはコンピュータによる一般的な信号処理として知られているので、そのような技術を用いればよい。

さて、このＮ個の周波数を、

とおく。この記法は、本節ｇ）においてのみ採用する記法であり、他の節では他の記法を採用する場合がある。また、正弦波を

とおくと、

と記述できる。これらによって表される各信号が満たすべき方程式は、

である。よって、これに対応する荷重積分式は、

である。この方程式は複素方程式であり、１個のωに対して２個の方程式を与える。一方、観測式としては全ての正弦波の混合信号に対する荷重積分

しか得られない。この式も複素方程式であり、同様に２個以上のωの組に対して２個の方程式を与える。これに対して、未知数は、１個の正弦波に対して、周波数

の１個と、実数である２種類の積分境界値である。したがって、この境界積分値の１＋１＝２個を足し、計３個であり、これに周波数Ｎ個を乗じて、合計３Ｎ個の未知数がある。ここで、未知数を周波数と振幅と位相で考えた場合も同様の結果になる。

さて、これらの未知数と荷重積分は一意に結ばれるため、ここで、２個の方程式が立てられる。

よって、推定に用いるωの個数をＭとすると、これらを未知数として解が得られるための必要条件は、下記式（２３）

となる。Ｎ＝１の場合は、この式（２３）から、Ｍ＞３／２と求まる。すなわち、２個の周波数が必要で、前節の２周波数連立法の結果と一致することが理解されよう。

なお、これらの方程式の構築は、コンピュータ上で実行することが好ましい。所定のωから方程式を得ることは上述のように機械的な作業であるので、汎用的なコンピュータ上でも容易に実行することが可能である。

本特許において、方程式を「立てる」と述べている部分があるが、本特許で「立てる」と述べている作業は、上述のごとく機械的に、決まった形の方程式を「準備する」作業であるので、汎用的なコンピュータで十分に実行することが可能である。

さて、標本値からのフーリエ係数の計算によってわずかに加わる誤差を無視すると、有限フーリエ変換は離散フーリエ変換で高速に実行することができる。このとき、例えば、２５６点のＦＦＴを使用することを前提とすれば、Ｍ（推定に用いるωの個数）＝１２７とすると（ここで、説明を簡単にするため、実の直流成分と最高周波数成分とは除く）、Ｎ＜（２Ｍ＋２）／３＝２５７／３＝８５．３３が得られる。したがって、計算上、８５個までの正弦波のパラメータを推定することができる。同様に、１２８点のＦＦＴでは４２個、６４点では２１個、３２点では１０個、１６点では５個、８点では２個、４点では１個の正弦波のパラメータを推定することができる。

以下、具体的な実施例を説明する。

本発明の実施においては、直接的な解法は困難であると考えられるので、推定精度を必要なだけ向上させるための反復法を用いるのが実用的であると考えられる。

「解法１」（反復法１）
ステップ１
まず、単一周波数を想定して全ての隣接３周波数荷重積分値から、周波数を推定する。

ステップ２
次に、主要な周波数成分の周波数と位相とを推定する。

ステップ３
上記の主要な周波数成分（周波数の摂動を含む）による荷重積分値の分布を算出する。

ステップ４
観測した荷重積分値の分布との誤差によって推定値を更新する。

ステップ５
上記のステップ３とステップ４を必要回数繰り返す。更新する量が一定値未満になった場合に所定の精度が得られたものと判断し、終了する。

「解法２」（反復法２）
この方法の処理の流れのフローチャートが図１に示されている。

ステップ１
まず、単一周波数を想定して全ての隣接３周波数荷重積分値から、周波数と振幅と位相とを推定する。この処理は、図１中、Ｓ１−１に相当する。

ステップ２
次に、ピークを検出することによって、主要な周波数成分を取り出し（図１中ステップＳ１−２に相当する）、これら周波数成分の周波数と振幅と位相とを推定する（図１中のＳ１−３に相当する）。

ステップ３
上記の主要な周波数成分の推定に用いた荷重積分値から他の主要な周波数成分の影響を差し引く。この処理は、図１中、Ｓ１−４に相当する。

ステップ４
修正された荷重積分値から、新たに周波数と振幅と位相とを推定する。この処理は図１中、Ｓ１−５に相当する。

ステップ５
上記のステップ２の後半（Ｓ１−３）と、ステップ３（Ｓ１−４）と、ステップ４（Ｓ１−５）を必要回数繰り返す。差し引く量が一定値未満になった場合に所定の精度が得られたものと判断し、終了する。この繰り返し判断処理は図１中、Ｓ１−６に相当する。

実際にこのアルゴリズムをコンピュータ上でシミュレーションした結果が図２に示されている。図２は、Ｎ＝３２による単一並びに複数正弦波信号に対する周波数の推定結果を表すグラフである。このグラフにおいて、横軸は与えられた周波数（右上がり４５度の周波数成分である）であり、縦軸が推定された周波数である。

図２（ａ）は１個の周波数を推定する場合であり、図２（ｂ）は、２個の周波数を推定する場合、そして、図２（４）は、４個の周波数を推定する場合である。

全ての隣接する３荷重積分で単一周波数を仮定して独立に推定し、ｄｅｔＡの値が所定の第１しきい値以上のものを大きな黒点で表し、第１しきい値より小さな第２しきい値以上のものを小さな黒点で示している。複数の周波数成分が交差する付近ではうねり状の誤差が生じている。

ｈ）基本値系列を用いる荷重積分の高精度化
一般にナイキスト（Ｎｙｑｕｉｓｔ）条件を満たして間隔Δｔで標本化された標本値の時系列

は、ｓｉｎｃ関数補間された連続関数

と等価である。よって区間［−Ｔ／２，Ｔ／２］の荷重積分は、

となる。ただし、ωＴ＝２ｎπ（ｎ：整数）のように記述される。上式（数８２）中の積分のｔ＞０をグラフで表したものが、図３に示されている。

この図３においては、ｓｉｎｃ関数の積分関数（正弦積分）と、複素正弦波荷重のｓｉｎｃ積分関数（ω＝２π／４）のグラフが示されている。実線がＳｉ（ｔ）であり、破線がＳｉ_ｕ（ｔ）の実部であり、一点鎖線はＳｉ_ｕ（ｔ）の虚部である。

すなわち、ｔ＜０でほぼ零、ｔ＝０では１／２、ｔ＞０では、ほぼ１の値をとる。しかし、原点付近では波打があり、そのｔ＞０における１との交差、及び、ｔ＜０における０との交差は、ｔの整数点のほぼ中間点に位置する。このことから、積分範囲を標本点の中間に取ることによって精度を格段に向上させることが可能である。

なお、本特許では、解析の対象である信号は、所定のサンプリング周期でサンプリング（標本化）することを前提としている。

正弦積分の整数点上での値は、原点（ｔ＝０）からそれぞれ、
０．００００００（ｔ＝０）
０．５８９４９０（ｔ＝１）
０．４５１４１２（ｔ＝２）
０．５３３０９４（ｔ＝４）
０．４７４９７０（ｔ＝５）
０．５２０１０８（ｔ＝６）
０．４８３２０６（ｔ＝７）
となる。中間点での値は、ｔ＝０．５からそれぞれ、
０．４３６３２８（ｔ＝０．５）
０．５１１９６１（ｔ＝１．５）
０．４９５２３７（ｔ＝２．５）
０．５０２５１９（ｔ＝３．５）
０．４９８４５２（ｔ＝４．５）
０．５０１０４７（ｔ＝５．５）
である。

このように、積分範囲を標本化メッシュに対して標本点の中間位置までシフトした場合の荷重分布と誤差分布のグラフが図４に示されている。

図４（１）は、積分境界を標本化メッシュの中間に取った場合の誤差分布を表すグラフである。その横軸は、サンプリング点を表し、縦軸は誤差を表す。また、図４（２）は、積分境界を標本化メッシュの中間に取った場合の荷重分布を表すグラフである。その横軸は、サンプリング点を表し、縦軸は荷重を表す。

この図４に示すように、荷重分布のばらつきは大きく減少し、境界に隣接する標本点に小さな誤差が残るのみとなる。この荷重は、境界の外側隣接点で０．５−０．４３６３２８＝０．０６３６７２となり、内側隣接点で０．５＋０．４３６３２８＝０．９３６３２８となる。上記の補償は−１、０及びＮ−１、Ｎの標本値に該当し、適用には前後に計２画素（サンプル）の拡張が必要である。補償については次に述べる。

補償処理
さて、このような積分範囲のシフト（積分範囲を標本点の中間位置までずらす）の影響を補償するために、以下のような処理を行う。まず、積分領域の標本化メッシュに対するシフト量をεとすると、積分境界値項は

のように位相シフト項

が乗じられる。よって、標本化間隔Δｔを用いて、補償量は、具体的には、

となることが理解されよう。

ｈ）推定法の性質についての検討
以上、種々の推定法についてそれぞれ検討を行ったが、推定法は以下のような性質を有していると想像される。

１．無雑音、単一周波数の正弦波で、積分によるフーリエ係数を用いる場合は、厳密解が得られる。

２．標本値を用いた離散フーリエ変換を利用しても、標本化密度を十分に高めることによって、厳密解にいくらでも近づけることが可能である。

３．無雑音、含まれる周波数の個数が既知、積分によるフーリエ係数を観測地に用いる場合は厳密解が得られる。

第３ 正弦波周波数推定の応用
以上述べたように、本実施の形態によれば、信号中に含まれる正弦波のパラメータを推定（周波数、振幅、位相）できる。このような正弦波の推定は、音声信号処理の基本処理とも言うべきものであり、その応用範囲は極めて広い。応用できると考えられる分野、用途、目的としては以下のようなものが考えられる。

・分析区間によらない厳密な周波数推定を測定。

・調和関係にない複数正弦波の周波数の測定
・一般調和解析
・音声、楽器音、異常音の時間周波数分析
・音の立ち上がり、立ち下がりの分析
・時間応答による残響分析
・周波数精度と周波数分解能の明確な概念的な分離
・周波数変調の復調、位相変調の微分復調
・ドップラ検出、反射波による振動検出、音源の振動検出
本発明は、音声信号の解析だけでなく、その他各種信号に応用できることは言うまでもない。また、本特許においては、正弦波のパラメータとは、周波数と振幅と位相を意味する。

なお、正弦波のパラメータという場合、周波数と振幅と位相とに加えて直流バイアスを意味する場合もある。

本特許においては、この直流バイアスを直接パラメータとして扱ってはいないが、周波数が直流という特別の場合の振幅として推定することができる。

第４ 移動フーリエ変換
時間積分型音源恒等式に必要な荷重積分、すなわち有限フーリエ変換には、窓関数が不要という顕著な特徴がある。これを数値計算法やハードウェア信号処理に利用することを検討する。

上述した特徴は重なりのある密な短時間フーリエ変換の計算の場合に特に顕著に発揮される特徴である。すなわち、この計算は、各標本点で複素指数関数

を乗じ、これらを区間幅Ｎで移動平均することと備えられる。荷重が「１」で一定であるため、この移動平均は、各標本点で、区間の前端と後端の差を累積メモリに加算するだけで実行可能である。これを、ここでは、移動（有限）フーリエ変換と称することにする。

このときに生じる分析区間に応じた初期位相のずれは、後で補償すれば問題はない。

４．１ 演算量の評価
移動フーリエ変換の乗算回数と加算回数（複素数）は、各標本点でＮ＋Ｎ＝２Ｎ回である。全データ長さをＭとすると、重なりなしのＦＦＴで乗算回数は

であり、完全に密な移動演算とすると、

である。これに対して上記の移動フーリエ変換法の乗算回数は、

であるので、ＦＦＴをそのまま用いるよりも少ない値となる。もし、窓関数が必要であれば、窓付きでＦＦＴを完全に密に適用するには、乗算回数

が必要である。これに対して上記の移動フーリエ変換法の場合は、

であり、やはり少なくなる。

次に、密な窓関数付きＦＦＴと、窓関数なしの移動フーリエ変換法との比較では、乗算回数では

倍の開きがある。Ｎ＝２５６の場合は、上記値は、２５６×８＝２０４８となり、２０４８倍の開きとなる。また、窓関数なしのＦＦＴと移動フーリエ変換法とでは

倍の開きがある。もし、Ｎ＝２５６の場合は、上記値は８となり、８倍の開きとなる。

すなわち、従来法（例えばフーリエ位相法）で密に遅延推定を行う場合に比べて３桁の高速化が可能となると期待される。また、音源恒等式を密にＦＦＴを用いて行う場合と比較しても、少なくとも１桁の高速化が期待される。

４．２ 移動平均の漸化式の計算法
一般に平均化の荷重分布を

とすると、区間［−Ｔ／２、Ｔ／２］の移動平均は、

であり、

である。ところが、一方、

である。ここで、

は区間内及びその両端で連続かつ微分可能である。したがって、

では、

となる。すなわち、

のような漸化式による計算が可能である。この漸化式の右辺の中括弧は単に積分の境界における信号値の差であり、区間内の値を用いずに容易に計算できる。

一方、離散値を用いる場合は、上記の漸化式は厳密な式となり、厳密解が得られる。その理由は、区間［０，Ｎ−１］の

の移動平均

は

であり、

とすると、

であり、

ならば

となるからである。なお、数値計算誤差の累積を防ぐ手段を講じることが好ましい。具体的には、漸化式の直流応答を有限なものとする等の手段を高じることが好ましい。このような手段は従来から知られている工夫であるから、当業者であれば容易に実行することができる。

４．３ 移動有限フーリエ変換の実施の一例
さて、これまで説明してきた事項に基づき、移動有限フーリエ変換を実行する構成のブロック図が図５に示されている。

この図５に示すように、入力データ（観測信号）

に

を乗じてＮ系統、Ｎ段の遅延部に入力する。

図５におけるＮ個の乗算手段１０は、入力した各データに上述した荷重を乗じる。乗算した結果は、遅延部２０に供給される。この遅延部２０は、単位遅延手段２２をＮ段接続して構成されている。

遅延部２０の先頭と終端の差が累積部３０入力となり、各時刻の累積値が各時刻のＳＴＦＴ出力を与える。ここで遅延部２０の先頭と終端の差と、それまでの累積値との加算が加算手段３２によって行われる。遅延部２０の先頭のデータすなわち入力しようとするデータは符号がそのままで加算手段３２に供給され、遅延部２０の終端のデータは符号が反転されてから加算手段３２に供給される。また、累積部３０の出力信号がさらに加算手段３２に供給される。このような構成によって、加算手段３２は先頭のデータから終端のデータを減算して「差」を求めてこの差と、累積部３０の現在の累積値とを「累積」して、その結果を累積部３０に再び格納する。累積部３０はいわゆるアキュムレータから構成される。この累積部３０の出力信号がＳＴＦＴ出力を与える。

図５に示す形態のままで、前端付近と後端付近の荷重補償も行うことが可能である。また、矩形以外の窓関数では、非零の窓関数の微分の個数だけ遅延部からの乗算と加算の段数が増えることに注意されたい。演算量の低減は、窓関数が不要（境界以外で一様に−１）であるために実現されていることに注意すべきである。

４．４ 分析区間の初期位相の補償
さて、ＳＴＦＴによって求められるのは、

である。これに対して、移動フーリエ変換で求められるのは、

であり、ｔ＋τを新たにｔと置き換えることによって、

のように互いに変換可能であることが理解されよう。この際の補償位相は分析の離散周波数ごとに異なることに注意するべきである。積分範囲における１／２標本点シフト（上述したようなサンプリング周期の１／２の期間のシフト）のために位相補償の乗算は元々必要である。したがって、これに含ませることによって、この演算は演算量の増大には結びつかない。

第５．参考発明
これまで、近似ではなく厳密にパラメータを求めることができる正弦波パラメータ推定方法を説明したが、この方法は、以下のようなハードウェアを用いて実施することが好適である。

５．１ 装置
上で述べた正弦波パラメータ推定方法を実施する為に、任意の荷重関数を導入した有限区間の積分形式で、前記信号の支配方程式を構成する手段と、前記方程式を解く手段と、を含むことを特徴とする正弦波パラメータ推定支援装置を用いることが好ましい。

このような手段は、例えば、デジタルシグナルプロセッサ等で実現することが好適である。例えば、ＦＦＴを実行するする半導体チップや、暗号化処理を実行する暗号チップ等が知られているので、このような半導体チップと同様に所定の処理を高速に実施する半導体チップを作成することが好ましい。

なお、種々のアプリケーションに対応する為に、種々のサンプリング周期や、様々な音源の数に対応する為には、ソフトウェアとそのソフトウェアを実行する汎用的なプロセッサとで実現することも好ましいが、速度的には遅くなってしまう可能性がある。

また、音源の数や、データの入出力部分や、サンプリング数の設定部分のみをソフトウェアで構成し、方程式を解く演算部分のみを半導体チップを用いて構成することも好適である。このように、一部をハードウェアで、他の一部をソフトウェアで構成すれば、種々のアプリケーションに対応できると共に高速なデータ処理（信号処理）を実現できる可能性がある。

５．２ 具体的な構成
また、参考発明としては、これまで述べた正弦波パラメータ推定方法を支援する装置において、前記解析信号に前記荷重を乗算するＮ個の乗算手段と、前記Ｎ個の乗算手段毎にそれぞれ設けられ、前記乗算手段の各出力信号をＮ周期遅延させるＮ個の遅延手段と、前記Ｎ個の遅延手段毎にそれぞれ設けられ、遅延後の出力信号と、遅延前の信号との差を求め、その差を累積値に累積するＮ個の加算手段と、前記加算手段毎にそれぞれ設けられ、前記加算手段の出力信号を累積するＮ個の累積手段と、を含むことを特徴とする正弦波パラメータ推定方法支援装置である。ここで、前記Ｎは正の整数である。

このような構成は、図５に示した通りである。

第６ 方程式の他の解法
上記第２章（ｃ）ＤＣ成分を含む連立法や、同（ｄ）ＤＣ成分を含む３周波数を用いた連立法等において、方程式の種々の解法を説明した。しかし、連立方程式の解法は、他にもあり、以下、他の解法を説明する。

ａ） 単一荷重周波数と窓関数を用いる方法
本節でも、荷重積分観測値の種類が増加する欠点を承知の上で、積分境界値を未知数から除去する方法を導くことにする。出発点となる荷重積分観測方程式

において、

かつ、この窓関数ｐ（ｔ）を原点対称に選ぶと、

である。ゆえに、荷重積分観測方程式は

となる。ただし、

である。特に

となるように、窓関数を選ぶと、荷重積分観測方程式は、

となり、

のように求められる。ｐ（ｔ）の具体例としてＨａｎｎの窓関数を利用することが考えられる。Ｈａｎｎの窓関数を利用し、

とすると（ここで、

である）、原点対象で

となる。また、

である。これらのグラフが図６に示されている。図６のグラフは、横軸が時間を表し、縦軸が荷重を表す。図６（ａ）は、窓関数を示し、図６（ｂ）は荷重関数（Ｔ＝８）を表す。図６（ａ）では、実線が

を表し、破線が

を表し、一点鎖線が

をそれぞれ表す。図６（ｂ）は、これらに、ｎ＝４のｃｏｓ正弦波荷重関数を乗じたものである。

Claims (6)

1. 解析の対象である信号中の正弦波のパラメータを推定する方法において、
   前記信号の微分を含む支配方程式を立てる工程と、
   前記微分を含む支配方程式に任意の荷重関数を導入し、有限区間の積分形式に置き換える工程と、
   前記積分形式に変換した支配方程式を解く工程と、
   を含み、
   前記積分形式の支配方程式を解く工程は、
   前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程
   を含むことを特徴とする正弦波パラメータ推定方法。
2. 解析の対象である信号中の正弦波のパラメータを推定する方法において、
   任意の荷重関数を導入した有限区間の積分形式で、前記信号の微分を含む支配方程式を構成する工程と、
   前記積分形式の支配方程式を解く工程と、
   を含み、
   前記積分形式の支配方程式を解く工程は、
   前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程
   を含むことを特徴とする正弦波パラメータ推定方法。
3. 請求項１又は２記載の正弦波パラメータ推定方法において、
   前記荷重関数は、前記有限期間の整数分の１の周期の三角関数であることを特徴とする正弦波パラメータ推定方法。
4. 請求項１又は２記載の正弦波パラメータ推定方法において、
   前記積分形式の支配方程式を解く工程は、
   異なる３個の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、
   前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程と、
   を含むことを特徴とする正弦波パラメータ推定方法。
5. 請求項１又は２記載の正弦波パラメータ推定方法において、
   前記積分形式の支配方程式を解く工程は、
   異なる２個の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、
   前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程と、
   を含むことを特徴とする正弦波パラメータ推定方法。
6. 請求項１又は２記載の正弦波パラメータ推定方法において、
   前記積分形式の支配方程式を解く工程は、
   ３個以上の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、
   前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する前記工程と、
   を含むことを特徴とする正弦波パラメータ推定方法。

Images