(57)【特許請求の範囲】
【請求項１】
入力された非同期信号の周波数解析を行うに際し、サンプリングされた解析対象信号と正弦波モデル信号との差の二乗和が最小値になるような正弦波モデル信号のパラメータである周波数ｆ’と位相φ’と振幅Ａ’についてそれぞれ初期値を求め、求めた各初期値から非線形方程式の最適解に収束させる信号解析方法であって、信号解析装置又はコンピュータを用いて、前記パラメータのうち周波数ｆ’と位相φ’について最急降下法を適用して収束させてから減速法を取り入れたニュートン法を適用して収束させ、その後、前記パラメータのうち振幅Ａ’について最急降下法を適用して収束させてから減速法を取り入れたニュートン法を適用して収束させることとし、
前記解析対象信号が１次元の計測値ｓ（ｘｎ）であり、データ個数をＮとすると、数式（１）の右辺が最小値になるような周波数ｆｘ’と位相φ’と振幅Ａ’を求めることを特徴とし、
前記解析対象信号が２次元の計測値ｓ（ｘｍ，ｙｎ）であり、各データ個数をＭ，Ｎとすると、数式（２）の右辺が最小値になるような周波数ｆｘ’，ｆｙ’と振幅Ａ’と位相φ’を求めることを特徴とし、
前記解析対象信号が３次元の計測値ｓ（ｘｌ，ｙｍ，ｚｎ）であり、各データ個数をＬ，Ｍ，Ｎとすると、数式（３）の右辺が最小値になるような周波数ｆｘ’，ｆｙ’，ｆｚ’と振幅Ａ’と位相φ’を求めることを特徴とし、
前記解析対象信号が４次元の計測値ｓ（ｘｋ，ｙｌ，ｚｍ，ｗｎ）であり、各データ個数をＫ，Ｌ，Ｍ，Ｎとすると、数式（４）の右辺が最小値になるような周波数ｆｘ’，ｆｙ’，ｆｚ’，ｆｗ’と振幅Ａ’と位相φ’を求めることを特徴とし、
前記解析対象信号がｎ次元（ｎは５以上の整数）の計測値ｓ（ｘ１ｎ１，ｘ２ｎ２，ｘ３ｎ３，・・・，ｘｎｎｎ）であり、各データ個数をＮ１，Ｎ２，Ｎ３，・・・，ＮＮとすると、数式（５）の右辺が最小値になるような周波数ｆｘ１’，ｆｘ２’，ｆｘ３’，・・・，ｆｘｎ’と振幅Ａ’と位相φ’を求めることを特徴とする信号解析方法。
【数１】
0005590547.tif 000028
【数２】
0005590547.tif 000029
【数３】
0005590547.tif 000030
【数４】
0005590547.tif 000031
【数５】
0005590547.tif 000032
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、入力された非同期信号の周波数解析を行う信号解析方法に関する。
【背景技術】
【０００２】
従来から、信号の計測、ディジタル信号の圧縮符号化、及び、将来における経済時系列信号変動の予測等、様々な分野において、入力された信号の周波数解析が行われている。いわゆる正弦波モデルを用いた信号解析方法としては、非特許文献１にも紹介されているように、様々な手法が提案されているが、その多くは、信号の周波数を事前に推定し、その推定した周波数における最適な振幅や初期位相を求めることによって行われている。【０００３】
信号は複数の周波数の正弦波の重ね合わせであり、信号を解析するには、当該信号を構成する様々な周波数の正弦波毎に振幅や初期位相等の周波数パラメータを推定することになる。このような周波数解析手法としては、フーリエ変換が広く知られている。フーリエ変換における周波数パラメータの推定は、周期的な信号の場合には、次式（１）に示すように、有限長の窓によるフーリエ変換式を用いて行うことができ、この次式（１）に基づいて、特定の周波数ｆ（＝ｎ／Ｔ：ｎは整数、Ｔは解析窓長）のパラメータ（振幅Ａ及び初期位相φ）を推定することができる。ただし、振幅Ａと初期位相φとＸ（ｆ）との関係は次式（２）である。【０００４】
【数１】
0005590547.tif 000002
【数２】
0005590547.tif 000003
【０００５】
このようなフーリエ変換によって推定することが可能な周波数パラメータは、解析窓長Ｔに依存し、１／Ｔの整数倍の周波数分しか計算することができない。したがって、上式（１）を用いてパラメータを推定している限り、周波数ｆは、原理的に離散的となり、非周期信号の周波数成分を特定することはできない。また、フーリエ変換を用いて周波数解析を行っている分野では、そのほとんどが、周期を持たない非周期信号を解析対象とするにもかかわらず、無限長の積分区間にわたる積分計算が極めて困難であることから、上式（１）を用いなければならないのが現状である。【０００６】
この問題を解決するために、窓関数の調整や解析区間を複数利用して見かけ上の周波数分解能を向上させる試みがある。【０００７】
例えば、非特許文献２に記載の短時間フーリエ変換（Short-Time Fourier Transform、Short-Term Fourier Transform；ＳＴＦＴ）は、音響信号を短時間の解析窓によって複数の信号に分割し、各信号に対して独立にフーリエ変換を行う手法である。【０００８】
また、非特許文献３乃至非特許文献６に記載のＡＢＳ（Analysis by Synthesis）法やＧＨＡ（Generalized Harmonic Analysis）等においては、最適な周波数を定めるための周波数を複数用意して振幅及び初期位相を求め、その中から最適な周波数と振幅と初期位相との組み合わせを選択している。【０００９】
【非特許文献１】東山三樹夫、小池恒彦、「高い周波数分析精度の信号分析手法」、日本音響学会論文誌、１９８８年、第５４巻、第８号
【非特許文献２】L. R. Rabiner and R. W. Schafer,“Digital Processing of Speech Signals”, Prentice-Hall, Englewood Cliffs, NJ, 1978年
【非特許文献３】T. F. Quatieriand R. J. McAulay, “Speechtransformations based on a sinesoidal representation”, IEEE Trans. Acoust. Speech Signal Process. ASSP 34, 1986年, p.1449-1464
【非特許文献４】E. B. Georgeand M. J. Smith, “Analysis-by-synthesis/overlap-addsinusoidal representation”, J. Audio Eng. Soc.40, 1992年, p.497-515
【非特許文献５】E. B. George and M. J. Smith, “Speech analysis/synthesis and modification using ananalysis-by-synthesis/overlap-add sinusoidal model”, IEEE Trans. Speech Audio Process. 5, 1997年, p.389-406
【非特許文献６】T. Terada, H.Nakajima, M. Tohyama and Y. Hirata, “Non-stationary waveform analysis and synthesis using generalizedharmonic analysis”, IEEE-SP, Int. Symp. TF/TSAnalysis, 1994年, p.429-432
【００１０】
しかしながら、上述した非特許文献２に記載の短時間フーリエ変換においては、短時間の解析窓長の逆数の整数倍の周波数でしか、周波数パラメータを推定することができないという問題がなおも存在する。【００１１】
また、上述した非特許文献３乃至非特許文献６に記載のＡＢＳ法やＧＨＡ等においては、元の信号に対して、予め用意された周波数と異なる周波数が存在する場合には、誤った組み合わせが複数検波されてしまい、解析精度が低下するという問題がある。【００１２】
これら窓関数の調整や解析区間を複数利用して見かけ上の周波数分解能を向上させる試みは、
１）周波数分解能を上げるためには解析窓長を長くしなければならないが、時間的な周波数変動をともなう信号では周波数が平均化され、時間的な変動を捉えることが困難となること
２）一方、解析窓長を短くすると、周波数分解能が低下すること
という問題が指摘されている。【００１３】
このように、周波数の変化を細かく捉えることと、周波数分解能を上げることとは、互いに相反する関係にあり、信号解析上の大きな問題となっていた。
【発明の開示】
【００１４】
本発明は、このような実情に鑑みてなされたものであり、周波数分解能が解析窓長に依存せずに高い周波数分解能を有し、驚くべき高精度に周波数解析を行うことができる信号解析方法を提供することを目的とする。
【００１５】
本発明の信号解析方法は、入力された非同期信号の周波数解析を行うに際し、サンプリングされた解析対象信号と正弦波モデル信号との差の二乗和が最小値になるような正弦波モデル信号のパラメータである周波数ｆ’と位相φ’と振幅Ａ’についてそれぞれ初期値を求め、求めた各初期値から非線形方程式の最適解に収束させる信号解析方法であって、信号解析装置又はコンピュータを用いて、前記パラメータのうち周波数ｆ’と位相φ’について最急降下法を適用して収束させてから減速法を取り入れたニュートン法を適用して収束させ、その後、前記パラメータのうち振幅Ａ’について最急降下法を適用して収束させてから減速法を取り入れたニュートン法を適用して収束させることとし、
前記解析対象信号が１次元の計測値ｓ（ｘｎ）であり、データ個数をＮとすると、数式（１）の右辺が最小値になるような周波数ｆｘ’と位相φ’と振幅Ａ’を求めることを特徴とし、
前記解析対象信号が２次元の計測値ｓ（ｘｍ，ｙｎ）であり、各データ個数をＭ，Ｎとすると、数式（２）の右辺が最小値になるような周波数ｆｘ’，ｆｙ’と振幅Ａ’と位相φ’を求めることを特徴とし、
前記解析対象信号が３次元の計測値ｓ（ｘｌ，ｙｍ，ｚｎ）であり、各データ個数をＬ，Ｍ，Ｎとすると、数式（３）の右辺が最小値になるような周波数ｆｘ’，ｆｙ’，ｆｚ’と振幅Ａ’と位相φ’を求めることを特徴とし、
前記解析対象信号が４次元の計測値ｓ（ｘｋ，ｙｌ，ｚｍ，ｗｎ）であり、各データ個数をＫ，Ｌ，Ｍ，Ｎとすると、数式（４）の右辺が最小値になるような周波数ｆｘ’，ｆｙ’，ｆｚ’，ｆｗ’と振幅Ａ’と位相φ’を求めることを特徴とし、
前記解析対象信号がｎ次元（ｎは５以上の整数）の計測値ｓ（ｘ１ｎ１，ｘ２ｎ２，ｘ３ｎ３，・・・，ｘｎｎｎ）であり、各データ個数をＮ１，Ｎ２，Ｎ３，・・・，ＮＮとすると、数式（５）の右辺が最小値になるような周波数ｆｘ１’，ｆｘ２’，ｆｘ３’，・・・，ｆｘｎ’と振幅Ａ’と位相φ’を求めることを特徴とする信号解析方法。
【数１】
0005590547.tif 000004
【数２】
0005590547.tif 000005
【数３】
0005590547.tif 000006
【数４】
0005590547.tif 000007
【数５】
0005590547.tif 000008
【００１８】
このような本発明にかかる信号解析方法、信号解析装置、及び信号解析プログラムが実装された装置は、非周期信号を想定した周波数解析を目的として定式化を行い、非周期信号のフーリエ変換式のパラメータを求める問題を非線形方程式の最適解を求める問題に置き換えることにより、周波数分解能が解析窓長に依存しない新たな周波数解析手法を実現することができる。
【００１９】
以上のような本発明においては、解析窓長と周波数分解能の制約を受けることなく高い周波数分解能を有し、従来の周波数解析手法における精度と比べて１０万〜１００億倍以上という驚くべき高精度に周波数解析を行うことができる。
【図面の簡単な説明】
【００２０】
【図１】本発明の実施の形態として示す信号解析装置の構成を示すブロック図である。
【図２】本周波数解析手法とＤＦＴとの違いを説明するための図である。
【図３】本周波数解析手法とＤＦＴとの違いを説明するための図であり、２次元の解析対象信号に各手法を適用した場合の具体例を示す図である。
【図４】本周波数解析手法とＤＦＴとＧＨＡとの違いを説明するための図であり、各手法の誤差を求めた結果を示す図である。
【図５】解析対象信号としてのピアノ音の具体例を示す図である。
【図６（Ａ）】ＤＦＴを用いて図５に示す信号を解析して得られたスペクトログラムを示す図である。
【図６（Ｂ）】本周波数解析手法を用いて図５に示す信号を解析して得られたスペクトログラムを示す図である。
【図７】本周波数解析手法とＤＦＴとの違いを説明するための図であり、２次元の解析対象信号に各手法を適用して得られるスペクトルの解析精度の具体例を示す図である。
【図８】本周波数解析手法とＪＰＥＧとＪＰＥＧ２０００との違いを説明するための図であり、２次元のビットマップ画像データを圧縮した場合の具体例を示す図である。
【図９】本周波数解析手法を適用した経済指標の予測例を示す図である。
【図１０】本周波数解析手法を適用して求めた予測データと実際のデータとの誤差の分布を示す図である。
【発明を実施するための最良の形態】
【００２１】
以下、本発明を適用した具体的な実施の形態について図面を参照しながら詳細に説明する。【００２２】
この実施の形態は、入力された信号の周波数解析を行う信号解析装置である。特に、この信号解析装置は、非線形方程式を解くことでフーリエ係数を推定することにより、周波数分解能が解析窓長に依存しない新たな周波数解析手法を実現するものである。【００２３】
［信号解析装置の構成］
信号解析装置は、例えばコンピュータ等から構成され、図１に示すように、各部を統括的に制御するＣＰＵ（Central Processing Unit）１１と、各種プログラムを含む各種情報を格納する読み取り専用のＲＯＭ（Read Only Memory）１２と、ワークエリアとして機能するＲＡＭ（Random Access Memory）１３と、各種情報を読み出し及び／又は書き込み可能に記憶する記憶部１４と、ユーザインターフェースとしての図示しない所定の操作デバイスを介した入力操作の処理及び制御を行う入力操作制御部１５と、各種情報を表示する表示部１６とを備える。【００２４】
ＣＰＵ１１は、記憶部１４等に格納されている各種アプリケーションプログラムをはじめとする各種プログラムを実行し、各部を統括的に制御する。【００２５】
ＲＯＭ１２は、各種プログラムをはじめとする各種情報を格納している。このＲＯＭ１２に格納されている情報は、ＣＰＵ１１の制御のもとに読み出される。【００２６】
ＲＡＭ１３は、ＣＰＵ１１が各種プログラムを実行する際のワークエリアとして機能し、ＣＰＵ１１の制御のもとに、各種情報を一時記憶するとともに、記憶している各種情報を読み出す。【００２７】
記憶部１４は、本発明にかかる信号解析プログラム等のアプリケーションプログラムの他、周波数解析の対象となる解析対象信号のデータをはじめとする各種情報を記憶する。この記憶部１４としては、例えば、ハードディスクや不揮発性メモリ等を用いることができる。また、記憶部１４には、本体に対して着脱可能とされるフレキシブルディスクやメモリカード等の記憶媒体に対して、各種情報の読み出し及び／又は書き込みを行うドライブ装置も含まれる。この記憶部１４に記憶されている各種情報は、ＣＰＵ１１の制御のもとに読み出される。【００２８】
入力操作制御部１５は、例えば、キーボード、マウス、キーパッド、赤外線リモートコントローラ、スティックキー、又はプッシュボタンといった、ユーザインターフェースとしての図示しない所定の操作デバイスを介した入力操作を受け付け、操作内容を示す制御信号をＣＰＵ１１に対して供給する。【００２９】
表示部１６は、例えば、液晶ディスプレイ（Liquid Crystal Display；ＬＣＤ）、プラズマ・ディスプレイ・パネル（Plasma Display Panel；ＰＤＰ）、有機エレクトロルミネッセンス（Organic ElectroLuminescent）ディスプレイ、又はＣＲＴ（Cathode Ray Tube）といった、各種表示デバイスであり、ＣＰＵ１１の制御のもとに各種情報を表示する。例えば、表示部１６は、ＣＰＵ１１によって信号解析プログラムが起動されると、その画面を表示し、入力された解析対象信号や周波数解析結果等を表示する。【００３０】
このような各部を備える信号解析装置は、ＣＰＵ１１の制御のもとに、信号解析プログラムを実行すると、ＣＰＵ１１の制御のもとに、入力された信号の周波数解析を行う。なお、周波数解析の対象となる信号、すなわち、上述した解析対象信号は、図示しない信号入力部を介して入力され、ＣＰＵ１１によって読み出し可能に記憶部１４に記憶される。例えば、信号解析装置は、任意のセンサによって計測された時系列信号を解析対象信号とする場合には、当該信号解析装置とセンサとを接続する所定のインターフェースを介して解析対象信号としての時系列信号を記憶部１４に記憶させることによって入力する。また、信号解析装置は、ディジタルカメラ等の撮像装置によって撮像された画像データの圧縮符号化を行う機能を有し、その圧縮符号化の際に周波数解析を行う場合には、当該信号解析装置と撮像装置とを接続する所定のインターフェースを介して解析対象信号としての画像データを記憶部１４に記憶させることによって入力する。さらに、信号解析装置は、株価等の将来における変動を予測するために用いる場合には、インターネット等を介して解析対象信号としての過去の株価データを記憶部１４に記憶させることによって入力する。その他、信号解析装置は、ユーザが作成した任意のデータを解析対象信号として記憶部１４に記憶させることによって入力するようにしてもよい。すなわち、信号入力部は、解析対象信号をＣＰＵ１１によって読み出し可能に入力させる機能を有する部位である。なお、信号入力部は、アナログ信号を入力した場合には、Ａ／Ｄ変換を行ってディジタル信号に変換する機能をあわせ持つことはいうまでもない。このとき、信号入力部は、必要に応じてアンチエイリアシングフィルタを含むＡ／Ｄ変換器としてもよい。信号解析装置は、ＣＰＵ１１の制御のもとに、このようにして入力された解析対象信号の周波数解析を行い、その解析結果等を、図示しない解析結果出力部を介して表示部１６に表示出力させたり、印刷装置やその他の機器に出力したりする。【００３１】
［周波数解析アルゴリズム］
以下、このような信号解析装置における周波数解析アルゴリズムについて詳述する。【００３２】
信号解析装置に適用する周波数解析手法、すなわち、本発明にて新たに提案する周波数解析手法（以下、本周波数解析手法という。）は、「周波数の変化を細かく捉えることと、周波数分解能を上げることとは、互いに相反する関係にある」という従来の周波数解析手法の問題を本質的に解決するために、非線形方程式の解を求める問題に着目した。すなわち、本周波数解析手法においては、次式（３）に示す非周期信号のフーリエ変換式の周波数パラメータを求める問題を非線形方程式の最適解を求める問題に置き換えた。【００３３】
【数３】
0005590547.tif 000009
【００３４】
具体的には、本周波数解析手法においては、次式（４）に示すように、解析対象信号ｘ（ｎ）と正弦波モデル信号との差の二乗和で表される非線形方程式の最適解として、この非線形方程式の右辺が最小値になるような周波数ｆ’、振幅Ａ’、及び初期位相φ’を求める。なお、次式（４）において、Ｌはフレーム長（解析窓長）であり、ｆｓはサンプリング周波数［Ｈｚ］である。このように、非線形方程式に周波数ｆ’をパラメータとして盛り込んだことは、いまだ事例がなく斬新である。換言すれば、本周波数解析手法においては、振幅Ａ’及び初期位相φ’のみならず、周波数ｆ’をも正確に求めることを目的とし、上式（３）に対して非線形方程式の解法を用いるものである。本周波数解析手法においては、このような最小二乗法によって非線形方程式の最適解を求める問題に帰着させることにより、解析窓の影響やエイリアシングの影響がなくなり、解析窓長が、１周期未満であってもよく、周期の整数倍でなくてもよく、さらには、不等間隔であってもよい等、柔軟な周波数解析処理を実現することが可能となる。【００３５】
【数４】
0005590547.tif 000010
【００３６】
なお、上式（４）は、解析対象信号ｘ（ｎ）が、センサによって計測された時系列信号等の１次元信号である場合を想定しているが、本周波数解析手法は、２次元の画像データ等の２次元信号や動画像データ等の３次元信号、ホログラムのような立体動画像データ等の４次元信号、さらにはそれ以上の次元の信号に対しても適用可能である。すなわち、本周波数解析手法は、１次元の解析対象信号をｓ（ｘｎ）、２次元の解析対象信号をｓ（ｘｍ，ｙｎ）、３次元の解析対象信号をｓ（ｘｌ，ｙｍ，ｚｎ）、４次元の解析対象信号をｓ（ｘｋ，ｙｌ，ｚｍ，ｗｎ）として一般化すると、解析対象信号が１次元、２次元、３次元、又は４次元の場合のそれぞれについて、次式（５）乃至次式（８）に示す非線形方程式の最適解を求めることになる。【００３７】
【数５】
0005590547.tif 000011
【数６】
0005590547.tif 000012
【数７】
0005590547.tif 000013
【数８】
0005590547.tif 000014
【００３８】
換言すれば、本周波数解析手法は、任意のｎ次元（ｎは１以上の整数）の解析対象信号をｓ（ｘ１ｎ１，ｘ２ｎ２，ｘ３ｎ３，・・・，ｘｎｎｎ）として一般化すると、次式（９）に示す非線形方程式の最適解を求めることになる。
【数９】
0005590547.tif 000015
【００３９】
なお、以下では、説明の便宜上、解析対象信号が１次元であり、上式（４）に示す非線形方程式の最適解を求めるものとする。【００４０】
さて、上式（４）に示す非線形方程式の最適解を実際に求めるにあたっては、以下のような方法をとることができる。【００４１】
本周波数解析手法においては、振幅Ａ’、周波数ｆ’、及び初期位相φ’のそれぞれについて適切な初期値を求め、これら初期値から非線形方程式の解法を用いて最適解に収束させる。この非線形問題では、上式（４）をコスト関数とする最小化問題とする。なお、適切な初期値は、離散フーリエ変換（Discrete Fourier Transform；以下、ＤＦＴという。）やウェーブレット変換等の任意の周波数変換を行ったり、フィルタリングを行うことによっておおよその見当をつけたりする等、既存の任意の方法を適用して求めることができる。【００４２】
まず、本周波数解析手法においては、上式（４）における正弦波モデル信号の位相を構成する周波数パラメータｆ’，φ’について、いわゆる最急降下法を適用し、周波数パラメータｆｍ’，φｍ’を次式（１０）及び次式（１１）によって求める。【００４３】
【数１０】
0005590547.tif 000016
【数１１】
0005590547.tif 000017
【００４４】
なお、上式（１０）及び上式（１１）においては、次式（１２）と略している。また、μｍは、いわゆる減速法に基づく重み係数であり、各漸化式によって求められるコスト関数を単調減少数列にするために、適時０〜１の値をとる。【００４５】
【数１２】
0005590547.tif 000018
【００４６】
周波数パラメータｆｍ’，φｍ’を求めることができれば、上式（４）における正弦波モデル信号の係数としての周波数パラメータＡ’を一意に求めることができるため、本周波数解析手法においては、次式（１３）によって周波数パラメータＡｍ’を収束させる。【００４７】
【数１３】
0005590547.tif 000019
【００４８】
本周波数解析手法においては、これら一連の計算を反復して行うことにより、振幅Ａ’、周波数ｆ’、及び初期位相φ’を高精度に収束させることができる。特に、本周波数解析手法においては、上式（４）における正弦波モデル信号の位相を構成する周波数パラメータｆ’，φ’と、係数としての周波数パラメータＡ’とを別個に求めることにより、計算を簡便に行うことができる。【００４９】
しかしながら、最急降下法は、比較的広い範囲から収束するものの、１回の反復では精度が低く、収束するまでに時間を要する。【００５０】
そこで、本周波数解析手法においては、最急降下法を適用して周波数パラメータｆｍ’，φｍ’をある程度まで収束させた後、さらに、いわゆるニュートン法を適用して高精度に収束させるのが望ましい。具体的には、本周波数解析手法においては、ニュートン法として、次式（１４）及び次式（１５）に示す漸化式によって周波数パラメータｆｍ’，φｍ’を求める。【００５１】
【数１４】
0005590547.tif 000020
【数１５】
0005590547.tif 000021
【００５２】
ただし、上式（１４）及び上式（１５）において、Ｊは次式（１６）とし、次式（１７）と略している。また、νｍもμｍと同様に減速法に基づく重み係数であり、適時０〜１の値をとる。【００５３】
【数１６】
0005590547.tif 000022
【数１７】
0005590547.tif 000023
【００５４】
本周波数解析手法においては、上式（１４）及び上式（１５）によって周波数パラメータｆｍ’，φｍ’を求めた後、最急降下法と同様に、上式（１３）によって周波数パラメータＡｍ’を収束させ、この一連の計算をさらに反復して行う。【００５５】
このように、本周波数解析手法においては、最急降下法とニュートン法とを組み合わせたハイブリッド型の解法を用いることにより、高速に且つ高精度に周波数パラメータＡ’，ｆ’，φ’を推定することができる。【００５６】
また、本周波数解析手法においては、解析対象信号ｘ（ｎ）が複合正弦波の場合であっても、逐次減算処理することにより、近似的にスペクトルパラメータを導出することができる。ここで、解析対象信号ｘ（ｎ）が複数の正弦波の和であり、次式（１８）のように表されているとする。【００５７】
【数１８】
0005590547.tif 000024
【００５８】
パーセヴァル（Parseval）の定理より、解析対象信号ｘ（ｎ）の周波数ｆｋと正弦波モデル信号の周波数パラメータｆ’とが全く一致しない場合、すなわち、次式（１９）である場合には、上式（４）に示す非線形方程式は次式（２０）となる。また、周波数パラメータｆ’，φ’の組が、周波数ｆｋ及び初期位相φｋの組のいずれかに一致する場合には、上式（４）に示す非線形方程式は次式（２１）となる。さらに、振幅Ａｊが周波数パラメータＡ’とも一致した場合には、解析対象信号から推定スペクトルに関する周波数成分を完全に消去することができる。そのため、最適解を求める問題は、周波数に対して独立であり、解析対象信号から順次個別に推定すれば、複数の正弦波で表される信号にも応用することができる。【００５９】
【数１９】
0005590547.tif 000025
【数２０】
0005590547.tif 000026
【数２１】
0005590547.tif 000027
【００６０】
すなわち、本周波数解析手法においては、解析対象信号ｘ（ｎ）が複合正弦波の場合であっても、逐次残差信号に対して同様に処理を行い、複数の正弦波を抽出することができる。【００６１】
ただし、複数のスペクトルの場合には、上式（３）からもわかるように、無限長を仮定しているため、複数のスペクトルの周波数同士が近い場合には、合成スペクトルの周期が解析窓長よりも数段長くなり、上式（１７）が満たされなくなるため、誤差が発生する場合もある。【００６２】
音声信号や音響信号等の信号を複合正弦波によって表現するためには、これまで多くのスペクトル数（正弦波の数）が必要であったが、本周波数解析手法においては、そのような信号であっても僅かなスペクトル数で誤差なく表現することができる。すなわち、信号をより少ないスペクトル数で表現可能であることは、情報圧縮の用途に有効であることは勿論のこと、それ以外にも、物理現象に関して波数表現をするようなアプローチに用いることにより、信号の本質的な特性も見抜くことができることを示している。このように、本周波数解析手法は、複合正弦波の解析対象信号に対しても有効であり、様々な分野に応用することが可能である。【００６３】
［本周波数解析手法の有効性］
以下、本周波数解析手法の有効性について具体的に説明する。【００６４】
従来、センサによって計測された信号を、計算機を用いて解析するためには、一般に、例えば図２に示すように、その信号の一部を切り出し、ＤＦＴといった代表的な調和型の周波数解析手法によって解析していた。【００６５】
しかしながら、解析対象信号を切り出して（窓掛けして）調和型の解析を行うことは、当該解析対象信号が、切り出した一部の信号が完全周期で繰り返される信号であると擬制して解析することに他ならず、当然ながらその周波数特性も、元のスペクトル特性とは完全に一致しないものとなる。また、この信号切り出しの影響を回避するために、窓関数によって影響を軽減することも試みられているが、切り出した一部の信号を完全周期として解析することには変わりなく、本質的な問題解決には至っていないのが現状である。【００６６】
これに対して、本周波数解析手法は、非周期信号を想定した周波数解析を目的として定式化を行っていることから、解析窓長と周波数分解能の制約を受けることなく、正確に周波数とそのパラメータとを推定することができる。【００６７】
また、本周波数解析手法は、周波数の推定についても、従来の離散的な探索から、動的且つ連続的な探索となるため、格段に精度を向上させることができ、革新的な手法である。特に、本周波数解析手法は、系が安定している場合には、現象の一部を観測して解析することによって正確な特性を推定することができるため、予測問題への応用が期待できる。【００６８】
さらに、多次元化に関しても、従来の周波数解析手法においては、解析対象信号が１次元である場合と同様に窓関数の影響を受ける。具体的には、例えば図３に示すように、２次元の解析対象信号（原画像データ）のうち図３左図中の破線部を切り出してＤＦＴを施して解析したとしても、その切り出した区間を完全周期とした画像が延々と繰り返されるものを解析したことと等価であることから、右上図に示すように、復元した画像データが窓関数の影響を受けたものとなり、本質的な特性を見出すことは困難である。【００６９】
これに対して、上式（６）を適用した本周波数解析手法においては、図３右下図に示すように、窓関数の影響が軽減され、本質的な周波数特性を見出し、切り出した区間の周囲の画像も完全に再現することができる。【００７０】
以上のように、本周波数解析手法は、非線形方程式の最適解を求めることにより、正弦波モデル信号の周波数ｆ’、振幅Ａ’、及び初期位相φ’を高速に且つ高精度に求めることができる。具体的な精度を立証するために、本願発明者は、ＤＦＴと、ＤＦＴの発展型のうち最も解析精度が高いといわれているＧＨＡ（Generalized Harmonic Analysis）とを比較対象として精度の検証を行った。【００７１】
なお、ＤＦＴやＧＨＡは、１つの解析窓長に見かけ上複数の窓長を持たせていることから、周波数分解能が解析窓長に依存するが、その分解周波数が有限長であり、解析対象信号の周波数が分解周波数以外の周波数となった場合には解析することができず、解析対象信号が正確に解析できる周波数と異なる場合には、最も近い分解周波数の他に、その周辺に小さなスペクトルの周波数（側帯波成分）が現れ、複数の周波数が出現してしまう。【００７２】
このような現象が本周波数解析手法においても生じるか否かについて、すなわち、本周波数解析手法の周波数分解能を検証するために、解析窓長を１秒（１０２４サンプル）とした１次元の非常に短い単一正弦波を解析し、各手法によって正弦波を１本抽出して元の信号との二乗誤差を調べた。その結果を図４に示す。【００７３】
図４に示すように、ＤＦＴにおいては、基本周波数の整数倍以外の周波数における解析精度の悪化がみられた。また、ＧＨＡにおいては、１Ｈｚ以上の周波数ではＤＦＴと比べて２〜５桁程度の精度向上がみられた。これに対して、本周波数解析手法においては、１Ｈｚ以上の周波数ではＤＦＴと比べて１０桁以上、ＧＨＡと比べて５桁以上の精度向上がみられた。すなわち、本周波数解析手法は、既存の周波数解析手法と比べて１０万〜１００億倍以上の精度向上がみられた。【００７４】
特に、１Ｈｚ以下の周波数を正確に推定することができるということは、解析窓長を超えた長い周期信号であっても解析可能であることを示しており、株価等の様々な変動要因が含まれる信号に対しても、そのスペクトル構造をより正確に推定可能であることを示している。【００７５】
このように、本周波数解析手法は、最も解析精度が高いといわれているＧＨＡと比べても驚くべき高精度に解析を行うことができる。国内外の様々な研究事例をみても、これほどの精度で信号を解析できる手法は存在せず、本周波数解析手法は、今後、正確な解析が必要とされる様々な分野への応用が期待できる手法であるといえる。【００７６】
［本周波数解析手法の具体的な応用例］
現在、様々な分野において周波数解析が利用されているが、本周波数解析手法は、これらの分野以外へも波及効果が期待できると考えられる。本周波数解析手法は、フーリエ係数パラメータを、フーリエ変換を用いずに高精度に求めることができることから、例えば周期性（波動性）をともなう音響・振動以外のこれまで様々な非定常的な変動として考えられていた現象を、周期性をともなう信号としてモデル化できる可能性もある。以下、本願発明者が行った本周波数解析手法の応用例について説明する。【００７７】
［計測分野への応用例］
近年、ピアノ弦がハンマーで叩かれた場合には、そのピアノ弦は、垂直な振動のみならず水平な振動も加わり、楕円軌道上で回転運動を行っていることが、ピアノ弦の近傍にレーザ計測機器を設置して当該ピアノ弦との距離変位を計測することにより、明らかになっている。また、近年では、１本のピアノ弦は、２つの基本周波数を有し、それらの基本周波数で交互に揺動しながら連続的に変化していることが光学的に実証されている。【００７８】
しかしながら、ピアノ音を計測してその音の周波数解析を行うことにより、かかるピアノ弦の基本周波数での小さな揺れを正確に捉えることは困難である。従来の周波数解析は、ＤＦＴや高速フーリエ変換（Fast Fourier Transform；ＦＦＴ）に代表されるスペクトラムアナライザーを用いて行われることが多かったが、一般に、窓関数の影響や周波数分解能の低さに起因して、ピアノ弦の小さな揺れを解析することはできなかった。このように、音のみの解析については、光学的な計測と同等な成功事例は未だ報告されていない。【００７９】
そこで、本願発明者は、本周波数解析手法を用いて、一般的に用いられているＣＤ（Compact Disc）程度の音質の観測音のみから、ピアノ弦のスペクトル特性の推定やそのスペクトル周辺の揺らぎを可視化することを試みた。【００８０】
具体的には、図５に示すように、ピアノ調律作業中における音名Ａ４（階名ラ、４４０Ｈｚ）のピアノ音を、ＣＤと同じサンプリング周波数４４．１ｋＨｚで録音した信号を解析対象信号とし、この解析対象信号をＤＦＴ及び本周波数解析手法のそれぞれによって解析し、そのスペクトログラムを評価した。解析区間は１秒であり、窓長のデータ点数を４０９６点（０．０９３ミリ秒）とし、シフト長を窓長の１／１６とした。この窓長を用いた場合、ＤＦＴの周波数分解幅は約１０．７６Ｈｚとなる。【００８１】
解析結果を図６(Ａ)及び図６(Ｂ)に示す。図６(Ａ)は、ＤＦＴによる解析結果であり、基音付近を拡大したものである。また、図６(Ｂ)は、本周波数解析手法による解析結果であり、同様に基音付近を拡大したものである。なお、これら図６(Ａ)及び図６(Ｂ)において、右側横軸は時間を示し、左側横軸は周波数を示し、縦軸は振幅を示している。【００８２】
図６(Ａ)においては上述した約１０．７６Ｈｚの周波数分解幅でスペクトルを表示している。この図６(Ａ)をみると、４４０Ｈｚ近傍に現れているスペクトルが最も大きいため、その周辺に正しいスペクトルがあることが窺える。また、４３０Ｈｚ近傍に波打ったようなエネルギの増減を呈しているスペクトルも確認できる。今回の調律時には、打鍵時以降に外力が作用しないため、このようなエネルギの増減は起こり得ないと考えると、４４０Ｈｚ近傍のスペクトル列周辺で周波数が揺れているものと予想できるが、周波数分解幅以下の小さな変動を解析することができず、正確にどのような現象が起きているのかを把握することができない。【００８３】
これに対して、図６(Ｂ)をみると、４４０Ｈｚ近傍に現れているスペクトルが連続的に揺らいでいる様子が明白に把握できる。これは、ピアノ弦が複数の楕円軌道上で回転運動を行っていることによる周波数の変化を捉えたものに他ならず、ＤＦＴでは捉えることができなかった周波数の揺れが確認できる。また、本周波数解析手法は、ＤＦＴと比べて、同じ解析窓長であっても正確に周波数（４４０Ｈｚ）を推定していることがわかる。【００８４】
このように、ＤＦＴにおいては、解析窓の影響のために周波数の揺れが存在する現象をスペクトルの単調減少関数で表現することができないが、本周波数解析手法においては、周波数変動をともなった単調減少スペクトルで表現することができる。そのため、本周波数解析手法は、複雑な周波数変動をともなう単調減少スペクトルの解析に有効である。【００８５】
この実験からも明らかなように、本周波数解析手法は、ＣＤ程度の音質の音の解析から、光学的計測と同等の成果を得ることができる画期的な手法であるといえる。【００８６】
［画像の符号化］
本周波数解析手法においては、窓の影響が少ないため、スペクトルに側帯波成分が全く出ない。すなわち、本周波数解析手法は、上述したように、複雑な解析対象信号であっても、僅かなスペクトル数で効率的に誤差なく表現することができる。【００８７】
図７に２次元信号をＤＦＴと本周波数解析手法とのそれぞれによって解析した具体例を示す。図７左図に示すように、本来であれば１本のスペクトルで表現される信号をＤＦＴによって変換した場合には、一般に、変換して得られるスペクトルも、図３右下図に示したような解析窓の影響に起因して、図７右上図に示すように、複数の側帯波成分が群れて現れてしまう。これは、画像の圧縮符号化において広く用いられている離散コサイン変換（Discrete Cosine Transform；ＤＣＴ）等においても同様である。これに対して、本周波数解析手法は、図７右下図に示すように、元の信号のスペクトルを正確に推定することができる。正確に解析できると、画像データを表現するスペクトル数を格段に削減できる可能性があり、画質を落とさずに飛躍的な情報圧縮効果を期待することができる。【００８８】
実際に、画像圧縮の一般的手法として広く知られているＪＰＥＧ（Joint Photographic Experts Group）とＪＰＥＧ２０００とを比較対象として、図８中央図に示す２次元のビットマップ画像データをソフトウェア圧縮した。なお、ＪＰＥＧは、ＤＣＴを基礎解析に利用している手法であり、ＪＰＥＧ２０００は、ウェーブレット変換を基礎解析に利用している手法であることは周知である。【００８９】
各周波数解析手法を適用して図８中央図に示す画像データを圧縮した後、同画像データに写り込んでいる背景部分を切り出した結果を図８左図に示し、同画像データに写り込んでいる人物の衣服部分を切り出した結果を図８右図に示す。すなわち、図８左図に示す結果は、変化が緩い部分画像（低周波数信号）の圧縮結果であり、図８右図に示す結果は、変化が激しい部分画像（高周波数信号）の圧縮結果である。【００９０】
図８左図をみると、５１７６バイトの原部分画像を圧縮した結果、ＪＰＥＧにおいては、７４０バイトとなり、原部分画像のサイズの１４．３％のサイズに圧縮された。また、ＪＰＥＧ２０００においては、３０６バイトとなり、原部分画像のサイズの５．９１％のサイズに圧縮された。これに対して、本周波数解析手法においては、４バイトとなり、原部分画像のサイズの０．０８％のサイズにまで圧縮された。なお、この本周波数解析手法による圧縮後の画像から原部分画像を復元するのに必要なスペクトル数は、たったの１本であった。【００９１】
一方、図８右図をみると、５１７６バイトの原部分画像を圧縮した結果、ＪＰＥＧにおいては、８２０バイトとなり、原部分画像のサイズの１５．４９％のサイズに圧縮された。また、ＪＰＥＧ２０００においては、３０７バイトとなり、原部分画像のサイズの５．９３％のサイズに圧縮された。これに対して、本周波数解析手法においては、２０バイトとなり、原部分画像のサイズの０．３９％のサイズにまで圧縮された。なお、この本周波数解析手法による圧縮後の画像から原部分画像を復元するのに必要なスペクトル数は、５本であった。【００９２】
このように、本周波数解析手法は、ＪＰＥＧやＪＰＥＧ２０００等の一般的な手法と比べて、圧縮後の情報量が約１／１００〜１／１０となり、他を遥かに凌ぐ圧縮率を実現することができる。また、本周波数解析手法は、低周波数信号から高周波数信号まで幅広く高精度に信号の特徴を特定することができる手法でもある。特に、本周波数解析手法においては、変化が緩い信号の解析を行う場合には、解析窓を超えた周期についても特定することができるため、背景等の画像を効率よく符号化することができる。したがって、本周波数解析手法は、ディジタル放送に移行しているテレビジョン放送に本周波数解析手法を適用することにより、携帯端末機を対象としたいわゆるワンセグであっても、ハイビジョン放送と同等若しくはそれ以上の画像品質を提供できる可能性がある。このように、次世代の画像／動画像符号化技術におけるコア技術となり得ることは極めて有意である。【００９３】
［経済時系列の予測］
大気圧等の平衡状態を保たせる作用のある状況下で、空間上の点において気圧の変化が生じると、その圧力変化が波動性をともなって周囲に伝搬する。かかる圧力を商品の売買における圧力に置き換え、経済指標も常に平衡状態を保たせる作用が働くと仮定すると、波動現象のように考えることができる。このような考え方に基づくと、株価等に代表される経済時系列等にも周期的な変動要因を仮定することができ、現実に、エリオット波動にともなう経済時系列の周期性に着目した研究も行われている。【００９４】
しかしながら、経済時系列の周期性に関しては、様々な周期性が複合的に組み合わされたモデルによって表されていると考えられている。しかも、実際の市場においては、その各周期性も時間的に変動する複雑なモデルとなる。【００９５】
そこで、本願発明者は、かかる困難な問題に対しても、本周波数解析手法を応用することを試みた。特に、本周波数解析手法は、解析窓への依存性が小さく、高い周波数分解能を有することから、解析対象信号の一部の区間を切り出して解析し、その後、延長描写を行うことにより、容易に慣性的で長期的な予測が可能となる。【００９６】
過去４年間（１９９７年７月〜２００１年６月）における実際の日経平均株価終値の推移を、図９中、太実線で示す。この株価データのうち、過去２年間（１９９７年７月〜１９９９年６月）を学習期間とし、この学習期間のデータを解析対象信号として本周波数解析手法を適用した。そして、学習期間以降の２年間（１９９９年７月〜２００１年６月）の株価がどのように推移するかを推定した。具体的には、過去２年間の学習期間の株価データを解析対象信号として本周波数解析手法を適用し、複数の正弦波を決定する周波数と振幅と初期位相とを求め、これらのパラメータが暫定的に変化しないものとして、それ以降２年間のデータを延長描写することにより、予測データを求めた。その結果を、図９中、細実線で示す。図９をみると、解析窓の影響が現れずに略正確に長期的な上下変動を捉えていることが確認できる。【００９７】
また、本周波数解析手法の予測可能時期を推定するために、１９９１年１月４日から２００５年８月１１日までの約１５年分の実際の日経平均株価終値のデータから、学習期間としての初期２年間のデータと評価のための終期２年間のデータとを除いた約１１年間の区間におけるデータを用いて、予測結果の検証を行った。具体的には、解析期間を一定として１９９３年２月１日から２００３年７月９日まで１日ずつ移動させながら図９と同様に本周波数解析手法を適用して予測データを求め、過去２年間のデータにおける各時点について、当該各時点以降２年間の実際の株価データと予測データとの誤差の分布を求めた。この結果を図１０に示す。【００９８】
図１０をみると、予測誤差が１０％以下の領域が多いことが窺える。なお、１９９７年前後の期間や２００１年前後の期間において誤差が大きくなっているのは、前者期間では、１９９７年７月に発生したアジア通貨危機により、後者期間では、２００１年９月に発生したアメリカ同時多発テロにより、株価に影響を与え、慣性的な変動が維持できなくなったためと考えられる。このような大事件が発生して突発的な変動をした直後においては誤差が大きくなる傾向にあるが、その後は再度誤差が小さくなることも確認できる。これは、解析期間を２年としているため、事件発生からある程度の期間経過までの予測においては事件発生前のデータの周期性を引きずらず、再度周期性が安定し出す頃には予測結果に正確な周期性が反映されるためであると考えられる。このことから、本周波数解析手法においては、２年程度の解析期間のスペクトル構成を正確に推定し、その時点でその後の２年間の変動を長期的に予測した場合であっても、予測しがたい事件が発生しなければ、予測誤差を小さく抑えることができることがわかった。【００９９】
このように、本周波数解析手法は、非常に高い周波数分解能を有することから、一部の区間を解析すれば、その信号の特性を本質的に見抜くことができる可能性がある。従来、この分野においては、いわゆるランダムウォークに代表される現象のモデル化、遺伝的アルゴリズム（Genetic Algorithm；ＧＡ）やニューラルネットワーク、さらにはエージェントアプローチ等に代表される人工知能的なアプローチ等、様々な方法が試みられているが、これほどの長期の予測が行われた事例は未だ報告されていない。この事例からもわかるように、本周波数解析手法は、他の手法を大きく凌駕する高い解析精度を有し、経済指標への応用分野だけでも、株価の中期・長期予測、金融政策における経済指標の影響予測、投資計画、政策決定等、様々なリスク管理に応用可能であるといえる。【０１００】
［本周波数解析手法の効果］
以上説明したように、本周波数解析手法は、非線形方程式の最適解を求めることにより、周波数分解能が解析窓長に依存せず、極めて高い信号解析精度を実現することができる。このような本周波数解析手法は、様々な工学系の分野において基本技術に関わるため、その利用価値が極めて大きい。例えば、本周波数解析手法は、高い周波数分解能を有するため、周波数解析を必要とする既存分野は勿論のこと、ＤＦＴに代表される調和形の周波数解析では限界があった応用分野に適用することにより、解析効果を大幅に改善することができる見通しがある。また、本周波数解析手法は、解析窓に依存しないことから、連続に変動する信号の一部を観測することにより、信号全体の特性を把握することができるため、予測問題への応用等も可能である。【０１０１】
なお、本発明は、上述した実施の形態に限定されるものではない。例えば、上述した実施の形態では、信号解析装置によってソフトウェアによる周波数解析を行うものとして説明したが、本発明は、本周波数解析手法のアルゴリズムを実装したＤＳＰ（Digital Signal Processor）等、積和演算を行うことが可能であればハードウェアによっても実現することができる。【０１０２】
このように、本発明は、その趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【図１】
0005590547.tif 000033
【図２】
0005590547.tif 000034
【図３】
0005590547.tif 000035
【図４】
0005590547.tif 000036
【図５】
0005590547.tif 000037
【図６（Ａ）】
0005590547.tif 000038
【図６（Ｂ）】
0005590547.tif 000039
【図７】
0005590547.tif 000040
【図８】
0005590547.tif 000041
【図９】
0005590547.tif 000042
【図１０】
0005590547.tif 000043