JP6090000B2

JP6090000B2 - 周波数解析装置

Info

Publication number: JP6090000B2
Application number: JP2013129682A
Authority: JP
Inventors: 勇登河地; 隼人北
Original assignee: Denso Wave Inc
Current assignee: Denso Wave Inc
Priority date: 2013-06-20
Filing date: 2013-06-20
Publication date: 2017-03-08
Anticipated expiration: 2033-06-20
Also published as: JP2015004573A

Description

本発明は、周波数解析装置に関し、特に、一般調和解析を用いる周波数解析装置に関する。

観測信号の周波数解析方法として、高速フーリエ変換（FFT ：Fast Fourier Transform）が知られている。FFTは観測信号が複数の調和的な正弦波に分解できると仮定して周波数解析を行うフーリエ変換を、コンピュータ上で高速に実行する方法である。しかし、周波数分解能が観測信号の観測時間長の逆数に依存することから、過渡的な観測時間の短い信号の解析は困難であるという問題点がある。

一方、非調和的な信号に対する高分解能の周波数解析方法として、N.Wienerによって提唱された一般調和解析（GHA：Generalized Harmonic Analysis）が知られている。一般調和解析は、観測信号と正弦波モデルを比較し、残差が最小になる正弦波モデルを求め、求めた正弦波モデルを観測信号から繰り返し減算していくことで周波数解析を行う方法である。観測信号と正弦波モデルを単純に比較する方法なので、過渡的で周期的でない観測時間の短い信号に対しても有効である。

正弦波モデルのパラメータには振幅・周波数・位相の３つがあり、これらのパラメータを数値探索によって求めるには膨大な計算量が必要となるが、近年のコンピュータの処理速度の向上などにより、GHAが注目されてきている。

しかし、GHAは観測信号から振幅の強い正弦波を一本ずつ繰り返し減算していくという解析方法であるため、近接周波数成分の抽出が難しいことが知られている（非特許文献１）。

Teruo Muraoka, Takahiro Miura, Tohru Ifukube, "Ambience sound recording utilizing dual MS(Mid-Side) microphone systems based upon Frequency dependent Spatial Cross Correlation(FSCC)", 122nd Convention of Audio Engineering Society 2007, 2007.5, Vienna, Austria.

観測信号に含まれている正弦波信号の周波数が互いに十分に離れていれば、一つの正弦波を観測信号から減算しても、観測信号に含まれている他の正弦波への影響は少ない。

しかし、観測信号に含まれている複数の正弦波の周波数が互いに近接していると、それら近接した正弦波は干渉し合うため、両者の中間付近の周波数の正弦波を観測信号から減算してしまう。この場合、減算後の観測信号である残差信号は、観測信号に含まれる正弦波が正しく抽出されていない意味のない信号になってしまう。そのため、通常のGHAでは、観測信号に含まれている複数の正弦波の周波数が互いに近接していると、それら複数の正弦波を正しく抽出することができないのである。

このように、通常のGHAは、過渡的な観測時間の短い信号であっても精度よく成分波形を抽出することができる特徴を有するものの、観測信号に含まれている複数の正弦波の周波数が互いに近接していると、それら複数の正弦波を正しく抽出することができないという問題がある。

また、すでに、説明したように、通常のGHAは、計算量が膨大になるという問題もある。

本発明は、この事情に基づいて成されたものであり、その目的とするところは、過渡的な観測時間の短い信号に含まれる近接周波数成分を精度よく決定でき、且つ、演算量を少なくすることもできる周波数解析装置を提供することにある。

その目的を達成するための本発明は、観測信号を取得する信号取得手段（１０）と、
周波数が互いに異なる複数の単一周波数波が合成された参照波であって、前記単一周波数波は、互いに同じ周波数の正弦関数と余弦関数の和であり、それら正弦関数および余弦関数にそれぞれ正弦振幅および余弦振幅が乗じられた式で記述されており、
前記観測信号と参照波との差に対応して変化する残差減少量が最大となるという条件に基づいて設定した連立一次方程式であって、前記正弦振幅および前記余弦振幅の合計数以上の式からなり、観測信号、観測時間長、正弦振幅、余弦振幅、および前記参照波を構成する単一周波数波の各周波数の関係を表す連立一次方程式に、
予め設定した観測時間長、信号取得手段が取得した観測信号、周波数探索範囲および周波数分解能に基づいて定まる一組の探索周波数の組み合わせを代入して、前記連立一次方程式を解くことで、前記一組の探索周波数の組み合わせに対して、前記正弦振幅および前記余弦振幅を演算し、
前記一組の探索周波数の組み合わせを前記周波数探索範囲および周波数分解能に基づいて変化させつつ前記連立一次方程式を解くことで、前記一組の探索周波数の組み合わせごとに前記正弦振幅および前記余弦振幅を決定する振幅決定手段（２１）と、
前記振幅決定手段において、一組の探索周波数の組み合わせごとに決定した前記正弦振幅、前記余弦振幅、およびそれら正弦振幅、余弦振幅を決定するために用いた一組の探索周波数の組み合わせを、予め設定した残差減少量算出式に代入することで、一組の探索周波数の組み合わせごとに、前記残差減少量を算出する残差減少量算出手段（２２）と、
前記残差減少量算出手段が一組の探索周波数の組み合わせごとに算出した残差減少量が最大となる一組の探索周波数の組み合わせと、残差減少量が最大となる一組の探索周波数の組み合わせから前記振幅決定手段で決定した正弦振幅、余弦振幅と、により定まる複数の単一周波数波を、前記観測信号に含まれる複数の周波数成分に決定する周波数成分決定手段（２３）と、を備えることを特徴とする周波数解析装置である。

本発明では、観測信号と参照波との差に対応して変化する残差減少量が最大となるという条件に基づいて連立一次方程式を設定している。そして、この連立一次方程式を解くことで、観測信号に含まれる単一周波数波の正弦振幅および余弦振幅を決定している。残差減少量が最大となるという条件は、通常のGHAの手法、すなわち、観測信号から残差が最小となる正弦波を減算することと同じである。よって、従来のGHAと同様、過渡的な観測時間の短い信号でも周波数成分の解析が可能である。

さらに、本発明では、観測信号と比較する参照波を、一つの単一周波数波とするのではなく、周波数が互いに異なる複数の単一周波数波が合成された波としている。

そして、上記連立一次方程式を解くことで、それら複数の単一周波数波の正弦振幅および余弦振幅をまとめて求めている。従来のGHAのように、観測信号から減算する正弦波が単一周波数の波であると、周波数が近接しており、互いに干渉している複数の波からなる観測信号に対して、正しい周波数の正弦波を減算することができない。しかし、本発明では、観測信号と比較する参照波も、複数の単一周波数波が合成された波としている。そのため、参照波も互いに干渉している複数の波を表すことができる。よって、観測信号に含まれている複数の正弦波の周波数が互いに近接していても、それら複数の正弦波を精度よく決定することができる。

しかも、本発明では、連立一次方程式を解くことで、正弦振幅および余弦振幅を決定することができ、数値探索を行うパラメータを周波数のみとすることができる。これにより、数値探索を行うパラメータが少なくなるので、演算量を少なくすることもできる。

また、請求項２記載の発明では、前記振幅決定手段において用いる連立一次方程式は、探索周波数の三角関数を前記観測時間長で定積分する定積分項を含んでおり、
前記探索周波数毎に予め前記定積分項を計算した計算値からなるテーブルを記憶したテーブル記憶部（４０）を備え、
前記振幅決定手段は、前記連立一次方程式に代入する前記一組の探索周波数の組み合わせに基づいて定まる定積分項の計算値を前記テーブルから抽出し、抽出した計算値を前記連立一次方程式に代入して前記連立一次方程式を解くことを特徴とする。

請求項２に係る発明では、連立一次方程式の定積分項が予め計算してあるので、演算量をさらに少なくすることができる。

本発明の実施形態となる周波数解析装置１の構成図実験例における１つ目の正弦波の周波数推定結果を示す図実験例における２つ目の正弦波の周波数推定結果を示す図実験例における１つ目、２つ目の正弦波の振幅推定結果を示す図実験例における１つ目、２つ目の正弦波の位相推定結果を示す図

（２波モデルによる理論説明）
GHAは、解析対象とする観測信号の解析区間内から残差エネルギーが最小となる正弦波を抽出し、その正弦波を観測信号から引き、残差信号に同様の処理を繰り返す解析手法である。

通常のGHAでは、原波形あるいは残差波形から１波ずつ正弦波を抽出している。これに対して、本発明では複数波の合成波である参照波を観測信号から減算する。

説明を簡単にするため、観測信号x(t)が２つの正弦波x₀(t)、x₁(t)の足し合わせである場合を考える。よって、観測信号x(t)は下記式１で表せる。

また、観測信号ｘ(t)と比較する正弦波モデルを参照波w(t)とする。参照波w(t)は、近接する周波数成分の数以上の単一周波数波を持つことが必要であり、ここでは、参照波w(t)を、それぞれは単一周波数の波である２つの正弦波w₀(t)とw₁(t)を足し合わせた波形とする。A₀,A₁を正弦振幅、B₀,B₁を余弦振幅とし、f₀,f₁を探索周波数とし、位相を正弦関数と余弦関数の振幅の組み合わせで表現すると、w₀(t)とw₁(t)は、下記式２、式３で表せ、また、参照波w(t)は式４で表せる。

前述したように、GHAは、観測信号から残差エネルギーが最小となる正弦波モデル（ここでは参照波）を抽出する手法である。よって、観測時間長をLとする下記残差エネルギーEの算出式（式５）により算出される残差エネルギーEが最小となるように参照波w(t)のパラメータを求めればよい。

式５の右辺において、観測信号x(t)だけで決まる項であるx²を除外し、式６で残差減少量eを定義する。残差エネルギーEが最小となる参照波w(t)のパラメータを求めるには、残差減少量eが最大となる参照波w(t)のパラメータを求めればよい。

式６の右辺第１項および第２項をそれぞれ展開し、文字による置き換えを行う。右辺第１項は、式２〜式４を用いて展開した後に下記式８に示す文字による置き換えを行うと、式７の最終行に示すように表現できる。

式６の右辺第２項を、式４を用いて展開すると、式９の第２行になる。式９におけるw₀ ²、w₁ ²、w₀w₁は式２、３より式１０で表すことができる。

式１０を式１１に示す文字により置き換えると、式９における右辺のように表現できる。

式７の最終行および式９の右辺より、残差減少量eは式１２で表すことができる。つまり、この式１２は残差減少量算出式である。

式１２は、各振幅A₀、A₁、B₀、B₁についてみると、上に凸の二次曲線である。それらの二次曲線がいずれも極大値となれば残差減少量eは最大となる。そこで、式１２を振幅A₀、A₁、B₀、B₁による偏微分する。各振幅により偏微分した式は式１３で表すことができる。

式１３がいずれもゼロとなれば、すなわち、式１４を満たせば、ある探索周波数f₀、f₁における残差減少量eは最大となる。

式１４と式１３から、式１５に示す、未知数４つ（A₀、A₁、B₀、B₁）の４元連立一次方程式が得られる。

式１５を解くことにより、ある探索周波数f₀、f₁において残差減少量eが最大となる参照波w(t)を決定することができる。

よって、探索周波数f₀、f₁の数値を変化させながら、各探索周波数において式１５を解き、残差減少量eが最大となる探索周波数f₀、f₁の組み合わせを探すことで、観測信号x(t)に含まれる２つの正弦波（x₀(t)、x₁(t)）を求めることができる。

この手法では、参照波w(t)を２波の正弦波（（w₀(t)、w₁(t)）で定義しているため、観測信号x(t)に含まれる２つの正弦波の周波数が近接していても、観測信号x(t)に含まれる２つの正弦波を精度よく求めることができる。

（計算量を少なくする処置）
式１５に示す連立一次方程式の右辺は、式８に示されるように、観測信号x(t)と三角関数の積分であるため、観測信号x(t)が得られるたびに計算する必要がある。これに対して、式１５に示す連立一次方程式の左辺に含まれる積分は、式１１に示されるように、すべて探索周波数f₀、f₁の関数であるため、探索周波数f₀、f₁、観測時間長Lを決めておけば、事前に計算結果のテーブルを作成しておくことができる。これにより、観測信号x(t)が得られた後の計算量を少なくすることができる。

（Ｍ波モデルへの拡張）
次に、上述した２波の参照波形をＭ波に拡張して、本発明手法を説明する。ここでのＭは、３以上の整数を取り得る。

観測信号x(t)が、式１６に示すように、N個の正弦波の足し合わせである場合を考える。

このとき、参照波形を構成する正弦波の数Mは、N≦Mであるとみなせる必要がある。よって、参照波w(t)は式１７で表すことができる。なお、みなせるとは、観測信号x(t)に含まれている正弦波のうち、無視できる小さい正弦波を除外して考えることを意味する。N≦Mである理由は、本発明手法は、観測信号において分離して抽出したい正弦波がN個以上あることに対応させて、参照波もN個以上の正弦波で構成することを特徴の一つとする手法だからである。なお、Mの具体的値は、分離して抽出したい正弦波の数が分かっている場合には、分離して抽出したい正弦波の数とする。分離して抽出したい正弦波の数が不明である場合、所望の結果が得られるまで徐々にMを大きくしてもよいし、また、十分に大きいMに設定することとしてもよい。

２波モデルで説明したと同様、式１７を展開し、文字による置き換えを行う。置き換えを行う文字を、下記式１８、式１９に示す。

文字による置き換えを行うと、残差減少量eを求める式、すなわち、残差減少量算出式は式２０で表すことができる。

この式２０を、２波モデルと同様、振幅A、Bにより偏微分し、偏微分した各式が０となるようにして連立一次方程式を立てると式２１になる。（k：0〜M-1）

式２１の連立一次方程式を解くことで、あるM個の探索周波数の組み合わせにおいて残差減少量eが最大となる参照波w(t)を決定することができる。M個の探索周波数の組み合わせを変化させつつ、各探索周波数の組み合わせにおいて残差減少量eを算出する。

そして、各探索周波数の組み合わせにおいて算出した残差減少量eの中で、残差減少量eが最大となる探索周波数の組み合わせを決定する。そして、その決定した探索周波数の組み合わせと、その探索周波数の組み合わせを使って求めた振幅A_k,B_kとから、参照波W(t)を構成する複数の正弦波W_i(t)を求めることができ、これら複数の正弦波W_i(t)は、観測信号x(t)に含まれている複数の正弦波を意味する。

以上のようにして、観測信号x(t)に含まれる正弦波の数が３以上であっても、本発明の手法では、参照波をM（≧N）波の正弦波の合成波としているので、観測信号x(t)に含まれるN個の正弦波の周波数が近接していても、それらの正弦波を精度よく求めることができる。

また、本発明の手法は、複数波モデルに拡張してはいるが、GHAの手法であるため、過渡的な観測時間の短い信号に対しても、観測信号に含まれる正弦波を精度よく抽出することができる。

（実施形態）
図１に本発明の実施形態となる周波数解析装置１の構成図を示す。この図１に示すように、周波数解析装置１は、信号取得部１０、信号処理部２０、信号出力部３０、記憶部４０を備える。

信号取得部１０は、無線受信機やマイク等、電波や音波などの観測信号x(t)を検出する信号検出装置（図示せず）と接続されており、信号検出装置から観測信号x(t)を取得する。そして、取得した観測信号x(t)を信号処理部２０へ供給する。なお、上記信号検出装置がこの信号取得部１０としての機能を備え、自身が検出した観測信号x(t)を信号処理部２０へ供給してもよい。この信号取得部１０は請求項の信号取得手段に相当する。

記憶部４０は、式２１の連立一次方程式において、左辺の定積分項、すなわち、式１８に示す項を、周波数探索範囲および周波数分解能に基づいて定まる探索周波数毎に計算した計算値からなるテーブルを記憶している。なお、周波数探索範囲、周波数分解能、観測時間長Lは予め定められている一定値とする。

信号処理部２０は、ＣＰＵ、ＲＯＭ、ＲＡＭ等を備えたコンピュータであり、ＣＰＵが、ＲＡＭの一時記憶機能を利用しつつＲＯＭに記憶されているプログラムを実行することで、信号処理部２０は、振幅決定手段２１、残差減少量算出手段２２、周波数成分決定手段２３として機能する。

振幅決定手段２１は、式１７で示した参照波の振幅A,Bを決定する。式１７から分かるように、参照波の式は、式２、式３等で示す単一周波数波が合成された合成波である。また、単一周波数は、正弦振幅Aと周波数fの正弦関数の積と、余弦振幅Bと周波数ｆの余弦関数の和で記述されている。

この参照波の振幅A,Bは、式２１に示した連立一次方程式を解くことで求める。式２１に示した連立一次方程式は、残差減少量eが最大となるという条件で設定しており、残差減少量eが最大となれば、観測信号x(t)から参照波w(t)を減算した残差エネルギーEは最小となる。

式２１に示した連立一次方程式に、一組の探索周波数の組み合わせを代入することで連立一次方程式を解く。式２１に示す連立一次方程式は、正弦振幅Aおよび余弦振幅Bの総数と同一数の一次方程式からなる。そのため、未知数であるそれら正弦振幅Aおよび余弦振幅Bは、式２１の連立一次方程式から解くことができる。

式２１の連立一次方程式を解くことにより、参照波を構成する単一周波数波の正弦振幅Aおよび余弦振幅Bを、一組の探索周波数の組み合わせ毎に決定する。そして、周波数探索範囲および周波数分解能により定まる全ての探索周波数の組み合わせに対して、式２１の連立一次方程式を解いて、全ての探索周波数の組み合わせに対して参照波を構成する単一周波数波の正弦振幅Aおよび余弦振幅Bを決定する。

残差減少量算出手段２２は、振幅決定手段２１が一組の探索周波数の組み合わせごとに決定した正弦振幅A、余弦振幅B、およびそれら正弦振幅A、余弦振幅Bを決定するために用いた一組の探索周波数の組み合わせを、式２０の残差減少量算出式に代入することで、一組の探索周波数の組み合わせごとに、残差減少量eを算出する。

周波数成分決定手段２３は、残差減少量算出手段２２が算出した残差減少量eの中での最大値を決定する。そして、その最大値を算出した一組の探索周波数の組み合わせ（たとえば、f₀,f₁,f₂）と、その一組の探索周波数の組み合わせから決定した正弦振幅A（たとえば、A₀,A₁,A₂）、余弦振幅B（たとえば、B₀,B₁,B₂）を、式１７の｛｝内に示す単一周波数波の式に代入する。これにより、参照波w(t)を構成するそれぞれの単一周波数波w_i(t)を決定するとともに、その単一周波数波w_i(t)を、観測信号x(t)に含まれる周波数成分とする。

信号出力部３０は、信号処理部２０で決定した周波数解析結果を出力する部分であり、たとえば、周波数解析結果を用いた演算処理を行う演算処理装置に、周波数解析結果を出力する。また、周波数解析結果を表示装置に出力してもよい。

（実験例）
次に、実験例を説明する。次に説明する実験例は、観測信号が２つの正弦波の足し合わせであるとして、本発明を適用したシミュレーション実験である。

この実験では、観測信号に含まれる１つ目の正弦波は、振幅1.0、周波数1000Hz、位相0degとしている。２つ目の正弦波は振幅0.5、位相0deg、周波数は0〜2000Hzの範囲を50Hz刻みで変化させている。本発明手法は参照波を２波モデル、すなわち、２つの単一周波数波の合成波とした。そして、本発明手法と通常のGHAによりそれぞれ観測信号を解析した。

図２に、１つ目の正弦波の周波数を、通常のGHA（抽出１回目）、および、本発明手法によりそれぞれ解析した結果を示す。また、図３に、２つ目の正弦波の周波数を、通常のGHA（抽出２回目）、および、本発明手法によりそれぞれ解析した結果の誤差を示す。

図２に示すように、通常のGHAでは、２つ目の正弦波の周波数が1000Hz付近になると、1つ目の正弦波の周波数推定結果が乱れてしまう。これに対して、本発明では、２つ目の正弦波の周波数がどのような値であっても、ほとんど誤差なく1つ目の正弦波の周波数を推定できている。

また、図３に示すように、通常のGHAは、２つ目の正弦波についても、２つ目の正弦波の周波数が1000Hz付近になると、周波数推定結果が乱れてしまう。これに対して、本発明手法では、２つ目の正弦波の周波数がどのような値であっても、２つ目の正弦波の周波数も、ほとんど誤差なく推定できている。なお、本発明手法でも、２つ目の正弦波の周波数を1000Hzにしたときは、推定周波数が非常に大きくなっている。

この理由は、２つ目の正弦波の周波数を1000Hzにすると、1つ目の正弦波の周波数と重なり、周波数1000Hz、振幅１．５の大きな１つの正弦波が合成され、それ以外の２つ目の正弦波がほとんど存在しなくなるためである。

図４には、振幅の推定結果を示している。この図４の横軸に示す２つ目の正弦波の周波数が1000Hzであるところを見ると、本発明手法では、２つ目の正弦波の周波数が1000Hzである場合、２つ目の正弦波の推定振幅がほとんど０になっている。このことからも、図３において、２つ目の正弦波の周波数を1000Hzとしたときに本発明手法が推定している２つ目の正弦波は、小さい波であることが分かる。

また、同じ図４から、本発明手法は、２つの正弦波の振幅をいずれも正しく推定できていることが分かる。さらに、図５から、本発明手法は、２つの正弦波の位相も正しく推定できていることが分かる。

１周波数解析装置、１０信号取得部（信号取得手段）２０信号処理部、２１振幅決定手段、２２残差減少量算出手段、２３周波数成分決定手段、３０信号出力部、４０記憶部（テーブル記憶部）

Claims

観測信号を取得する信号取得手段（１０）と、
周波数が互いに異なる複数の単一周波数波が合成された参照波であって、前記単一周波数波は、互いに同じ周波数の正弦関数と余弦関数の和であり、それら正弦関数および余弦関数にそれぞれ正弦振幅および余弦振幅が乗じられた式で記述されており、
前記観測信号と参照波との差に対応して変化する残差減少量が最大となるという条件に基づいて設定した連立一次方程式であって、前記正弦振幅および前記余弦振幅の合計数以上の式からなり、観測信号、観測時間長、正弦振幅、余弦振幅、および前記参照波を構成する単一周波数波の各周波数の関係を表す連立一次方程式に、
予め設定した観測時間長、信号取得手段が取得した観測信号、周波数探索範囲および周波数分解能に基づいて定まる一組の探索周波数の組み合わせを代入して、前記連立一次方程式を解くことで、前記一組の探索周波数の組み合わせに対して、前記正弦振幅および前記余弦振幅を演算し、
前記一組の探索周波数の組み合わせを前記周波数探索範囲および周波数分解能に基づいて変化させつつ前記連立一次方程式を解くことで、前記一組の探索周波数の組み合わせごとに前記正弦振幅および前記余弦振幅を決定する振幅決定手段（２１）と、
前記振幅決定手段において、一組の探索周波数の組み合わせごとに決定した前記正弦振幅、前記余弦振幅、およびそれら正弦振幅、余弦振幅を決定するために用いた一組の探索周波数の組み合わせを、予め設定した残差減少量算出式に代入することで、一組の探索周波数の組み合わせごとに、前記残差減少量を算出する残差減少量算出手段（２２）と、
前記残差減少量算出手段が一組の探索周波数の組み合わせごとに算出した残差減少量が最大となる一組の探索周波数の組み合わせと、残差減少量が最大となる一組の探索周波数の組み合わせから前記振幅決定手段で決定した正弦振幅、余弦振幅と、により定まる複数の単一周波数波を、前記観測信号に含まれる複数の周波数成分に決定する周波数成分決定手段（２３）と、
を備えることを特徴とする周波数解析装置。
請求項１において、
前記振幅決定手段において用いる連立一次方程式は、探索周波数の三角関数を前記観測時間長で定積分する定積分項を含んでおり、
前記探索周波数毎に予め前記定積分項を計算した計算値からなるテーブルを記憶したテーブル記憶部（４０）を備え、
前記振幅決定手段は、前記連立一次方程式に代入する前記一組の探索周波数の組み合わせに基づいて定まる定積分項の計算値を前記テーブルから抽出し、抽出した計算値を前記連立一次方程式に代入して前記連立一次方程式を解くことを特徴とする周波数解析装置。