(57)【要約】【目的】簡単な処理により一般調和解析を行え、該解析結果を用いて音楽信号の周波数分析を行う「一般調和解析装置および周波数分析装置」を提供することである。【構成】一般調和解析により音響信号の周波数分析を行う周波数分析装置である。周波数分析に際して、所定速度でサンプリングした音響信号のサンプリングデータを入力されてフーリエ変換処理により複数の周波数ポイントのフーリエ係数を演算し、かつ、各周波数ポイントのそれぞれにおいて、一般調和解析による残差エネルギーが最小となる正弦波の振幅値を前記各周波数ポイントのフーリエ係数を用いて演算する。ついで、各周波数ポイントにおける音響信号の推定エネルギーを前記各周波数ポイントの振幅値を用いて演算し、推定エネルギーが極大値を示す周波数を決定し、該決定した各周波数における音響信号のパワーまたは位相を該周波数における前記振幅値を用いて演算する。【選択図】図４
2010185682.tif 000002
【特許請求の範囲】
【請求項１】
音響信号の一般調和解析を行う一般調和解析装置において、
音響信号を所定速度でサンプリングしたサンプリングデータを入力されてフーリエ変換処理により複数の周波数ポイントのフーリエ係数を演算するフーリエ変換処理部、
前記周波数ポイントのそれぞれにおいて、一般調和解析による残差エネルギーが最小となる正弦波の振幅値を前記各周波数ポイントのフーリエ係数を用いて演算する振幅計算部、
各周波数ポイントにおける音響信号の推定エネルギーを前記各周波数ポイントの振幅値を用いて演算する推定エネルギー演算部、
推定エネルギーが極大値を示す周波数を決定する周波数決定部、
前記決定した周波数の振幅値リストを出力する振幅出力部、
を備えたことを特徴とする一般調和解析装置。
【請求項２】
前記フーリエ変換処理部は、音響信号の1観測区間に含まれるサンプリングデータ数により定まる周波数ポイントの数を整数倍して各周波数ポイントにおけるフーリエ係数を高速演算する高速フーリエ変換処理部である、
ことを特徴とする請求項１記載の一般調和解析装置。
【請求項３】
一般調和解析により音響信号の周波数分析を行う周波数分析装置において、
音響信号を所定速度でサンプリングしたサンプリングデータを入力されてフーリエ変換処理により複数の周波数ポイントのフーリエ係数を演算するフーリエ変換処理部、
前記周波数ポイントのそれぞれにおいて、一般調和解析による残差エネルギーが最小となる正弦波の振幅値を前記各周波数ポイントのフーリエ係数を用いて演算する振幅計算部、
各周波数ポイントにおける音響信号の推定エネルギーを前記各周波数ポイントの振幅値を用いて演算する推定エネルギー演算部、
推定エネルギーが極大値を示す周波数を決定する周波数決定部、
前記決定した各周波数における音響信号のパワーまたは位相を該周波数における前記振幅値を用いて演算して出力する演算部、
を備えたことを特徴とする周波数分析装置。
【請求項４】
前記フーリエ変換処理部は、音響信号の1観測区間に含まれるサンプリングデータ数により定まる周波数ポイントの数を整数倍して各周波数ポイントにおけるフーリエ係数を高速演算する高速フーリエ変換処理部である、
ことを特徴とする請求項３記載の周波数分析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は一般調和解析装置および周波数分析装置に関わり、特に、音響信号の一般調和解析を行う一般調和解析装置および一般調和解析により音響信号の周波数分析を行う周波数分析装置に関する。
【背景技術】
【０００２】
周波数分析方法として高速フーリエ変換(Fast Fourier Transform：FFT)が知られている。しかし、この周波数分析方法は周期的で調和的な信号を分析する理論であり、信号を観測した観測区間の外でも観測区間内の信号が周期的に繰り返すことを仮定している。したがって、FFTによる周波数分析方法では、観測区間によっては実際と異なった多数の周波数成分が抽出される問題がある。また、FFTでは、周波数分解能が観測区間の長さと反比例する(観測区間が長くなるほど分解能が向上する) ことから、過渡的な信号の分析を行う場合でも観測区間長を余り短くすることができない。さらに、観測区間の基本周波数（観測区間周期の逆数）の整数倍の周波数の振幅レベルを求めることを基本としていることから、基本周波数と一次高調波間の低周波信号の分析が困難となる問題がある。
【０００３】
一方、調和的でない信号に拡張したフーリエ解析の理論として一般調和解析（Generalized Harmonic Analysis：GHA）が知られている。このGHAは、観測区間内において原波形から残差エネルギーが最小となる最も優勢な(主要な)正弦波成分を抽出し、残差成分にも同様の処理を繰り返すという解析方法である。このため、GHAによれば、定常的でない僅かな周波数変動に対しても正確な周波数分析の抽出が可能であり、観測区間長と周波数の分解能は互いに独立して自由な設定が可能であるなどの特徴があるが、演算量が膨大である。すなわち、一般調和解析（GHA）は高速フーリエ変換（FFT）などの調和解析よりも精度の高い周波数分析が可能であるが、観測区間内の原波形から残差エネルギーが最小となる優勢な正弦波を抽出しなければならないことから、特定すべき正弦波の周波数、振幅、位相の３変数について繰返し計算により評価する必要があり、演算量が極端に増加するという問題がある。
【０００４】
このため、一般調和解析による精度の高い周波数分析によって符号化効率を大幅に向上させながら演算量を抑える符号化装置が提案されている(特許文献１参照)。この従来技術は、一般調和解析によって信号の周波数成分抽出を行う際、フレームをFFTして最もレベルの高い周波数を１個ないし複数個選択することによって、抽出の対象となる周波数帯域をあらかじめ限定するものである。すなわち、従来技術は一般調和解析する帯域を限定することにより一般調和解析による演算量の増加を抑えるもので、一般調和解析自体を簡単な処理により行うようにしたものではない。しかも、従来技術は抽出対象となる周波数帯域全域の一般調和解析を行わず、限定された帯域においてのみ処理するもので、全体の状況あるいはレベルの低い帯域の状況を把握できない問題がある。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００４−６９９０６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
以上から本発明の目的は、一般調和解析を簡単な処理により行なえるようにすることである。
本発明の別の目的は、一般調和解析（GHA）による精度の高い周波数分析によって分析性能を大幅に向上させながらも、演算量の増加を抑えることができるようにすることである。
本発明の別の目的は、周波数帯域全域の一般調和解析を実行しても演算量の増加を抑えることができるようにすることである。
【課題を解決するための手段】
【０００７】
本発明は、音響信号の一般調和解析を行う一般調和解析装置および一般調和解析により音響信号の周波数分析を行う周波数分析装置である。
・一般調和解析装置
本発明の一般調和解析装置は、音響信号を所定速度でサンプリングしたサンプリングデータを入力されてフーリエ変換処理により複数の周波数ポイントのフーリエ係数を演算するフーリエ変換処理部、前記周波数ポイントのそれぞれにおいて、一般調和解析による残差エネルギーが最小となる正弦波の振幅値を前記各周波数ポイントのフーリエ係数を用いて演算する振幅計算部、各周波数ポイントにおける音響信号の推定エネルギーを前記各周波数ポイントの振幅値を用いて演算する推定エネルギー演算部、推定エネルギーが極大値を示す周波数を決定する周波数決定部、前記決定した周波数の振幅値リストを出力する振幅出力部、を備えている。
前記フーリエ変換処理部は、音響信号の1観測区間に含まれるサンプリングデータ数により定まる周波数ポイントの数を整数倍して各周波数ポイントにおけるフーリエ係数を高速演算する高速フーリエ変換処理部である。
・周波数分析装置
本発明の周波数分析装置は、音響信号を所定速度でサンプリングしたサンプリングデータを入力されてフーリエ変換処理により複数の周波数ポイントのフーリエ係数を演算するフーリエ変換処理部、前記周波数ポイントのそれぞれにおいて、一般調和解析による残差エネルギーが最小となる正弦波の振幅値を前記各周波数ポイントのフーリエ係数を用いて演算する振幅計算部、各周波数ポイントにおける音響信号の推定エネルギーを前記各周波数ポイントの振幅値を用いて演算する推定エネルギー演算部、推定エネルギーが極大値を示す周波数を決定する周波数決定部、前記決定した各周波数における音響信号のパワーまたは位相を該周波数における前記振幅値を用いて演算して出力する演算部を備えている。
前記フーリエ変換処理部は、音響信号の1観測区間に含まれるサンプリングデータ数により定まる周波数ポイントの数を整数倍して各周波数ポイントにおけるフーリエ係数を高速演算する高速フーリエ変換処理部である。
【発明の効果】
【０００８】
本発明によれば、一般調和解析による残差エネルギーが最小となる正弦波の振幅値を各周波数ポイントのフーリエ係数を用いて演算し、各周波数ポイントにおける音響信号の推定エネルギーを前記各周波数ポイントの振幅値を用いて演算し、該推定エネルギーが極大値を示す周波数を決定し、該決定した各周波数における振幅値リストを出力するようにしたから、一般調和解析を簡単な処理により行なうことができるようになった。又、本発明によれば周波数帯域全域の一般調和解析を実行しても演算量の増加を抑えることができる。
また、本発明によれば、各周波数における音響信号のパワーまたは位相を該周波数における前記振幅値を用いて演算して出力するようにしたから、一般調和解析（GHA）による精度の高い周波数分析によって分析性能を大幅に向上させながらも、演算量の増加を抑えることができる。
又、本発明によれば、高速フーリエ変換処理部により、音響信号の1観測区間に含まれるサンプリングデータ数により定まる周波数ポイントの数を整数倍して各周波数ポイントにおけるフーリエ係数を高速演算するようにしたから、前記振幅演算に使用するフーリエ係数を高速に演算することができ、一般調和解析処理を高速に行なうことができる。また、周波数ポイントの数を整数倍しているため、高精度の周波数分析を行なうことができる。
【図面の簡単な説明】
【０００９】
【図１】本発明に係る一般調和解析装置の全体の処理フローである。
【図２】所定周波数の正弦波の一般調和解析（GHA）における推定エネルギー値を演算する演算処理フローである。
【図３】あらかじめ設定された周波数範囲の中から一般調和解析（GHA）による残差エネルギーが極小、換言すれば推定エネルギーが極大となる周波数を求める演算処理フローである。
【図４】第１実施例の一般調和解析装置の構成図である。
【図５】第１実施例の周波数分析装置の構成図である。
【図６】本発明によるクラッシク音楽信号(44.1KHz、16ビットのWavデータ)の周波数分析結果である。
【図７】DFT部を備えた一般調和解析装置の変形例である。
【発明を実施するための形態】
【００１０】
（A）実施例
・処理フロー
図１は、本発明に係る一般調和解析装置の全体の処理フローを説明するためのフローチャート、図２は所定周波数の正弦波の一般調和解析（GHA）における推定エネルギー値を演算する演算処理フロー、図３はあらかじめ設定された周波数範囲の中から一般調和解析（GHA）による残差エネルギーが極小、換言すれば推定エネルギーが極大となる周波数を求める演算処理フローである。
図1において,一般調和解析装置の演算処理部は、デジタル化された観測信号(例えば音響信号)D(t)から任意の観測区間で抜き出した原波形を入力として、初期化処理(ステップ１０１)、高速フーリエ変換処理及び推定エネルギー演算処理(ステップ１０２)、推定エネルギー値保存処理(ステップ１０３)、推定エネルギー値EF(fk)が極大になる周波数fkの抽出処理(ステップ１０４)、抽出周波数の正弦波の振幅リスト出力処理(ステップ１０５)により原波形に含まれる主要な正弦波の周波数、振幅、位相をリスト化して出力する。なお、高速フーリエ変換処理では、周波数ポイントの数Nsを整数倍(例えば8倍)して8×Ns(=Nf)個の各周波数ポイントfk(k＝0〜Nf-1)おけるフーリエ係数ak、bkを高速演算する。
観測信号D(t)のデジタル化サンプリング周期Ts、観測区間長t0は分析すべき正弦波の周波数及び過渡応答分析などの分析目的から適切に決定すれば良い。
【００１１】
図２において,一般調和解析装置の演算処理部は、観測区間で決まる基本周波数f0の1/8の周波数(=f0/8)の整数倍毎にフーリエ係数を求める拡張された高速フーリエ変換処理（FFT）と(ステップ２０１〜２０２)、このフーリエ係数を利用して一般調和解析（GHA）による残差エネルギーが最小となる周波数を算出するために必要な推定エネルギー値計算処理（ステップ２０３〜２０７）を実行する。
図３において、一般調和解析装置の演算処理部は、図２の処理フローで求められた一般調和解析（GHA）における予め定められた周波数範囲内の複数の推定エネルギー値を配列データとして入力し、一般調和解析（GHA）による残差エネルギーが極小、換言すれば、推定エネルギーが極大となる周波数を求める(ステップ３０１〜３０９)。
【００１２】
以下、本発明の図１〜図３の処理フローにおいて使用される数式について説明する。
原波形を時間tの関数としてD(t)と表現すると一般にD(t)は次式で表すことが出来る。
【数１】
2010185682.tif 000003
ただし、fn：任意の周波数、An、Bn：任意の振幅、  n=1,2,3,・・・である。
各正弦波の位相φnはφn=tan-1(Bn/An)と表すことができ、２つ振幅An、Bn により位相情報を演算できる。
【００１３】
図２に示す高速フーリエ変換処理(ステップ２０１〜２０２)において、原波形の観測区間の周期t0で決まる基本周波数f0の1/8の周波数(=f0/8)毎にフーリエ係数を求める。任意の周波数fkでのフーリエ係数はフーリエ変換の理論より次式
【数２】
2010185682.tif 000004
【数３】
2010185682.tif 000005
となる。高速フーリエ変換（FFT）の理論より複素数表現をした場合、フーリエ係数は次式
【数４】
2010185682.tif 000006
となる。ここで任意の周波数fkを基本周波数f0の1/8を用いて次式
【数５】
2010185682.tif 000007
に示すように表現する。観測区間のデータ数をNsとする(4)式のckは(6)式
【数６】
2010185682.tif 000008
により演算される。(6)式から高速フーリエ変換（FFT）により8Ns個ポイントの各周波数fkのフーリエ係数ak、bkは次式
【数６ａ】
2010185682.tif 000009
により求めることができ、求めたフーリエ係数ak、bkは配列データとして保存する(ステップ２０２)。
【００１４】
前記周波数fkの正弦波での一般調和解析（GHA）による残差エネルギーE(fk,Ak,Bk)は次式
【数７】
2010185682.tif 000010
で与えられる。上式の右辺を展開すると、
【数８】
2010185682.tif 000011
ここで、θk=2πfkt0、Ks=sin(2πfkt0)/(2πfkt0) と置くと(8)式の残差エネルギーE(fk,Ak,Bk)は次式
【数９】
2010185682.tif 000012
で表現できる。
【００１５】
残差エネルギーE(fk,Ak,Bk)を振幅Ak及びBkで偏微分したときのそれぞれの偏微分値がゼロとなる振幅Ak、Bkにおいて一般調和解析（GHA）による残差エネルギーE(fk,Ak,Bk)が最小となる。振幅Ak、Bkの偏微分値はそれぞれ次式で表わされる。
【数１０】
2010185682.tif 000013
【数１１】
2010185682.tif 000014
(10)，(11)式をそれぞれ０にした連立方程式より振幅Ak、Bkを求めれば、解が周波数fkの正弦波において、一般調和解析（GHA）による残差エネルギーE(fk,Ak,Bk)が最小となる振幅値となる。この振幅値Ak、Bkは次式
【数１２】
2010185682.tif 000015
【数１３】
2010185682.tif 000016
で与えられる(ステップ２０４、２０５)。
【００１６】
(12)式、(13)式で示されるように周波数fkの正弦波での一般調和解析（GHA）による残差エネルギーE(fk,Ak,Bk)が最小となる振幅Ak、Bkを周波数fkのフーリエ係数ak、bkを利用することにより複雑な繰り返し計算をすることなく簡単に求めることが出来る。これが本発明の最大の特徴である。
(9)式の右辺の第１項は原波形のエネルギー、第２項以降は任意の周波数fkに関して推定したエネルギー（推定エネルギー）と考えることが出来る。ここでその推定エネルギー値をEf(fk)とすると、次式
【数１４】
2010185682.tif 000017
で計算することができる(ステップ２０６)。
【００１７】
(14)式を(9)式に代入することにより、周波数fkの正弦波での一般調和解析（GHA）による残差エネルギーE(fk,Ak,Bk)は次式
【数１５】
2010185682.tif 000018
となる。
kを0，1，2、・・・Nf-1(=8Ns-1)として全周波数fk(k=0〜Nf-1)の推定エネルギー値Ef(fk)を(14)式により演算しフーリエ係数同様に配列データとして保存する(ステップ２０３〜２０７)。残差エネルギーE(fk,Ak,Bk)が極小となる周波数fkは、(15)式より右辺第2項の推定エネルギー値Ef(fk)が極大となる周波数fkである。
【００１８】
図３に示す処理フローでは、原波形からの残差エネルギーE(fk,Ak,Bk)が極小となる周波数fk、換言すれば、推定エネルギー値Ef(fk)が極大となる周波数fkを求める。すなわち、ステップ３０１〜３０９において、推定エネルギー値Ef(k)（k=0〜8Ns-1）の配列データから極大になる周波数を決定する。Ef(k)はf0/8刻みのk番目の周波数fkにおける推定エネルギー値であるが、フローではkの代わりにnを用いている。
推定エネルギー値Ef(n)が前の値Ef(n-1)より大きくあるいは等しく、かつ、次の値Ef(n＋1)より小さくなるときの周波数fnを求めれば、該周波数fnは残差エネルギーE(fn,An,Bn)が極小となる周波数となる(ステップ３０３〜３０７)。そして、極大となる周波数fnが求まり毎にNdをカウントアップし、Ndにより残差エネルギーが最小になる周波数の数を表示言する。
【００１９】
・構成
図４は第１実施例の一般調和解析装置の構成図である。
周波数ポイント整数倍FFT部１は、音響信号を所定速度でサンプリングしたサンプリングデータを入力され、1観測区間(フレーム区間)に含まれるサンプリングデータ数Nsにより定まる周波数ポイントの数を整数倍(例えば8倍)してNf(=8×Ns)個の各周波数ポイントfk（k＝0〜Nf-1）におけるフーリエ係数を高速演算し、演算したフーリエ係数ak，bkをフーリエ係数保存部２に保存する(図2のステップ２０１、２０２)。
一般調和解析振幅計算部３は、(12),(13)式により、前記周波数ポイントfk（k＝0〜Nf-1）のそれぞれにおいて、一般調和解析による残差エネルギーが最小となる正弦波の振幅値Ak，Bkを前記各周波数ポイントのフーリエ係数ak，bkを用いて演算し、振幅記憶部４に保存する(ステップ２０４，２０５)。推定エネルギー計算部５は、各周波数ポイントfk（k＝0〜Nf-1）における音響信号の推定エネルギーEkを前記各周波数ポイントの振幅値Ak，Bk及びフーリエ係数ak，bkを用いて演算して推定エネルギー記憶部６に保存する(ステップ２０６)。
エネルギー極大周波数計算部７は、推定エネルギーが極大値を示す周波数が、残差エネルギーが最小となる最も主要な正弦波成分の周波数であると見なして、該周波数F1〜FNdを算出してエネルギー極大周波数記憶部８に保存する(ステップ３０１〜３０９)。振幅リスト出力部９は、周波数F1〜FNdの正弦波の振幅値AmK，BmK（K=1〜Nd）を振幅記憶部４より読み出してリストにして出力する。
【００２０】
図５は第１実施例の周波数分析装置の構成図であり、図４の一般調和解析装置と同一部分には同一符号を付している。異なる点は、振幅リスト出力部９に替えて、パワー・位相計算部１０、パワー・位相記憶/出力部１１を設けた点である。パワー・位相計算部１０は、各周波数F1〜FNdの正弦波の振幅値AmK，BmK（K=1〜Nd）を用いて次式
PK＝(AmK2＋BmK2)1/2
によりパワーPKを計算し、又、次式
φK＝tan-1（BmK/AmK）
を決定する位相φKを計算し、パワー・位相記憶/出力部１１は各周波数F1〜FNdの正弦波の振幅値AmK，BmK、パワーPK、φKを保存すると共に適宜出力する。なお、推定エネルギー値Ef(k)を出力することもできる。
図示しない表示装置は、出力された各周波数のパワーPK(あるいは推定エネルギーEf(k))、φKを表示する。
【００２１】
図６は本発明によるクラッシク音楽信号(44.1KHz、16ビットのWavデータ)の周波数分析結果であり、（A）はWavデータ波形、（B）は本発明の周波数分析結果であり、横軸が周波数、縦軸は推定エネルギー値Ef(k)である。本例では、FFT対象区間のポイント数は2000点(45.36msec)で分析している。この分析結果から、信号波形に含まれる周波数成分が正確に抽出されていることがわかる。又、本発明では、一般調和解析を使用しているが高速処理が可能なため、リアルタイムでの周波数分析が可能で、演奏を聴きながら周波数分析結果を見ることができる。
【００２２】
（B）変形例
実施例では、ステップ２０２において、周波数ポイント整数倍FFT部1が(6a)，(6b)を用いて周波数ポイントを整数倍し、各周波数ポイントにおけるフーリエ係数をFFTにより演算しているが、必ずしもフーリエ係数をFFTにより演算する必要はなく、離散フーリエ変換（ＤＦＴ）によりフーリエ係数ak、bkを求めることができる。
図７は、周波数ポイント整数倍FFT部1に替えてDFT部２１を設けた一般調和解析装置の変形例であり、他の構成は図４の実施例の構成と同じである。なお、図示しないが、周波数分析装置も同様に変形することが可能である。
以上本発明によれば、一般調和解析を簡単な処理により行なうことができる。又、本発明によれば周波数帯域全域の一般調和解析を実行しても演算量の増加を抑えることができる。また、本発明によれば、各周波数における音響信号のパワーまたは位相を該周波数における前記振幅値を用いて演算して出力するようにしたから、一般調和解析（GHA）による精度の高い周波数分析によって分析性能を大幅に向上させながらも、演算量の増加を抑えることができる。
【符号の説明】
【００２３】
１  周波数ポイント整数倍FFT部
２  フーリエ係数保存部
３  一般調和解析振幅計算部
４  振幅記憶部
５  推定エネルギー計算部
６  推定エネルギー記憶部
７  エネルギー極大周波数計算部
８  エネルギー極大周波数記憶部
９  振幅リスト出力部
【図１】
2010185682.tif 000019
【図２】
2010185682.tif 000020
【図３】
2010185682.tif 000021
【図４】
2010185682.tif 000022
【図５】
2010185682.tif 000023
【図６】
2010185682.tif 000024
【図７】
2010185682.tif 000025