(57)【要約】【課題】周波数や振幅に変動があっても、正弦成分の周波数、その振幅及び位相を正確に検出する。【解決手段】入力信号に窓関数を乗算してフーリエ変換する。そして、フーリエ変換された各周波数ビンに対して、近接した２つの周波数ビンの振幅比と位相差を算出する。さらに、算出した振幅比と位相差から、周波数ビンのビン周波数と周波数変動パラメータと振幅変動パラメータを検出し、検出した各周波数ビンのビン周波数と周波数変動パラメータ及び振幅変動パラメータから正弦成分を統合して、入力信号に含まれる正弦成分の周波数、周波数変化率、振幅、振幅変化率、及び位相を得る。【選択図】図１
2018025628.tif 000002
【特許請求の範囲】
【請求項１】
入力信号に窓関数を乗算してフーリエ変換を行うフーリエ変換部と、
前記フーリエ変換部がフーリエ変換して出力する各周波数ビンに対して、近接した２つの周波数ビンの振幅比と位相差を算出する振幅比・位相差算出部と、
前記振幅比・位相差算出部が算出した振幅比と位相差から、周波数ビンのビン周波数と周波数変動パラメータと振幅変動パラメータを検出する検出部と、
前記検出部が検出した各周波数ビンのビン周波数、周波数変動パラメータ及び振幅変動パラメータから正弦成分を統合して、前記入力信号に含まれる正弦成分の周波数、周波数変化率、振幅、振幅変化率、及び位相を得る正弦成分統合部とを備える
信号中に含まれる周波数や振幅が変動する正弦成分の抽出装置。
【請求項２】
前記周波数変動パラメータは、周波数変動率と前記窓関数の幅とで定義したパラメータであり、
前記振幅変動パラメータは、振幅変動率と前記窓関数の幅とで定義したパラメータである
請求項１に記載の信号中に含まれる周波数や振幅が変動する正弦成分の抽出装置。
【請求項３】
前記窓関数は、窓の中点に対して対称な実数値関数で窓の両端で値が小さくなる特性を持つ関数であり、窓長を有限に設定することで、前記フーリエ変換部が短時間フーリエ変換を行う
請求項１または２に記載の信号中に含まれる周波数や振幅が変動する正弦成分の抽出装置。
【請求項４】
入力信号に窓関数を乗算してフーリエ変換するフーリエ変換処理と、
前記フーリエ変換処理によりフーリエ変換して出力する各周波数ビンに対して、近接した２つの周波数ビンの振幅比と位相差を算出する振幅比・位相差算出処理と、
前記振幅比・位相差算出処理により算出した振幅比と位相差から、周波数ビンのビン周波数と周波数変動パラメータと振幅変動パラメータを検出する検出処理と、
前記検出処理により検出した各周波数ビンのビン周波数と周波数変動パラメータ及び振幅変動パラメータから正弦成分を統合して、前記入力信号に含まれる正弦成分の周波数、周波数変化率、振幅、振幅変化率、及び位相を得る正弦成分統合処理とを含む
信号中に含まれる周波数や振幅が変動する正弦成分の抽出方法。
【請求項５】
入力信号に窓関数を乗算してフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順によりフーリエ変換して出力する各周波数ビンに対して、近接した２つの周波数ビンの振幅比と位相差を算出する振幅比・位相差算出手順と、
前記振幅比・位相差算出手順により算出した振幅比と位相差から、周波数ビンのビン周波数と周波数変動パラメータと振幅変動パラメータを検出する検出手順と、
前記検出手順により検出した各周波数ビンのビン周波数、周波数変動パラメータ及び振幅変動パラメータから正弦成分を統合して、前記入力信号に含まれる正弦成分の周波数、周波数変化率、振幅、振幅変化率、及び位相を得る正弦成分統合手順とを、
コンピュータに実装して実行させる
プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、入力信号に含まれる周波数や振幅が変動する正弦成分の周波数、その振幅及び位相を抽出する正弦成分の抽出装置及び正弦成分の抽出方法、ならびにその正弦成分の抽出方法を適用したプログラムに関する。
【背景技術】
【０００２】
音声信号などの入力信号に含まれる正弦成分の周波数、振幅、及び位相を検出する技術として、本出願の発明者は、先に特許文献１に記載の正弦成分抽出装置及び正弦成分抽出方法を提案した。
この特許文献１の技術は、入力信号に窓関数を乗算してフーリエ変換し、フーリエ変換して得た特定の周波数ビンの出力と、その周波数ビンの近傍の周波数ビンの出力との振幅比を算出する。そして、その算出した振幅比の値から、入力信号に含まれる正弦成分の周波数を検出し、検出した正弦成分の周波数から、正弦成分の振幅及び位相を検出するものである。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１５−１６５３００号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
先に提案した特許文献１に記載の正弦信号の検出処理は、高い精度で正弦信号の周波数、振幅及び位相を検出できるが、正弦信号の周波数や振幅が分析区間内で一定と見なせることが前提である。ところが、実際の音声信号などの入力信号では、分析区間内で周波数や振幅が大きく変動することがある。このため、周波数や振幅が分析区間内で変動しても、正確に正弦成分を検出できる方法の開発が望まれていた。
【０００５】
本発明は、分析区間内で周波数や振幅に変動があっても、正弦成分の周波数と、その振幅及び位相を正確に検出することができる正弦成分抽出装置、正弦成分抽出方法及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００６】
本発明の信号中に含まれる正弦成分抽出装置は、フーリエ変換部と、振幅比・位相差算出部と、パラメータの検出部と、正弦成分統合部とを備える。
フーリエ変換部は、入力信号に窓関数を乗算してフーリエ変換する。
振幅比・位相差算出部は、フーリエ変換部がフーリエ変換して出力する各周波数ビンに対して、近接した２つの周波数ビンの振幅比と位相差を算出する。
検出部は、振幅比・位相差算出部が算出した振幅比と位相差から、周波数ビンのビン周波数と周波数変動パラメータと振幅変動パラメータを検出する。
正弦成分統合部は、検出部が検出した各周波数ビンのビン周波数、周波数変動パラメータ及び振幅変動パラメータから正弦成分を統合して、入力信号に含まれる正弦成分の周波数、周波数変化率、振幅、振幅変化率、及び位相を得る。
【０００７】
本発明の正弦成分抽出方法は、入力信号に窓関数を乗算してフーリエ変換するフーリエ変換処理と、フーリエ変換処理によりフーリエ変換して出力する各周波数ビンに対して、近接した２つの周波数ビンの振幅比と位相差を算出する振幅比・位相差算出処理と、振幅比・位相差算出処理により算出した振幅比と位相差から、周波数ビンのビン周波数と周波数変動パラメータと振幅変動パラメータを検出する検出処理と、検出処理により検出した各周波数ビンのビン周波数、周波数変動パラメータ及び振幅変動パラメータから正弦成分を統合して、入力信号に含まれる正弦成分の周波数、周波数変化率、振幅、振幅変化率、及び位相を得る正弦成分統合処理とを含む。
【０００８】
本発明のプログラムは、入力信号に窓関数を乗算してフーリエ変換するフーリエ変換手順と、フーリエ変換手順によりフーリエ変換して出力する各周波数ビンに対して、近接した２つの周波数ビンの振幅比と位相差を算出する振幅比・位相差算出手順と、振幅比・位相差算出手順により算出した振幅比と位相差から、周波数ビンのビン周波数と周波数変動パラメータと振幅変動パラメータを検出する検出手順と、検出手順により検出した各周波数ビンのビン周波数、周波数変動パラメータ及び振幅変動パラメータから正弦成分を統合して、入力信号に含まれる正弦成分の周波数、周波数変化率、振幅、振幅変化率、及び位相を得る正弦成分統合手順とを、コンピュータに実装して実行させるものである。
【発明の効果】
【０００９】
本発明によれば、入力信号に含まれる正弦成分の周波数、周波数変化率、振幅変化率に依存して、フーリエ変換して得た各周波数ビンの出力の振幅と位相が変化することを利用して、その正弦成分の周波数、周波数変化率、振幅変化率を正確に求めることができる。また、本発明によると、フーリエ変換して得た周波数ビンの出力と、検出した周波数、周波数変化率、振幅変化率を利用して、検出した正弦成分の振幅と位相を抽出することができる。
【図面の簡単な説明】
【００１０】
【図１】本発明の実施の形態例である正弦成分抽出装置の構成を示すブロック図である。
【図２】本発明の実施の形態例である正弦成分抽出装置における正弦成分抽出処理を示すフローチャートである。
【図３】窓関数の特性例（項の数４の窓関数４種に対するＷ（ξ，０，０）の絶対値のデシベル）を示す図である。
【図４】窓関数による周波数ビンの応答特性の例（Ｗ（ξ，０，０），Ｗ（ξ，２０，０），Ｗ（ξ，０，１．５），Ｗ（ξ，２０，１．５）の絶対値のデシベル）を示す図である。
【図５】窓関数による周波数ビンの応答特性の例（Ｗ（ξ，０，０），Ｗ（ξ，２０，０），Ｗ（ξ，０，１．５），Ｗ（ξ，２０，１．５）の偏角）を示す図である。
【図６】第ｋ番目のビンの出力に対する第ｋ＋１番目のビンの出力の振幅の比の例を示す図である。
【図７】第ｋ番目のビンの出力に対する第ｋ＋１番目のビンの出力の位相の差の例を示す図である。
【図８】関数Ｗξ の概形を示す図である。
【図９】関数Ｗαの概形を示す図である。
【図１０】関数Ｗβの概形を示す図である。
【図１１】本発明の適用例（その１）を示すブロック図である。
【図１２】本発明の適用例（その２）を示すブロック図である。
【発明を実施するための形態】
【００１１】
［１．実施の形態例の装置構成］
図１に、本発明の実施の形態例（以下、「本例」という）の正弦成分抽出装置の構成例を示す。まず、図１を参照して正弦成分抽出装置の構成を説明するが、正弦成分抽出装置内の各部で行われる具体的な処理については後述する。
図１に示すように、本例の正弦成分抽出装置は、信号入力端子１１、窓関数入力端子１２、乗算器１３、離散フーリエ変換回路１４、記憶部１５、振幅比・位相差算出部１６、パラメータ検出部１７、正弦成分統合部１８、及び出力端子１９を備える。
【００１２】
入力端子１１は、信号を入力する端子であり、この入力端子１１に入力される入力信号が、乗算器１３に供給される。
乗算器１３は、入力端子１１から入力される入力信号に、窓関数入力端子１２から供給される窓関数をずらしながら乗算する。この窓関数入力端子１２から入力される窓関数は後段の離散フーリエ変換回路１４におけるフーリエ変換用として利用される窓関数である。なお、窓関数は実数値の偶関数である。
乗算器１３において窓関数が乗算された入力信号は、フーリエ変換部である離散フーリエ変換回路１４に供給される。離散フーリエ変換回路１４は、窓関数が乗算された入力信号を離散フーリエ変換する。なお、ここで用いられる離散フーリエ変換は、高速フーリエ変換である。
【００１３】
離散フーリエ変換回路１４がフーリエ変換した信号は、記憶部１５に供給され記憶される。記憶部１５には、離散フーリエ変換回路１４によりフーリエ変換された信号が、フーリエ変換によって生成された周波数ビンごとに記憶される。
【００１４】
そして、記憶部１５に記憶された信号は、振幅比・位相差算出部１６に供給され、各周波数ビンに対して近接した２つの周波数ビンの振幅比と位相差が算出される。すなわち、振幅比・位相差算出部１６は、それぞれの周波数ビンに対して、その前及び後の周波数ビンとの振幅比ｒ−ｉ，ｒ＋ｉ、及びその前及び後の周波数ビンとの位相比ｐ−ｉ，ｐ＋ｉを算出する。
【００１５】
次に、検出部１７は、各周波数ビンについて、その前及び後の周波数ビンとの振幅比及び位相比から、ビン周波数ξ、周波数変動パラメータα、及び振幅変動パラメータβを検出する。
さらに、正弦成分統合部１８は、検出部１７で検出された各周波数ビンでのビン周波数ξ、周波数変動パラメータα、及び振幅変動パラメータβを統合して、入力信号に含まれる正弦信号の周波数ｆ（ｔ）、周波数変化率δｆ（ｔ）、振幅Ａ（ｔ）、振幅変化率δＡ（ｔ）、及び位相φ（ｔ）を得る。
【００１６】
この正弦成分統合部１８で得られた正弦信号の周波数ｆ（ｔ）、周波数変化率δｆ（ｔ）、振幅Ａ（ｔ）、振幅変化率δＡ（ｔ）、及び位相φ（ｔ）は、出力端子１９から出力される。本例の場合には、正弦信号の周波数や振幅に変動があっても、その変動に追随した正弦信号の周波数ｆ（ｔ）及び振幅Ａ（ｔ）が検出される。また、入力信号に複数の正弦信号が含まれる場合でも、それぞれの正弦信号の周波数が互いに十分離れているときには、それぞれの正弦信号が検出可能である。
振幅比・位相差算出部１６、検出部１７、及び正弦成分統合部１８での具体的な処理の詳細は、数式を用いて後述する。
【００１７】
［２．実施の形態例の処理］
次に、図２のフローチャートを参照して、本例の正弦成分抽出装置で行われる処理の流れを説明する。
まず、離散フーリエ変換回路１４は、窓関数が乗算された入力信号の短時間フーリエ変換処理を行う（ステップＳ２１）。このステップＳ２１が短時間フーリエ変換処理となる。
ここで、ステップＳ２１の短時間フーリエ変換処理のために乗算された窓関数は、窓の中点に対して対称な実数値関数である。
【００１８】
ステップＳ２１で離散フーリエ変換回路１４が短時間フーリエ変換処理を行った後、振幅比・位相差算出部１６及び検出部１７は、全ての周波数ビンについてステップＳ２３及びステップＳ２４を行ったか否かを判断する（ステップＳ２２）。
ここで、全ての周波数ビンについてステップＳ２３及びステップＳ２４を行っていないと判定された場合には（ステップＳ２２のＮｏ）、振幅比・位相差算出部１６が、周波数ビンとその前及び後の周波数ビンとの振幅比ｒ−ｉ，ｒ＋ｉを算出する振幅算出処理を行うと共に位相比ｐ−ｉ，ｐ＋ｉを算出する位相差算出処理を行う（ステップＳ２３）。このステップＳ２３の処理が振幅比・位相差算出処理になる。
【００１９】
その後、検出部１７は、ビン周波数ξ、周波数変動パラメータα、振幅変動パラメータβを検出する処理を行い（ステップＳ２４）、ステップＳ２２の判断に戻る。ステップＳ２４の処理が、各種パラメータの検出処理になる。
【００２０】
また、ステップＳ２２で、ステップＳ２３及びステップＳ２４が終了したと判定された場合は（ステップＳ２２のＹｅｓ）、正弦成分統合部１８は、各周波数ビンで検出されたパラメータから、正弦成分を統合して、正弦信号の周波数ｆ（ｔ）、周波数変化率δｆ（ｔ）、振幅Ａ（ｔ）、振幅変化率δＡ（ｔ）、及び位相φ（ｔ）を得る（ステップＳ２５）。そして、これらのデータを出力端子１９から出力して、処理を終了する。このステップＳ２５の処理が正弦成分統合処理となる。
【００２１】
［３．具体的な検出状態の例］
次に、図２で示した本例の処理を行う具体的な例について説明する。
［３−１．周波数や振幅が変動する正弦信号の短時間フーリエ変換］
まず、入力端子１１に得られる時刻ｔにおける入力信号ｘ（ｔ）が、以下の［数１］式で示すように、振幅Ａ（ｔ）と位相φ（ｔ）を持つ正弦成分であるとする。また、瞬時周波数ｆ（ｔ）は、［数２］式のように位相φ（ｔ）の微分φ′（ｔ）で示すことができる。
【００２２】
【数１】
2018025628.tif 000003
【００２３】
【数２】
2018025628.tif 000004
【００２４】
但し、時刻ｔの近傍で振幅Ａ（ｔ），位相φ（ｔ），瞬時周波数ｆ（ｔ）は、次の［数３］式〜［数６］式に示すように変動しているものとする。
【００２５】
【数３】
2018025628.tif 000005
【００２６】
【数４】
2018025628.tif 000006
【００２７】
【数５】
2018025628.tif 000007
【００２８】
【数６】
2018025628.tif 000008
【００２９】
ここで、δＡ（ｔ）を時刻ｔでの振幅変動率、δｆ（ｔ）を時刻ｔでの周波数変動率と呼ぶ。この正弦信号に対して、窓関数ｗ（λ）を用いて、次のように短時間フーリエ変換を行う。Ｌは、窓関数ｗ（λ）を適用する時間幅である。この短時間フーリエ変換が、離散フーリエ変換回路１４で実行される。
【００３０】
【数７】
2018025628.tif 000009
【００３１】
ここで、周波数ｆ（ｔ）に対して、周波数ビンｋにおけるビン周波数ξを［数８］式で定義し、周波数変動δｆ（ｔ）に対して、周波数変動パラメータαを［数９］式で定義し、振幅変動率δＡ（ｔ）に対して、振幅変動パラメータβを［数１０］式で定義する。
【００３２】
【数８】
2018025628.tif 000010
【００３３】
【数９】
2018025628.tif 000011
【００３４】
【数１０】
2018025628.tif 000012
【００３５】
これらを用いて、関数Ｗ（ξ，α，β）を次式のように定義する。
【００３６】
【数１１】
2018025628.tif 000013
【００３７】
この関数Ｗ（ξ，α，β）を用いることで、［数７］式は、以下のように表すことができる。
【００３８】
【数１２】
2018025628.tif 000014
【００３９】
ここで、［数７］式を［数１２］式のように置き換えることができる点を、証明する。
τ＝Ｌλとすると、［数７］式は以下のように変形できる。
【００４０】
【数１３】
2018025628.tif 000015
【００４１】
さらに、［数１４］式のように変形でき、Ｘ（ｋ，ｔ）がｘ（ｔ）Ｗ（ξ，α，β）で示されることが証明される。
【００４２】
【数１４】
2018025628.tif 000016
【００４３】
周波数変化率の絶対値｜δｆ（ｔ）｜は、音声の基本周波数は、おおよそ８００Ｈｚ／ｓが上限であり、窓関数の幅Ｌが６４ｍｓのとき、周波数変動パラメータαの絶対値の上限は１０程度である。したがって、基本周波数のｎ倍高調波の｜α｜の上限は１０ｎ程度になる。
振幅変化率の絶対値｜δＡ（ｔ）｜は、音声では、おおよそ２０ｓ−１が上限である。窓関数の幅Ｌが６４ｍｓのとき、振幅変動パラメータβの絶対値の上限は１．３程度である。
【００４４】
［３−２．関数Ｗ（ξ，α，β）のパラメータαとβに対する依存性］
次に、短時間フーリエ変換に用いる窓関数の例として、以下の［数１５］式の窓関数について、周波数変動パラメータαと振幅変動パラメータβに対する依存性を説明する。なお、窓関数が他の形の関数で定義される場合も、関数Ｗ（ξ，α，β）は［数１１］式で定義され、周波数変動パラメータαと振幅変動パラメータβに対する依存性を持つ。
【００４５】
【数１５】
2018025628.tif 000017
【００４６】
図３は、この窓関数に対する関数Ｗの絶対値｜Ｗ（ξ，０，０）｜のデシベル表示を示す。
図３では、項の数ｎが４の窓関数のMinimun 4-term（以下、窓関数ｔ４ｄと称する）の例、両端で一次微分係数が連続の4-Term with Contimuous First Dirivative（以下、窓関数ｔ４ｃ１と称する）の例、３次微分係数が連続の4-Term with Contimuous Third Dirivative（以下、窓関数ｔ４ｃ３と称する）の例、５次微分係数が連続の4-Term with Contimuous Fifth Dirivative（以下、窓関数ｔ４ｃ5と称する）の例、のそれぞれに対する関数Ｗの絶対値｜Ｗ（ξ，０，０）｜のデシベル表示を示す。
この図３から、それぞれの窓関数のメインローブは、ξが−４から＋４の範囲にあり、その外側の信号は大きく減衰することがわかる。
【００４７】
ここで、Ｗ（ξ，α，β）を、周波数変動パラメータαと振幅変動パラメータβに関して、テーラ展開をするために、次の関数を定義する。
【００４８】
【数１６】
2018025628.tif 000018
【００４９】
【数１７】
2018025628.tif 000019
【００５０】
ここで、［数１６］式の関数ｓｍ（ξ，α）の計算方法の例を、以下に示す。
【００５１】
【数１８】
2018025628.tif 000020
【００５２】
［数１８］式は、さらに次の［数１９］式に続く。
【００５３】
【数１９】
2018025628.tif 000021
【００５４】
また、関数ｓ２ｍ（ξ，α）及び関数ｓ２ｍ＋１（ξ，α）は、以下のように展開できる。
【００５５】
【数２０】
2018025628.tif 000022
【００５６】
【数２１】
2018025628.tif 000023
【００５７】
［数１７］式の関数ｓ（ξ，α，β）を用いると、Ｗ（ξ，α，β）は以下のように表される。
【００５８】
【数２２】
2018025628.tif 000024
【００５９】
［数２１］式は、以下のように証明される。
【００６０】
【数２３】
2018025628.tif 000025
【００６１】
［３−２．窓関数による周波数ビンの応答特性］
次に、窓関数による周波数ビンの応答特性について説明する。
先に提案した特許文献１に記載してあるように、周波数や振幅が変化しない、すなわち、周波数変化率δｆ（ｔ）＝０かつ振幅変化率δＡ（ｔ）＝０のときには、［数１５］式で定義される窓関数が、実数値で偶関数であることから、Ｗ（ξ，０，０）は実数になる。この結果、［数１２］式より、“ｋ”に依らずにＸ（ｋ，ｔ）の位相はｘ（ｔ）の位相に等しく、振幅はｘ（ｔ）のＷ（ξ，０，０）倍である。このことを利用して、Ｘ（ｋ＋ｉ，ｔ）とＸ（ｋ，ｔ）の振幅の比からξが求まり、［数８］式から正弦信号の周波数を算出することができる。さらに、周波数ｆ（ｔ）からＷ（ｆ（ｔ）Ｌ−ｋ，０，０）を求め、Ｘ（ｋ，ｔ）からｘ（ｔ）の振幅Ａ（ｔ）と位相φ（ｔ）を求めることができる。
【００６２】
しかし、周波数や振幅が変化するときには、Ｗ（ξ，α，β）は一般には実数にならず、Ｘ（ｋ，ｔ）の位相はｘ（ｔ）の位相に等しくはない。ここで、窓関数Ｗ（ξ，α，β）の絶対値のデシベルや偏角のαやβに対する影響について述べる。
【００６３】
例えば、δｆ（ｔ）＝１５５４に対するα＝２０と、δＡ（ｔ）＝２３．４４に対するβ＝１．５について、窓関数の両端で１次微分係数が連続の窓関数ｔ４ｃ１に対するＷ（ξ，０，０），Ｗ（ξ，α，０），Ｗ（ξ，０，β），Ｗ（ξ，α，β）の絶対値のデシベルを図４に、これらの偏角を図５に示す。
【００６４】
周波数変化率δｆ（ｔ）≠０あるいは振幅変化率δＡ（ｔ）≠０の場合、本発明者が既に提案した方法（特許文献１）では、正弦信号の周波数ｆ（ｔ）の推定値や振幅Ａ（ｔ）と位相φ（ｔ）の推定値に誤差が生じる。
ここで、本発明では、周波数や振幅が変動する正弦信号ｘ（ｔ）の時刻ｔにおける周波数ｆ（ｔ）、周波数変化率δｆ、振幅Ａ（ｔ）、振幅変化率δＡ、及び位相φ（ｔ）を求める方法を開発した。次に、本発明の具体的な方法の詳細を説明する。
【００６５】
［３−３．周波数ビン間の出力の比と位相の差から正弦信号を求める処理］
まず、正弦信号の周波数ｆ（ｔ）をメインローブに含む周波数ビンをｋとする。このとき、窓関数に段落［００４６］で示した項の数が４の窓関数を用いると、ξ＝ｆ（ｔ）Ｌ−ｋと置くとき、ξは−４≦ξ≦４の範囲にある。
目的とする正弦信号の近傍に他の正弦信号がない場合は、［数１２］式から、近傍の周波数ビンｋ＋ｉ（ｉ＝±１，±２，・・・）での値及び周波数ビンｋでの値は、［数２４］式及び［数２５］式に示す条件が成り立つ。
【００６６】
【数２４】
2018025628.tif 000026
【００６７】
【数２５】
2018025628.tif 000027
【００６８】
近傍の周波数ビンｋ＋ｉでの値Ｘ（ｋ＋ｉ，ｔ）の周波数ビンｋでの値Ｘ（ｋ，ｔ）に対する振幅の比のデシベルをｒｉ、位相の差をｐｉとすると、振幅の比ｒｉ、位相差ｐｉは、それぞれ［数２６］式、［数２７］式に示すようになる。ここで、∠ｘは、複素数ｘの偏角を表す。
【００６９】
【数２６】
2018025628.tif 000028
【００７０】
【数２７】
2018025628.tif 000029
【００７１】
ここで、［数２４］式及び［数２５］式から、次の式が成り立つ。
【００７２】
【数２８】
2018025628.tif 000030
【００７３】
【数２９】
2018025628.tif 000031
【００７４】
これら［数２８］式及び［数２９］式の右辺は信号に依存せず、窓関数ｗ（ｔ）だけで定まる。そこで、これら［数２８］式及び［数２９］式の右辺に対して、以下のようにξ，α，βとｉの関数Ｗａ（ξ，α，β，ｉ）とＷｐ（ξ，α，β，ｉ）を定義する。
【００７５】
【数３０】
2018025628.tif 000032
【００７６】
【数３１】
2018025628.tif 000033
【００７７】
これらの関数を用いて、観測値ｒｉ，ｐｉに対して、以下の非線形連立方程式を解くことによって、ξ，α，βを求めることができる。ここで、ｉ＝±１，±２，・・・である。
【００７８】
【数３２】
2018025628.tif 000034
【００７９】
【数３３】
2018025628.tif 000035
【００８０】
これらのξ，α，βから、Ｗ（ξ，α，β）を計算することができ、［数１２］式からｘ（ｔ）を以下のように求めることができる。
【００８１】
【数３４】
2018025628.tif 000036
【００８２】
ξ，α，βからｆ（ｔ），δｆ（ｔ），δＡ（ｔ）を求め、［数１］式よりｘ（ｔ）の振幅Ａ（ｔ）、位相φ（ｔ）を求めることができる。
ここで、Ｗａ（ξ，α，β，ｉ），Ｗｐ（ξ，α，β，ｉ）のαやβに対する影響について考察する。
図６は、第ｋ番目のビンの出力に対する第ｋ＋１番目のビンの出力の振幅の比のデシベルを示す。図７は、第ｋ番目のビンの出力に対する第ｋ＋１番目のビンの出力の位相の差を示す。
この図６は、δｆ（ｔ）＝１５５４に対するα＝２０と、δＡ（ｔ）＝２３．４４に対するβ＝１．５について、窓関数の両端で１次微分係数が連続の窓関数ｔ４ｃ１に対するＷａ（ξ，０，０，１），Ｗａ（ξ，α，０，１），Ｗａ（ξ，０，β，１），Ｗａ（ξ，α，β，１）のデシベルを示す。また、図７は、窓関数ｔ４ｃ１に対するＷｐ（ξ，０，０，１），Ｗｐ（ξ，α，０，１），Ｗp（ξ，０，β，１），Ｗｐ（ξ，α，β，１）の位相差を示す。
【００８３】
［３−４．周波数ビン間の出力の比から正弦信号を求める手順の詳細処理］
入力信号ｘ（ｔ）の中に複数の正弦信号があっても、それらの周波数が互いに十分離れ、窓関数のメインローブ中に複数の正弦信号が含まれていない場合は、それぞれの正弦信号について、例えば、周波数ビン間の振幅の比や位相差から、以下の手順１−５で正弦信号を求めることができる。
・手順１：入力信号ｘ（ｔ）の短時間フーリエ変換Ｘ（ｋ，ｔ）（ｋ＝０，１，・・・）を求める。
・手順２：各周波数ビンｋに対して、以下の処理（ａ），（ｂ），（ｃ）を行う。
（ａ）［数２６］式を用いて振幅比ｒｉ，ｒ−ｉを求め、［数２６］式を用いて位相差ｐｉ，ｐ−ｉを求める。
（ｂ）振幅比ｒｉ，ｒ−ｉ、位相差ｐｉ，ｐ−ｉに対して、非線形方程式［数３２］式，［数３３］式の解の近似値ξ，α，βを求める。
（ｃ）求めた近似値ξ，α，βを初期値として、非線形連立方程式［数３２］式，［数３３］式を、例えばニュートン法などを用いて解き、ξ，α，βを求める。
・手順３：隣り合う一連の連続のビンには同一の正弦成分が検出されるので、これらの値を統合し、周波数ｆ（ｔ），周波数変動パラメータα，振幅変動パラメータβを定める。
・手順４：これらのξ，α，βから、Ｗ（ξ，α，β）を計算することができ、［数３４］式から入力信号ｘ（ｔ）を求める。
・手順５：［数１］式から入力信号ｘ（ｔ）の振幅Ａ（ｔ）と位相φ（ｔ）を求める。
【００８４】
［３−５．ξ，α，βの近似値を求める処理の例］
上述した手順２でｉとして１のみを用いた例を以下に示す。手順２の（ｂ）で、振幅比ｒ１，ｒ−１及び位相差ｐ１，ｐ−１から近似値ξ，α，βを求める方法としては、窓関数の特性を用いる方法、表探索法、ニューラルネットワークを用いる方法など、様々な方法が考えられる。ここでは、窓関数の特性を用いる方法について説明する。
まず、以下の関数を定義する。
【００８５】
【数３５】
2018025628.tif 000037
【００８６】
［数２６］式で定義されるＸ（ｋ＋１，ｔ）とＸ（ｋ，ｔ）の振幅比のデシベルｒ１と、Ｘ（ｋ−１，ｔ）とＸ（ｋ，ｔ）の振幅比のデシベルｒ−１との差ｒ１−ｒ−１を、ｒξと表す。
【００８７】
【数３６】
2018025628.tif 000038
【００８８】
ここで、［数２４］式，［数２６］式，［数２８］式より、以下の式が成り立つ。
【００８９】
【数３７】
2018025628.tif 000039
【００９０】
図８は、α＝２０，β＝１．５とした場合の、窓関数ｔ４ｃ１に対するＷξ（ξ，０，０），Ｗξ（ξ，α，０），Ｗξ（ξ，０，β），Ｗξ（ξ，α，β）の例を示す。図８に示すように、［数３６］式はαにあまり依存せず、βにごく僅かしか依存せず、ほぼξに直線的に依存している。
そこで、Ｗξ（ξ，α，β）を、ξで接する直線で以下のように近似する。
【００９１】
【数３８】
2018025628.tif 000040
【００９２】
【数３９】
2018025628.tif 000041
【００９３】
そして、［数３７］式の近似解を求める関数を、［数３９］式から以下の［数４０］式のように定義する。これにより、ξの近似度を向上させることができる。
【００９４】
【数４０】
2018025628.tif 000042
【００９５】
次にαの近似値を求めるために［数４１］式の関数を定義する。
【００９６】
【数４１】
2018025628.tif 000043
【００９７】
ここで、［数２７］式で定義されるＸ（ｋ＋１，ｔ）とＸ（ｋ，ｔ）の位相の差ｐ１と、Ｘ（ｋ−１，ｔ）とＸ（ｋ，ｔ）の位相の差ｐ−１との和ｐ１＋ｐ−１をｐαとしたとき、次の［数４２］式で示される。
【００９８】
【数４２】
2018025628.tif 000044
【００９９】
そして、［数２７］式，［数２９］式より、以下の［数４３］式が成り立つ。
【０１００】
【数４３】
2018025628.tif 000045
【０１０１】
図９は、α＝２，β＝１．５とした場合の、窓関数ｔ４ｃ１に対するＷα（０，α，０），Ｗα（ξ，α，０），Ｗα（０，α，β），Ｗα（ξ，α，β）の例を示す。図９に示すように、｜ｐα｜＜０．３２ならばＷα（ξ，α，β）はξにあまり依存せず、βにはごく僅かしか依存せず、αに直線的に依存している。そこで、Ｗα（ξ，α，β）を、αで接する直線で以下の［数４４］式，［数４５］式のように近似する。
【０１０２】
【数４４】
2018025628.tif 000046
【０１０３】
【数４５】
2018025628.tif 000047
【０１０４】
ここで、［数４３］式の近似解を求める関数を、［数４５］式から、以下のように定義する。
【０１０５】
【数４６】
2018025628.tif 000048
【０１０６】
次に、βについて、以下の関数を定義する。
【０１０７】
【数４７】
2018025628.tif 000049
【０１０８】
ここで、［数２７］式で定義されるＸ（ｋ＋１，ｔ）とＸ（ｋ，ｔ）の位相の差ｐ１と、Ｘ（ｋ−１，ｔ）とＸ（ｋ，ｔ）の位相の差ｐ−１との差ｐ１−ｐ−１をｐβとしたとき、ｐβは次の［数４８］式で示される。
【０１０９】
【数４８】
2018025628.tif 000050
【０１１０】
［数２７］式，［数２９］式より、以下の式が成り立つ。
【０１１１】
【数４９】
2018025628.tif 000051
【０１１２】
図１０は、ξ＝２，α１＝１０，α２＝２０とした場合の、窓関数ｔ４ｃ１に対するＷβ（０，０，β），Ｗβ（ξ，０，β），Ｗβ（０，α１，β），Ｗβ（ξ，α１，β），Ｗβ（０，α２，β），Ｗβ（ξ，α２，β）の例を示す。図１０では違いが分かりにくいが、Ｗβ（ξ，α，β）は、αが大きくなるにつれて、ξに大きく依存するようになる。また、Ｗβ（ξ，α，β）は、βに直線的に依存している。そこで、Ｗβ（ξ，α，β）を、βで接する直線で以下の［数５０］式，［数５１］式のように近似する。
【０１１３】
【数５０】
2018025628.tif 000052
【０１１４】
【数５１】
2018025628.tif 000053
【０１１５】
ここで、［数４９］式の近似解を求める関数を、［数５１］式から、以下のように定義する。
【０１１６】
【数５２】
2018025628.tif 000054
【０１１７】
そして、［数４０］式，［数４６］式，［数５２］式で定義される近似関数を用いて、［数３６］式，［数４２］式，［数４８］式で定義される観測値から、ξ，α，βの近似値を計算する方法の例を、次式に示す。
【０１１８】
【数５３】
2018025628.tif 000055
【０１１９】
この方法で求めた近似値の誤差を、表１から表９に示す。表１〜表９の各欄において、カッコ内は誤差の絶対値の最大値である。
【０１２０】
【表１】
2018025628.tif 000056
【０１２１】
【表２】
2018025628.tif 000057
【０１２２】
【表３】
2018025628.tif 000058
【０１２３】
【表４】
2018025628.tif 000059
【０１２４】
【表５】
2018025628.tif 000060
【０１２５】
【表６】
2018025628.tif 000061
【０１２６】
【表７】
2018025628.tif 000062
【０１２７】
【表８】
2018025628.tif 000063
【０１２８】
【表９】
2018025628.tif 000064
【０１２９】
これらの表で、｜ξ｜≦０．５の場合とは、観測している正弦成分の周波数がｋ番目のビンの中心周波数に最も近い場合である。０．５＜｜ξ｜＜１．５の場合とは、観測している正弦成分の周波数が、ｋ番目のビンの隣のｋ−１ビン又はｋ＋１ビンにある場合であり、１．５＜｜ξ｜≦２．５の場合とは、ｋ番目のビンから２つ離れたｋ−２ビン又はｋ＋２ビンにある場合である。
｜α｜≦１０、１０＜｜α｜≦３０、３０＜｜α｜≦５０は、信号の基本周波数の周波数変化がそれぞれこの範囲に入る場合を表している。信号の振幅変化については、表の中でβ０と記載した｜β｜≦０．４の範囲、表の中でβ１と記載した０．４＜｜β｜≦１．２の範囲、β２と記載した１．２＜｜β｜≦２の範囲にそれぞれ収まる場合を表している。
【０１３０】
第三近似では値が求まらない場合があるため、ξとαについては第二近似ξv，αvを、βについては第二近似が第一近似と変わらないことから、第一近似β^を近似値とするのが妥当である。（ξv，αv，β^の「v」 や「^」は、本来は［数５３］式のように、それぞれ各記号ξ，α，βの上に記載するものであるが、表記上の制約のため、右にずらして記載している。）
なお、ここでは、観測値ｒξ，ｐα，ｐβを与えるξ，α，βの組み合わせが、一つの組み合わせに定まらない場合がある。
【０１３１】
［３−５．パラメータξ，α，βを求める方法の例］
手順２（ｃ）を実行する際には、非線形連立方程式の［数３２］式，［数３３］式を解く必要があるが、このとき、上述したようにパラメータξ，α，βの近似値を初期値として、ニュートン法を用いる方法がある。
以下、ニュートン法を用いてパラメータξ，α，βを求める方法の一例を説明する。
既に述べたように、［数３７］式，［数４３］式，［数４９］式からなる連立方程式は、解が複数存在する場合がある。
そこで、［数３２］式，［数３３］式において、ｉを＋１と−１とした場合を組み合わせ、以下の式の中から３式を選び、ニュートン法を用いる。
【０１３２】
【数５４】
2018025628.tif 000065
【０１３３】
この［数５４］式から選んだ３つの式を、ｆｉ（ξ，α，β）＝ｖｉ，ｉ＝１，２，３とする。
ｆｉ（ξ＋δξ，α＋δα，β＋δβ）を、ξ，α，βの周りで１次の項まで展開する。
【０１３４】
【数５５】
2018025628.tif 000066
【０１３５】
現在の近似値ξ，α，βに対して、ｆｉ（ξ＋δξ，α＋δα，β＋δβ）がｖｉとなるδξ，δα，δβを以下の連立方程式から求める。次式において、ｉ＝１，２，３である。
【０１３６】
【数５６】
2018025628.tif 000067
【０１３７】
このδξ，δα，δβを用いて、近似値ξ，α，βを以下のように更新する。
【０１３８】
【数５７】
2018025628.tif 000068
【０１３９】
以下の式で定義する誤差が減少する間、方程式の組み合わせを変えながら更新を繰り返す。
【０１４０】
【数５８】
2018025628.tif 000069
【０１４１】
このニュートン法を用いて求めた誤差を、表１０から表１８に示す。表１０〜表１８の各欄において、カッコ内は誤差の絶対値の最大値である。また、ニュートン法の繰り返し回数の平均値と最大値を、表１９から表２１に示す。
【０１４２】
【表１０】
2018025628.tif 000070
【０１４３】
【表１１】
2018025628.tif 000071
【０１４４】
【表１２】
2018025628.tif 000072
【０１４５】
【表１３】
2018025628.tif 000073
【０１４６】
【表１４】
2018025628.tif 000074
【０１４７】
【表１５】
2018025628.tif 000075
【０１４８】
【表１６】
2018025628.tif 000076
【０１４９】
【表１７】
2018025628.tif 000077
【０１５０】
【表１８】
2018025628.tif 000078
【０１５１】
【表１９】
2018025628.tif 000079
【０１５２】
【表２０】
2018025628.tif 000080
【０１５３】
【表２１】
2018025628.tif 000081
【０１５４】
表１２、表１５、表１８から１．５＜｜ξ｜≦３の全ての場合、推定値の誤差が大きい。表１１、表１４、表１７から０．５＜｜ξ｜≦１．５かつ  ３０＜｜α｜≦  ５０のβ１、β２の場合も、推定誤差が大きい。これら以外の場合は、ニュートン法によって十分な精度のξ，α，βの推定値を得ることができた。
なお、非線形連立方程式の［数３２］式，［数３３］式を解く方法は、ここまで説明した推定方法以外を適用することも可能で、適用する方法によっては、より広い範囲で正確な推定値を得ることができると考えられる。
【０１５５】
［４．一実施の形態による効果］
以上説明したように、本発明の実施の形態例によると、周波数や振幅が変動する正弦信号を含む信号を短時間フーリエ変換し、各正弦信号の周波数に最も近い中心周波数を持つ周波数ビンの出力に対する近傍のビンの出力の振幅や比や位相の差から、窓関数を用いてその正弦信号の振幅Ａ（ｔ）、振幅変動率δＡ（ｔ）、位相φ（ｔ）、瞬時周波数ｆ（ｔ）、周波数変動率δｆ（ｔ）を求めることができる。
例えば、１２５Ｈｚから１５６．２５Ｈｚに変化（α：６．２８に相当）すると共に、振幅が１．０から１．５に変化（β：０．５から０．３３に相当）する正弦信号について、本実施の形態例の正弦成分抽出装置での処理を実行するプログラムを作成したところ、十分な精度で正弦信号の周波数を検出できることが確かめられた。
したがって、本発明の実施の形態例によると、入力信号から周波数や振幅が変動する正弦成分を抽出できるようになり、基本周波数の推定、雑音除去、調波信号の分析合成、音声や楽曲などの信号分離、波動の到来方向の推定などの様々な分野に適用することができる。
【０１５６】
［５．適用例］
図１１は、本発明の正弦成分抽出装置を音解析に適用した例である。この適用例は、本発明者が先に提案した特許文献１でも記載した例であるが、本発明にも適用が可能であるので、再度同じ例を記載する。
すなわち、図１１に示すように、マイクロフォン１と音入力装置２とを接続する。そして、音入力装置２は、マイクロフォン１から入力された信号による音信号を生成する。音入力装置２により生成された音信号は、音声解析装置３に供給される。
また、音入力装置２により生成された音信号は、正弦成分抽出装置４に供給される。この正弦成分抽出装置４は、図１に示す構成の装置であり、音信号に含まれる正弦成分として、音信号に含まれる正弦成分やその高調波の正弦成分を抽出する。そして、正弦成分抽出装置４で抽出された正弦成分の周波数，振幅，位相の情報が、音声解析装置３に供給される。音声解析装置３は、正弦成分抽出装置４から供給される正弦成分の周波数，振幅，位相の情報に基づいて、入力した音信号の解析を行い、解析結果を出力端子５に出力する。
【０１５７】
図１２は、本発明の正弦成分抽出装置を音源分離に適用した例である。
すなわち、図１２に示すように、音入力装置２は、マイクロフォン１から入力された信号に基づいて音信号を生成し、その生成した音信号を音源分離装置６に供給する。音源分離装置６は、マイクロフォン１に入力された音が、複数の正弦成分からなる音である場合に、その複数の音源成分を分離するものである。
また、音入力装置２が生成した音信号は、正弦成分抽出装置４に供給される。この正弦成分抽出装置４は、音信号に含まれる正弦成分を抽出し、抽出された正弦成分の周波数，振幅，位相の情報を音源分離装置６に供給する。音源分離装置６は、正弦成分抽出装置４から供給される正弦成分の周波数，振幅，位相の情報に基づいて、入力信号に含まれるそれぞれの音源からの正弦成分を分離し、個別の出力端子７，８に出力する。例えば入力信号が複数の楽曲の演奏信号である場合には、図１２の構成を利用して、１つ１つの楽器ごとの演奏信号を分離することができる。
【０１５８】
これら図１１，図１２は、正弦成分抽出装置の適用例のごく一部の例であり、その他の各種信号処理にも適用可能である。音信号や楽曲の演奏信号に適用することも一例であり、その他の信号に含まれる周波数成分の抽出に適用が可能である。
【０１５９】
また、図１に示した正弦成分抽出装置の構成は、好適な一例を示したものであり、本発明の正弦成分抽出装置及び抽出方法は、図１に示す構成のものに限定されない。例えば、実施の形態例のフーリエ変換処理を行うプログラムを実装したコンピュータ装置を用意して、そのコンピュータ装置での演算処理で、同様の正弦成分抽出処理を行うようにしてもよい。コンピュータ装置に実装するプログラムとしては、例えば図２に示した手順を実行するものとする。
【符号の説明】
【０１６０】
１…マイクロフォン、２…音入力装置、３…音声解析装置、４…正弦成分抽出装置、５，７，８…出力端子、６…音源分離装置、１１…信号入力端子、１２…窓関数入力端子、１３…乗算器、１４…離散フーリエ変換回路、１５…記憶部、１６…振幅比・位相差算出部、１７…検出部、１８…正弦成分統合部、１９…出力端子
【図１】
2018025628.tif 000082
【図２】
2018025628.tif 000083
【図３】
2018025628.tif 000084
【図４】
2018025628.tif 000085
【図５】
2018025628.tif 000086
【図６】
2018025628.tif 000087
【図７】
2018025628.tif 000088
【図８】
2018025628.tif 000089
【図９】
2018025628.tif 000090
【図１０】
2018025628.tif 000091
【図１１】
2018025628.tif 000092
【図１２】
2018025628.tif 000093