Top > Search of Japanese Patents > (In Japanese)信号解析方法 > Specification

Specification :(In Japanese)信号解析方法

Country (In Japanese)日本国特許庁(JP)
Gazette (In Japanese)特許公報(B2)
Patent Number P5590547
Date of registration Aug 8, 2014
Date of issue Sep 17, 2014
Title of the invention, or title of the device (In Japanese)信号解析方法
IPC (International Patent Classification) G01R  23/16        (2006.01)
G10L  25/18        (2013.01)
FI (File Index) G01R 23/16 D
G10L 25/18
Number of claims or invention 1
Total pages 25
Application Number P2009-533139
Date of filing Sep 17, 2008
International application number PCT/JP2008/066689
International publication number WO2009/038056
Date of international publication Mar 26, 2009
Application number of the priority 2007243858
Priority date Sep 20, 2007
Claim of priority (country) (In Japanese)日本国(JP)
Date of request for substantive examination Sep 13, 2011
Patentee, or owner of utility model right (In Japanese)【識別番号】305060567
【氏名又は名称】国立大学法人富山大学
Inventor, or creator of device (In Japanese)【氏名】広林 茂樹
Representative (In Japanese)【識別番号】100105809、【弁理士】、【氏名又は名称】木森 有平
Examiner (In Japanese)【審査官】吉田 久
Document or reference (In Japanese)国際公開第2006/010002(WO,A2)
特開2007-207387(JP,A)
特開2007-193783(JP,A)
特開2005-275588(JP,A)
特開昭61-239169(JP,A)
特開2000-181504(JP,A)
特開2003-227854(JP,A)
特開2001-116780(JP,A)
広林 茂樹 Shigeki Hirobayashi,柴野 洋平 Yohei Shibano,山淵 龍夫 Tatsuo Yamabuchi,“高分解能の周波数解析法を用いたスペクトルサブトラクションの改善 Improvement of the Spectral Subtraction Method by High-resolution Frequency Analysis”,電気学会論文誌C 電子・情報・システム部門誌,社団法人電気学会,2005年 1月 1日,Vol.125, No.1,p.147-148
東山 三樹夫 Mikio Tohyama,小池 恒彦 Tsunehiko Koike,“高い周波数分析精度の信号分析手法 High resolution frequency analysis”,日本音響学会誌,社団法人日本音響学会,1998年 8月 1日,VOL.54,NO.8,p.568-574
Field of search G01R 23/00~23/20
G10L 25/00~25/93
CiNii
Scope of claims (In Japanese)【請求項1】
入力された非同期信号の周波数解析を行うに際し、サンプリングされた解析対象信号と正弦波モデル信号との差の二乗和が最小値になるような正弦波モデル信号のパラメータである周波数f’と位相φ’と振幅A’についてそれぞれ初期値を求め、求めた各初期値から非線形方程式の最適解に収束させる信号解析方法であって、信号解析装置又はコンピュータを用いて、前記パラメータのうち周波数f’と位相φ’について最急降下法を適用して収束させてから減速法を取り入れたニュートン法を適用して収束させ、その後、前記パラメータのうち振幅A’について最急降下法を適用して収束させてから減速法を取り入れたニュートン法を適用して収束させることとし、
前記解析対象信号が1次元の計測値s(xn)であり、データ個数をNとすると、数式(1)の右辺が最小値になるような周波数fx’と位相φ’と振幅A’を求めることを特徴とし、
前記解析対象信号が2次元の計測値s(xm,yn)であり、各データ個数をM,Nとすると、数式(2)の右辺が最小値になるような周波数fx’,fy’と振幅A’と位相φ’を求めることを特徴とし、
前記解析対象信号が3次元の計測値s(xl,ym,zn)であり、各データ個数をL,M,Nとすると、数式(3)の右辺が最小値になるような周波数fx’,fy’,fz’と振幅A’と位相φ’を求めることを特徴とし、
前記解析対象信号が4次元の計測値s(xk,yl,zm,wn)であり、各データ個数をK,L,M,Nとすると、数式(4)の右辺が最小値になるような周波数fx’,fy’,fz’,fw’と振幅A’と位相φ’を求めることを特徴とし、
前記解析対象信号がn次元(nは5以上の整数)の計測値s(x1n1,x2n2,x3n3,・・・,xnnn)であり、各データ個数をN1,N2,N3,・・・,NNとすると、数式(5)の右辺が最小値になるような周波数fx1’,fx2’,fx3’,・・・,fxn’と振幅A’と位相φ’を求めることを特徴とする信号解析方法。
【数1】
JP0005590547B2_000028t.gif
【数2】
JP0005590547B2_000029t.gif
【数3】
JP0005590547B2_000030t.gif
【数4】
JP0005590547B2_000031t.gif
【数5】
JP0005590547B2_000032t.gif
Detailed description of the invention (In Japanese)【技術分野】
【0001】
本発明は、入力された非同期信号の周波数解析を行う信号解析方法に関する。
【背景技術】
【0002】
従来から、信号の計測、ディジタル信号の圧縮符号化、及び、将来における経済時系列信号変動の予測等、様々な分野において、入力された信号の周波数解析が行われている。いわゆる正弦波モデルを用いた信号解析方法としては、非特許文献1にも紹介されているように、様々な手法が提案されているが、その多くは、信号の周波数を事前に推定し、その推定した周波数における最適な振幅や初期位相を求めることによって行われている。
【0003】
信号は複数の周波数の正弦波の重ね合わせであり、信号を解析するには、当該信号を構成する様々な周波数の正弦波毎に振幅や初期位相等の周波数パラメータを推定することになる。このような周波数解析手法としては、フーリエ変換が広く知られている。フーリエ変換における周波数パラメータの推定は、周期的な信号の場合には、次式(1)に示すように、有限長の窓によるフーリエ変換式を用いて行うことができ、この次式(1)に基づいて、特定の周波数f(=n/T:nは整数、Tは解析窓長)のパラメータ(振幅A及び初期位相φ)を推定することができる。ただし、振幅Aと初期位相φとX(f)との関係は次式(2)である。
【0004】
【数1】
JP0005590547B2_000002t.gif
【数2】
JP0005590547B2_000003t.gif

【0005】
このようなフーリエ変換によって推定することが可能な周波数パラメータは、解析窓長Tに依存し、1/Tの整数倍の周波数分しか計算することができない。したがって、上式(1)を用いてパラメータを推定している限り、周波数fは、原理的に離散的となり、非周期信号の周波数成分を特定することはできない。また、フーリエ変換を用いて周波数解析を行っている分野では、そのほとんどが、周期を持たない非周期信号を解析対象とするにもかかわらず、無限長の積分区間にわたる積分計算が極めて困難であることから、上式(1)を用いなければならないのが現状である。
【0006】
この問題を解決するために、窓関数の調整や解析区間を複数利用して見かけ上の周波数分解能を向上させる試みがある。
【0007】
例えば、非特許文献2に記載の短時間フーリエ変換(Short-Time Fourier Transform、Short-Term Fourier Transform;STFT)は、音響信号を短時間の解析窓によって複数の信号に分割し、各信号に対して独立にフーリエ変換を行う手法である。
【0008】
また、非特許文献3乃至非特許文献6に記載のABS(Analysis by Synthesis)法やGHA(Generalized Harmonic Analysis)等においては、最適な周波数を定めるための周波数を複数用意して振幅及び初期位相を求め、その中から最適な周波数と振幅と初期位相との組み合わせを選択している。
【0009】

【非特許文献1】東山三樹夫、小池恒彦、「高い周波数分析精度の信号分析手法」、日本音響学会論文誌、1988年、第54巻、第8号
【非特許文献2】L. R. Rabiner and R. W. Schafer,“Digital Processing of Speech Signals”, Prentice-Hall, Englewood Cliffs, NJ, 1978年
【非特許文献3】T. F. Quatieriand R. J. McAulay, “Speechtransformations based on a sinesoidal representation”, IEEE Trans. Acoust. Speech Signal Process. ASSP 34, 1986年, p.1449-1464
【非特許文献4】E. B. Georgeand M. J. Smith, “Analysis-by-synthesis/overlap-addsinusoidal representation”, J. Audio Eng. Soc.40, 1992年, p.497-515
【非特許文献5】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
【非特許文献6】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
【0010】
しかしながら、上述した非特許文献2に記載の短時間フーリエ変換においては、短時間の解析窓長の逆数の整数倍の周波数でしか、周波数パラメータを推定することができないという問題がなおも存在する。
【0011】
また、上述した非特許文献3乃至非特許文献6に記載のABS法やGHA等においては、元の信号に対して、予め用意された周波数と異なる周波数が存在する場合には、誤った組み合わせが複数検波されてしまい、解析精度が低下するという問題がある。
【0012】
これら窓関数の調整や解析区間を複数利用して見かけ上の周波数分解能を向上させる試みは、
1)周波数分解能を上げるためには解析窓長を長くしなければならないが、時間的な周波数変動をともなう信号では周波数が平均化され、時間的な変動を捉えることが困難となること
2)一方、解析窓長を短くすると、周波数分解能が低下すること
という問題が指摘されている。
【0013】
このように、周波数の変化を細かく捉えることと、周波数分解能を上げることとは、互いに相反する関係にあり、信号解析上の大きな問題となっていた。
【発明の開示】
【0014】
本発明は、このような実情に鑑みてなされたものであり、周波数分解能が解析窓長に依存せずに高い周波数分解能を有し、驚くべき高精度に周波数解析を行うことができる信号解析方法を提供することを目的とする。
【0015】
本発明の信号解析方法は、入力された非同期信号の周波数解析を行うに際し、サンプリングされた解析対象信号と正弦波モデル信号との差の二乗和が最小値になるような正弦波モデル信号のパラメータである周波数f’と位相φ’と振幅A’についてそれぞれ初期値を求め、求めた各初期値から非線形方程式の最適解に収束させる信号解析方法であって、信号解析装置又はコンピュータを用いて、前記パラメータのうち周波数f’と位相φ’について最急降下法を適用して収束させてから減速法を取り入れたニュートン法を適用して収束させ、その後、前記パラメータのうち振幅A’について最急降下法を適用して収束させてから減速法を取り入れたニュートン法を適用して収束させることとし、
前記解析対象信号が1次元の計測値s(xn)であり、データ個数をNとすると、数式(1)の右辺が最小値になるような周波数fx’と位相φ’と振幅A’を求めることを特徴とし、
前記解析対象信号が2次元の計測値s(xm,yn)であり、各データ個数をM,Nとすると、数式(2)の右辺が最小値になるような周波数fx’,fy’と振幅A’と位相φ’を求めることを特徴とし、
前記解析対象信号が3次元の計測値s(xl,ym,zn)であり、各データ個数をL,M,Nとすると、数式(3)の右辺が最小値になるような周波数fx’,fy’,fz’と振幅A’と位相φ’を求めることを特徴とし、
前記解析対象信号が4次元の計測値s(xk,yl,zm,wn)であり、各データ個数をK,L,M,Nとすると、数式(4)の右辺が最小値になるような周波数fx’,fy’,fz’,fw’と振幅A’と位相φ’を求めることを特徴とし、
前記解析対象信号がn次元(nは5以上の整数)の計測値s(x1n1,x2n2,x3n3,・・・,xnnn)であり、各データ個数をN1,N2,N3,・・・,NNとすると、数式(5)の右辺が最小値になるような周波数fx1’,fx2’,fx3’,・・・,fxn’と振幅A’と位相φ’を求めることを特徴とする信号解析方法。
【数1】
JP0005590547B2_000004t.gif
【数2】
JP0005590547B2_000005t.gif
【数3】
JP0005590547B2_000006t.gif
【数4】
JP0005590547B2_000007t.gif
【数5】
JP0005590547B2_000008t.gif

【0018】
このような本発明にかかる信号解析方法、信号解析装置、及び信号解析プログラムが実装された装置は、非周期信号を想定した周波数解析を目的として定式化を行い、非周期信号のフーリエ変換式のパラメータを求める問題を非線形方程式の最適解を求める問題に置き換えることにより、周波数分解能が解析窓長に依存しない新たな周波数解析手法を実現することができる。
【0019】
以上のような本発明においては、解析窓長と周波数分解能の制約を受けることなく高い周波数分解能を有し、従来の周波数解析手法における精度と比べて10万~100億倍以上という驚くべき高精度に周波数解析を行うことができる。
【図面の簡単な説明】
【0020】
【図1】本発明の実施の形態として示す信号解析装置の構成を示すブロック図である。
【図2】本周波数解析手法とDFTとの違いを説明するための図である。
【図3】本周波数解析手法とDFTとの違いを説明するための図であり、2次元の解析対象信号に各手法を適用した場合の具体例を示す図である。
【図4】本周波数解析手法とDFTとGHAとの違いを説明するための図であり、各手法の誤差を求めた結果を示す図である。
【図5】解析対象信号としてのピアノ音の具体例を示す図である。
【図6(A)】DFTを用いて図5に示す信号を解析して得られたスペクトログラムを示す図である。
【図6(B)】本周波数解析手法を用いて図5に示す信号を解析して得られたスペクトログラムを示す図である。
【図7】本周波数解析手法とDFTとの違いを説明するための図であり、2次元の解析対象信号に各手法を適用して得られるスペクトルの解析精度の具体例を示す図である。
【図8】本周波数解析手法とJPEGとJPEG2000との違いを説明するための図であり、2次元のビットマップ画像データを圧縮した場合の具体例を示す図である。
【図9】本周波数解析手法を適用した経済指標の予測例を示す図である。
【図10】本周波数解析手法を適用して求めた予測データと実際のデータとの誤差の分布を示す図である。

【発明を実施するための最良の形態】
【0021】
以下、本発明を適用した具体的な実施の形態について図面を参照しながら詳細に説明する。
【0022】
この実施の形態は、入力された信号の周波数解析を行う信号解析装置である。特に、この信号解析装置は、非線形方程式を解くことでフーリエ係数を推定することにより、周波数分解能が解析窓長に依存しない新たな周波数解析手法を実現するものである。
【0023】
[信号解析装置の構成]
信号解析装置は、例えばコンピュータ等から構成され、図1に示すように、各部を統括的に制御するCPU(Central Processing Unit)11と、各種プログラムを含む各種情報を格納する読み取り専用のROM(Read Only Memory)12と、ワークエリアとして機能するRAM(Random Access Memory)13と、各種情報を読み出し及び/又は書き込み可能に記憶する記憶部14と、ユーザインターフェースとしての図示しない所定の操作デバイスを介した入力操作の処理及び制御を行う入力操作制御部15と、各種情報を表示する表示部16とを備える。
【0024】
CPU11は、記憶部14等に格納されている各種アプリケーションプログラムをはじめとする各種プログラムを実行し、各部を統括的に制御する。
【0025】
ROM12は、各種プログラムをはじめとする各種情報を格納している。このROM12に格納されている情報は、CPU11の制御のもとに読み出される。
【0026】
RAM13は、CPU11が各種プログラムを実行する際のワークエリアとして機能し、CPU11の制御のもとに、各種情報を一時記憶するとともに、記憶している各種情報を読み出す。
【0027】
記憶部14は、本発明にかかる信号解析プログラム等のアプリケーションプログラムの他、周波数解析の対象となる解析対象信号のデータをはじめとする各種情報を記憶する。この記憶部14としては、例えば、ハードディスクや不揮発性メモリ等を用いることができる。また、記憶部14には、本体に対して着脱可能とされるフレキシブルディスクやメモリカード等の記憶媒体に対して、各種情報の読み出し及び/又は書き込みを行うドライブ装置も含まれる。この記憶部14に記憶されている各種情報は、CPU11の制御のもとに読み出される。
【0028】
入力操作制御部15は、例えば、キーボード、マウス、キーパッド、赤外線リモートコントローラ、スティックキー、又はプッシュボタンといった、ユーザインターフェースとしての図示しない所定の操作デバイスを介した入力操作を受け付け、操作内容を示す制御信号をCPU11に対して供給する。
【0029】
表示部16は、例えば、液晶ディスプレイ(Liquid Crystal Display;LCD)、プラズマ・ディスプレイ・パネル(Plasma Display Panel;PDP)、有機エレクトロルミネッセンス(Organic ElectroLuminescent)ディスプレイ、又はCRT(Cathode Ray Tube)といった、各種表示デバイスであり、CPU11の制御のもとに各種情報を表示する。例えば、表示部16は、CPU11によって信号解析プログラムが起動されると、その画面を表示し、入力された解析対象信号や周波数解析結果等を表示する。
【0030】
このような各部を備える信号解析装置は、CPU11の制御のもとに、信号解析プログラムを実行すると、CPU11の制御のもとに、入力された信号の周波数解析を行う。なお、周波数解析の対象となる信号、すなわち、上述した解析対象信号は、図示しない信号入力部を介して入力され、CPU11によって読み出し可能に記憶部14に記憶される。例えば、信号解析装置は、任意のセンサによって計測された時系列信号を解析対象信号とする場合には、当該信号解析装置とセンサとを接続する所定のインターフェースを介して解析対象信号としての時系列信号を記憶部14に記憶させることによって入力する。また、信号解析装置は、ディジタルカメラ等の撮像装置によって撮像された画像データの圧縮符号化を行う機能を有し、その圧縮符号化の際に周波数解析を行う場合には、当該信号解析装置と撮像装置とを接続する所定のインターフェースを介して解析対象信号としての画像データを記憶部14に記憶させることによって入力する。さらに、信号解析装置は、株価等の将来における変動を予測するために用いる場合には、インターネット等を介して解析対象信号としての過去の株価データを記憶部14に記憶させることによって入力する。その他、信号解析装置は、ユーザが作成した任意のデータを解析対象信号として記憶部14に記憶させることによって入力するようにしてもよい。すなわち、信号入力部は、解析対象信号をCPU11によって読み出し可能に入力させる機能を有する部位である。なお、信号入力部は、アナログ信号を入力した場合には、A/D変換を行ってディジタル信号に変換する機能をあわせ持つことはいうまでもない。このとき、信号入力部は、必要に応じてアンチエイリアシングフィルタを含むA/D変換器としてもよい。信号解析装置は、CPU11の制御のもとに、このようにして入力された解析対象信号の周波数解析を行い、その解析結果等を、図示しない解析結果出力部を介して表示部16に表示出力させたり、印刷装置やその他の機器に出力したりする。
【0031】
[周波数解析アルゴリズム]
以下、このような信号解析装置における周波数解析アルゴリズムについて詳述する。
【0032】
信号解析装置に適用する周波数解析手法、すなわち、本発明にて新たに提案する周波数解析手法(以下、本周波数解析手法という。)は、「周波数の変化を細かく捉えることと、周波数分解能を上げることとは、互いに相反する関係にある」という従来の周波数解析手法の問題を本質的に解決するために、非線形方程式の解を求める問題に着目した。すなわち、本周波数解析手法においては、次式(3)に示す非周期信号のフーリエ変換式の周波数パラメータを求める問題を非線形方程式の最適解を求める問題に置き換えた。
【0033】
【数3】
JP0005590547B2_000009t.gif

【0034】
具体的には、本周波数解析手法においては、次式(4)に示すように、解析対象信号x(n)と正弦波モデル信号との差の二乗和で表される非線形方程式の最適解として、この非線形方程式の右辺が最小値になるような周波数f’、振幅A’、及び初期位相φ’を求める。なお、次式(4)において、Lはフレーム長(解析窓長)であり、fsはサンプリング周波数[Hz]である。このように、非線形方程式に周波数f’をパラメータとして盛り込んだことは、いまだ事例がなく斬新である。換言すれば、本周波数解析手法においては、振幅A’及び初期位相φ’のみならず、周波数f’をも正確に求めることを目的とし、上式(3)に対して非線形方程式の解法を用いるものである。本周波数解析手法においては、このような最小二乗法によって非線形方程式の最適解を求める問題に帰着させることにより、解析窓の影響やエイリアシングの影響がなくなり、解析窓長が、1周期未満であってもよく、周期の整数倍でなくてもよく、さらには、不等間隔であってもよい等、柔軟な周波数解析処理を実現することが可能となる。
【0035】
【数4】
JP0005590547B2_000010t.gif

【0036】
なお、上式(4)は、解析対象信号x(n)が、センサによって計測された時系列信号等の1次元信号である場合を想定しているが、本周波数解析手法は、2次元の画像データ等の2次元信号や動画像データ等の3次元信号、ホログラムのような立体動画像データ等の4次元信号、さらにはそれ以上の次元の信号に対しても適用可能である。すなわち、本周波数解析手法は、1次元の解析対象信号をs(xn)、2次元の解析対象信号をs(xm,yn)、3次元の解析対象信号をs(xl,ym,zn)、4次元の解析対象信号をs(xk,yl,zm,wn)として一般化すると、解析対象信号が1次元、2次元、3次元、又は4次元の場合のそれぞれについて、次式(5)乃至次式(8)に示す非線形方程式の最適解を求めることになる。
【0037】
【数5】
JP0005590547B2_000011t.gif
【数6】
JP0005590547B2_000012t.gif
【数7】
JP0005590547B2_000013t.gif
【数8】
JP0005590547B2_000014t.gif

【0038】
換言すれば、本周波数解析手法は、任意のn次元(nは1以上の整数)の解析対象信号をs(x1n1,x2n2,x3n3,・・・,xnnn)として一般化すると、次式(9)に示す非線形方程式の最適解を求めることになる。
【数9】
JP0005590547B2_000015t.gif

【0039】
なお、以下では、説明の便宜上、解析対象信号が1次元であり、上式(4)に示す非線形方程式の最適解を求めるものとする。
【0040】
さて、上式(4)に示す非線形方程式の最適解を実際に求めるにあたっては、以下のような方法をとることができる。
【0041】
本周波数解析手法においては、振幅A’、周波数f’、及び初期位相φ’のそれぞれについて適切な初期値を求め、これら初期値から非線形方程式の解法を用いて最適解に収束させる。この非線形問題では、上式(4)をコスト関数とする最小化問題とする。なお、適切な初期値は、離散フーリエ変換(Discrete Fourier Transform;以下、DFTという。)やウェーブレット変換等の任意の周波数変換を行ったり、フィルタリングを行うことによっておおよその見当をつけたりする等、既存の任意の方法を適用して求めることができる。
【0042】
まず、本周波数解析手法においては、上式(4)における正弦波モデル信号の位相を構成する周波数パラメータf’,φ’について、いわゆる最急降下法を適用し、周波数パラメータfm’,φm’を次式(10)及び次式(11)によって求める。
【0043】
【数10】
JP0005590547B2_000016t.gif
【数11】
JP0005590547B2_000017t.gif

【0044】
なお、上式(10)及び上式(11)においては、次式(12)と略している。また、μmは、いわゆる減速法に基づく重み係数であり、各漸化式によって求められるコスト関数を単調減少数列にするために、適時0~1の値をとる。
【0045】
【数12】
JP0005590547B2_000018t.gif

【0046】
周波数パラメータfm’,φm’を求めることができれば、上式(4)における正弦波モデル信号の係数としての周波数パラメータA’を一意に求めることができるため、本周波数解析手法においては、次式(13)によって周波数パラメータAm’を収束させる。
【0047】
【数13】
JP0005590547B2_000019t.gif

【0048】
本周波数解析手法においては、これら一連の計算を反復して行うことにより、振幅A’、周波数f’、及び初期位相φ’を高精度に収束させることができる。特に、本周波数解析手法においては、上式(4)における正弦波モデル信号の位相を構成する周波数パラメータf’,φ’と、係数としての周波数パラメータA’とを別個に求めることにより、計算を簡便に行うことができる。
【0049】
しかしながら、最急降下法は、比較的広い範囲から収束するものの、1回の反復では精度が低く、収束するまでに時間を要する。
【0050】
そこで、本周波数解析手法においては、最急降下法を適用して周波数パラメータfm’,φm’をある程度まで収束させた後、さらに、いわゆるニュートン法を適用して高精度に収束させるのが望ましい。具体的には、本周波数解析手法においては、ニュートン法として、次式(14)及び次式(15)に示す漸化式によって周波数パラメータfm’,φm’を求める。
【0051】
【数14】
JP0005590547B2_000020t.gif
【数15】
JP0005590547B2_000021t.gif

【0052】
ただし、上式(14)及び上式(15)において、Jは次式(16)とし、次式(17)と略している。また、νmもμmと同様に減速法に基づく重み係数であり、適時0~1の値をとる。
【0053】
【数16】
JP0005590547B2_000022t.gif
【数17】
JP0005590547B2_000023t.gif

【0054】
本周波数解析手法においては、上式(14)及び上式(15)によって周波数パラメータfm’,φm’を求めた後、最急降下法と同様に、上式(13)によって周波数パラメータAm’を収束させ、この一連の計算をさらに反復して行う。
【0055】
このように、本周波数解析手法においては、最急降下法とニュートン法とを組み合わせたハイブリッド型の解法を用いることにより、高速に且つ高精度に周波数パラメータA’,f’,φ’を推定することができる。
【0056】
また、本周波数解析手法においては、解析対象信号x(n)が複合正弦波の場合であっても、逐次減算処理することにより、近似的にスペクトルパラメータを導出することができる。ここで、解析対象信号x(n)が複数の正弦波の和であり、次式(18)のように表されているとする。
【0057】
【数18】
JP0005590547B2_000024t.gif

【0058】
パーセヴァル(Parseval)の定理より、解析対象信号x(n)の周波数fkと正弦波モデル信号の周波数パラメータf’とが全く一致しない場合、すなわち、次式(19)である場合には、上式(4)に示す非線形方程式は次式(20)となる。また、周波数パラメータf’,φ’の組が、周波数fk及び初期位相φkの組のいずれかに一致する場合には、上式(4)に示す非線形方程式は次式(21)となる。さらに、振幅Ajが周波数パラメータA’とも一致した場合には、解析対象信号から推定スペクトルに関する周波数成分を完全に消去することができる。そのため、最適解を求める問題は、周波数に対して独立であり、解析対象信号から順次個別に推定すれば、複数の正弦波で表される信号にも応用することができる。
【0059】
【数19】
JP0005590547B2_000025t.gif
【数20】
JP0005590547B2_000026t.gif
【数21】
JP0005590547B2_000027t.gif

【0060】
すなわち、本周波数解析手法においては、解析対象信号x(n)が複合正弦波の場合であっても、逐次残差信号に対して同様に処理を行い、複数の正弦波を抽出することができる。
【0061】
ただし、複数のスペクトルの場合には、上式(3)からもわかるように、無限長を仮定しているため、複数のスペクトルの周波数同士が近い場合には、合成スペクトルの周期が解析窓長よりも数段長くなり、上式(17)が満たされなくなるため、誤差が発生する場合もある。
【0062】
音声信号や音響信号等の信号を複合正弦波によって表現するためには、これまで多くのスペクトル数(正弦波の数)が必要であったが、本周波数解析手法においては、そのような信号であっても僅かなスペクトル数で誤差なく表現することができる。すなわち、信号をより少ないスペクトル数で表現可能であることは、情報圧縮の用途に有効であることは勿論のこと、それ以外にも、物理現象に関して波数表現をするようなアプローチに用いることにより、信号の本質的な特性も見抜くことができることを示している。このように、本周波数解析手法は、複合正弦波の解析対象信号に対しても有効であり、様々な分野に応用することが可能である。
【0063】
[本周波数解析手法の有効性]
以下、本周波数解析手法の有効性について具体的に説明する。
【0064】
従来、センサによって計測された信号を、計算機を用いて解析するためには、一般に、例えば図2に示すように、その信号の一部を切り出し、DFTといった代表的な調和型の周波数解析手法によって解析していた。
【0065】
しかしながら、解析対象信号を切り出して(窓掛けして)調和型の解析を行うことは、当該解析対象信号が、切り出した一部の信号が完全周期で繰り返される信号であると擬制して解析することに他ならず、当然ながらその周波数特性も、元のスペクトル特性とは完全に一致しないものとなる。また、この信号切り出しの影響を回避するために、窓関数によって影響を軽減することも試みられているが、切り出した一部の信号を完全周期として解析することには変わりなく、本質的な問題解決には至っていないのが現状である。
【0066】
これに対して、本周波数解析手法は、非周期信号を想定した周波数解析を目的として定式化を行っていることから、解析窓長と周波数分解能の制約を受けることなく、正確に周波数とそのパラメータとを推定することができる。
【0067】
また、本周波数解析手法は、周波数の推定についても、従来の離散的な探索から、動的且つ連続的な探索となるため、格段に精度を向上させることができ、革新的な手法である。特に、本周波数解析手法は、系が安定している場合には、現象の一部を観測して解析することによって正確な特性を推定することができるため、予測問題への応用が期待できる。
【0068】
さらに、多次元化に関しても、従来の周波数解析手法においては、解析対象信号が1次元である場合と同様に窓関数の影響を受ける。具体的には、例えば図3に示すように、2次元の解析対象信号(原画像データ)のうち図3左図中の破線部を切り出してDFTを施して解析したとしても、その切り出した区間を完全周期とした画像が延々と繰り返されるものを解析したことと等価であることから、右上図に示すように、復元した画像データが窓関数の影響を受けたものとなり、本質的な特性を見出すことは困難である。
【0069】
これに対して、上式(6)を適用した本周波数解析手法においては、図3右下図に示すように、窓関数の影響が軽減され、本質的な周波数特性を見出し、切り出した区間の周囲の画像も完全に再現することができる。
【0070】
以上のように、本周波数解析手法は、非線形方程式の最適解を求めることにより、正弦波モデル信号の周波数f’、振幅A’、及び初期位相φ’を高速に且つ高精度に求めることができる。具体的な精度を立証するために、本願発明者は、DFTと、DFTの発展型のうち最も解析精度が高いといわれているGHA(Generalized Harmonic Analysis)とを比較対象として精度の検証を行った。
【0071】
なお、DFTやGHAは、1つの解析窓長に見かけ上複数の窓長を持たせていることから、周波数分解能が解析窓長に依存するが、その分解周波数が有限長であり、解析対象信号の周波数が分解周波数以外の周波数となった場合には解析することができず、解析対象信号が正確に解析できる周波数と異なる場合には、最も近い分解周波数の他に、その周辺に小さなスペクトルの周波数(側帯波成分)が現れ、複数の周波数が出現してしまう。
【0072】
このような現象が本周波数解析手法においても生じるか否かについて、すなわち、本周波数解析手法の周波数分解能を検証するために、解析窓長を1秒(1024サンプル)とした1次元の非常に短い単一正弦波を解析し、各手法によって正弦波を1本抽出して元の信号との二乗誤差を調べた。その結果を図4に示す。
【0073】
図4に示すように、DFTにおいては、基本周波数の整数倍以外の周波数における解析精度の悪化がみられた。また、GHAにおいては、1Hz以上の周波数ではDFTと比べて2~5桁程度の精度向上がみられた。これに対して、本周波数解析手法においては、1Hz以上の周波数ではDFTと比べて10桁以上、GHAと比べて5桁以上の精度向上がみられた。すなわち、本周波数解析手法は、既存の周波数解析手法と比べて10万~100億倍以上の精度向上がみられた。
【0074】
特に、1Hz以下の周波数を正確に推定することができるということは、解析窓長を超えた長い周期信号であっても解析可能であることを示しており、株価等の様々な変動要因が含まれる信号に対しても、そのスペクトル構造をより正確に推定可能であることを示している。
【0075】
このように、本周波数解析手法は、最も解析精度が高いといわれているGHAと比べても驚くべき高精度に解析を行うことができる。国内外の様々な研究事例をみても、これほどの精度で信号を解析できる手法は存在せず、本周波数解析手法は、今後、正確な解析が必要とされる様々な分野への応用が期待できる手法であるといえる。
【0076】
[本周波数解析手法の具体的な応用例]
現在、様々な分野において周波数解析が利用されているが、本周波数解析手法は、これらの分野以外へも波及効果が期待できると考えられる。本周波数解析手法は、フーリエ係数パラメータを、フーリエ変換を用いずに高精度に求めることができることから、例えば周期性(波動性)をともなう音響・振動以外のこれまで様々な非定常的な変動として考えられていた現象を、周期性をともなう信号としてモデル化できる可能性もある。以下、本願発明者が行った本周波数解析手法の応用例について説明する。
【0077】
[計測分野への応用例]
近年、ピアノ弦がハンマーで叩かれた場合には、そのピアノ弦は、垂直な振動のみならず水平な振動も加わり、楕円軌道上で回転運動を行っていることが、ピアノ弦の近傍にレーザ計測機器を設置して当該ピアノ弦との距離変位を計測することにより、明らかになっている。また、近年では、1本のピアノ弦は、2つの基本周波数を有し、それらの基本周波数で交互に揺動しながら連続的に変化していることが光学的に実証されている。
【0078】
しかしながら、ピアノ音を計測してその音の周波数解析を行うことにより、かかるピアノ弦の基本周波数での小さな揺れを正確に捉えることは困難である。従来の周波数解析は、DFTや高速フーリエ変換(Fast Fourier Transform;FFT)に代表されるスペクトラムアナライザーを用いて行われることが多かったが、一般に、窓関数の影響や周波数分解能の低さに起因して、ピアノ弦の小さな揺れを解析することはできなかった。このように、音のみの解析については、光学的な計測と同等な成功事例は未だ報告されていない。
【0079】
そこで、本願発明者は、本周波数解析手法を用いて、一般的に用いられているCD(Compact Disc)程度の音質の観測音のみから、ピアノ弦のスペクトル特性の推定やそのスペクトル周辺の揺らぎを可視化することを試みた。
【0080】
具体的には、図5に示すように、ピアノ調律作業中における音名A4(階名ラ、440Hz)のピアノ音を、CDと同じサンプリング周波数44.1kHzで録音した信号を解析対象信号とし、この解析対象信号をDFT及び本周波数解析手法のそれぞれによって解析し、そのスペクトログラムを評価した。解析区間は1秒であり、窓長のデータ点数を4096点(0.093ミリ秒)とし、シフト長を窓長の1/16とした。この窓長を用いた場合、DFTの周波数分解幅は約10.76Hzとなる。
【0081】
解析結果を図6(A)及び図6(B)に示す。図6(A)は、DFTによる解析結果であり、基音付近を拡大したものである。また、図6(B)は、本周波数解析手法による解析結果であり、同様に基音付近を拡大したものである。なお、これら図6(A)及び図6(B)において、右側横軸は時間を示し、左側横軸は周波数を示し、縦軸は振幅を示している。
【0082】
図6(A)においては上述した約10.76Hzの周波数分解幅でスペクトルを表示している。この図6(A)をみると、440Hz近傍に現れているスペクトルが最も大きいため、その周辺に正しいスペクトルがあることが窺える。また、430Hz近傍に波打ったようなエネルギの増減を呈しているスペクトルも確認できる。今回の調律時には、打鍵時以降に外力が作用しないため、このようなエネルギの増減は起こり得ないと考えると、440Hz近傍のスペクトル列周辺で周波数が揺れているものと予想できるが、周波数分解幅以下の小さな変動を解析することができず、正確にどのような現象が起きているのかを把握することができない。
【0083】
これに対して、図6(B)をみると、440Hz近傍に現れているスペクトルが連続的に揺らいでいる様子が明白に把握できる。これは、ピアノ弦が複数の楕円軌道上で回転運動を行っていることによる周波数の変化を捉えたものに他ならず、DFTでは捉えることができなかった周波数の揺れが確認できる。また、本周波数解析手法は、DFTと比べて、同じ解析窓長であっても正確に周波数(440Hz)を推定していることがわかる。
【0084】
このように、DFTにおいては、解析窓の影響のために周波数の揺れが存在する現象をスペクトルの単調減少関数で表現することができないが、本周波数解析手法においては、周波数変動をともなった単調減少スペクトルで表現することができる。そのため、本周波数解析手法は、複雑な周波数変動をともなう単調減少スペクトルの解析に有効である。
【0085】
この実験からも明らかなように、本周波数解析手法は、CD程度の音質の音の解析から、光学的計測と同等の成果を得ることができる画期的な手法であるといえる。
【0086】
[画像の符号化]
本周波数解析手法においては、窓の影響が少ないため、スペクトルに側帯波成分が全く出ない。すなわち、本周波数解析手法は、上述したように、複雑な解析対象信号であっても、僅かなスペクトル数で効率的に誤差なく表現することができる。
【0087】
図7に2次元信号をDFTと本周波数解析手法とのそれぞれによって解析した具体例を示す。図7左図に示すように、本来であれば1本のスペクトルで表現される信号をDFTによって変換した場合には、一般に、変換して得られるスペクトルも、図3右下図に示したような解析窓の影響に起因して、図7右上図に示すように、複数の側帯波成分が群れて現れてしまう。これは、画像の圧縮符号化において広く用いられている離散コサイン変換(Discrete Cosine Transform;DCT)等においても同様である。これに対して、本周波数解析手法は、図7右下図に示すように、元の信号のスペクトルを正確に推定することができる。正確に解析できると、画像データを表現するスペクトル数を格段に削減できる可能性があり、画質を落とさずに飛躍的な情報圧縮効果を期待することができる。
【0088】
実際に、画像圧縮の一般的手法として広く知られているJPEG(Joint Photographic Experts Group)とJPEG2000とを比較対象として、図8中央図に示す2次元のビットマップ画像データをソフトウェア圧縮した。なお、JPEGは、DCTを基礎解析に利用している手法であり、JPEG2000は、ウェーブレット変換を基礎解析に利用している手法であることは周知である。
【0089】
各周波数解析手法を適用して図8中央図に示す画像データを圧縮した後、同画像データに写り込んでいる背景部分を切り出した結果を図8左図に示し、同画像データに写り込んでいる人物の衣服部分を切り出した結果を図8右図に示す。すなわち、図8左図に示す結果は、変化が緩い部分画像(低周波数信号)の圧縮結果であり、図8右図に示す結果は、変化が激しい部分画像(高周波数信号)の圧縮結果である。
【0090】
図8左図をみると、5176バイトの原部分画像を圧縮した結果、JPEGにおいては、740バイトとなり、原部分画像のサイズの14.3%のサイズに圧縮された。また、JPEG2000においては、306バイトとなり、原部分画像のサイズの5.91%のサイズに圧縮された。これに対して、本周波数解析手法においては、4バイトとなり、原部分画像のサイズの0.08%のサイズにまで圧縮された。なお、この本周波数解析手法による圧縮後の画像から原部分画像を復元するのに必要なスペクトル数は、たったの1本であった。
【0091】
一方、図8右図をみると、5176バイトの原部分画像を圧縮した結果、JPEGにおいては、820バイトとなり、原部分画像のサイズの15.49%のサイズに圧縮された。また、JPEG2000においては、307バイトとなり、原部分画像のサイズの5.93%のサイズに圧縮された。これに対して、本周波数解析手法においては、20バイトとなり、原部分画像のサイズの0.39%のサイズにまで圧縮された。なお、この本周波数解析手法による圧縮後の画像から原部分画像を復元するのに必要なスペクトル数は、5本であった。
【0092】
このように、本周波数解析手法は、JPEGやJPEG2000等の一般的な手法と比べて、圧縮後の情報量が約1/100~1/10となり、他を遥かに凌ぐ圧縮率を実現することができる。また、本周波数解析手法は、低周波数信号から高周波数信号まで幅広く高精度に信号の特徴を特定することができる手法でもある。特に、本周波数解析手法においては、変化が緩い信号の解析を行う場合には、解析窓を超えた周期についても特定することができるため、背景等の画像を効率よく符号化することができる。したがって、本周波数解析手法は、ディジタル放送に移行しているテレビジョン放送に本周波数解析手法を適用することにより、携帯端末機を対象としたいわゆるワンセグであっても、ハイビジョン放送と同等若しくはそれ以上の画像品質を提供できる可能性がある。このように、次世代の画像/動画像符号化技術におけるコア技術となり得ることは極めて有意である。
【0093】
[経済時系列の予測]
大気圧等の平衡状態を保たせる作用のある状況下で、空間上の点において気圧の変化が生じると、その圧力変化が波動性をともなって周囲に伝搬する。かかる圧力を商品の売買における圧力に置き換え、経済指標も常に平衡状態を保たせる作用が働くと仮定すると、波動現象のように考えることができる。このような考え方に基づくと、株価等に代表される経済時系列等にも周期的な変動要因を仮定することができ、現実に、エリオット波動にともなう経済時系列の周期性に着目した研究も行われている。
【0094】
しかしながら、経済時系列の周期性に関しては、様々な周期性が複合的に組み合わされたモデルによって表されていると考えられている。しかも、実際の市場においては、その各周期性も時間的に変動する複雑なモデルとなる。
【0095】
そこで、本願発明者は、かかる困難な問題に対しても、本周波数解析手法を応用することを試みた。特に、本周波数解析手法は、解析窓への依存性が小さく、高い周波数分解能を有することから、解析対象信号の一部の区間を切り出して解析し、その後、延長描写を行うことにより、容易に慣性的で長期的な予測が可能となる。
【0096】
過去4年間(1997年7月~2001年6月)における実際の日経平均株価終値の推移を、図9中、太実線で示す。この株価データのうち、過去2年間(1997年7月~1999年6月)を学習期間とし、この学習期間のデータを解析対象信号として本周波数解析手法を適用した。そして、学習期間以降の2年間(1999年7月~2001年6月)の株価がどのように推移するかを推定した。具体的には、過去2年間の学習期間の株価データを解析対象信号として本周波数解析手法を適用し、複数の正弦波を決定する周波数と振幅と初期位相とを求め、これらのパラメータが暫定的に変化しないものとして、それ以降2年間のデータを延長描写することにより、予測データを求めた。その結果を、図9中、細実線で示す。図9をみると、解析窓の影響が現れずに略正確に長期的な上下変動を捉えていることが確認できる。
【0097】
また、本周波数解析手法の予測可能時期を推定するために、1991年1月4日から2005年8月11日までの約15年分の実際の日経平均株価終値のデータから、学習期間としての初期2年間のデータと評価のための終期2年間のデータとを除いた約11年間の区間におけるデータを用いて、予測結果の検証を行った。具体的には、解析期間を一定として1993年2月1日から2003年7月9日まで1日ずつ移動させながら図9と同様に本周波数解析手法を適用して予測データを求め、過去2年間のデータにおける各時点について、当該各時点以降2年間の実際の株価データと予測データとの誤差の分布を求めた。この結果を図10に示す。
【0098】
図10をみると、予測誤差が10%以下の領域が多いことが窺える。なお、1997年前後の期間や2001年前後の期間において誤差が大きくなっているのは、前者期間では、1997年7月に発生したアジア通貨危機により、後者期間では、2001年9月に発生したアメリカ同時多発テロにより、株価に影響を与え、慣性的な変動が維持できなくなったためと考えられる。このような大事件が発生して突発的な変動をした直後においては誤差が大きくなる傾向にあるが、その後は再度誤差が小さくなることも確認できる。これは、解析期間を2年としているため、事件発生からある程度の期間経過までの予測においては事件発生前のデータの周期性を引きずらず、再度周期性が安定し出す頃には予測結果に正確な周期性が反映されるためであると考えられる。このことから、本周波数解析手法においては、2年程度の解析期間のスペクトル構成を正確に推定し、その時点でその後の2年間の変動を長期的に予測した場合であっても、予測しがたい事件が発生しなければ、予測誤差を小さく抑えることができることがわかった。
【0099】
このように、本周波数解析手法は、非常に高い周波数分解能を有することから、一部の区間を解析すれば、その信号の特性を本質的に見抜くことができる可能性がある。従来、この分野においては、いわゆるランダムウォークに代表される現象のモデル化、遺伝的アルゴリズム(Genetic Algorithm;GA)やニューラルネットワーク、さらにはエージェントアプローチ等に代表される人工知能的なアプローチ等、様々な方法が試みられているが、これほどの長期の予測が行われた事例は未だ報告されていない。この事例からもわかるように、本周波数解析手法は、他の手法を大きく凌駕する高い解析精度を有し、経済指標への応用分野だけでも、株価の中期・長期予測、金融政策における経済指標の影響予測、投資計画、政策決定等、様々なリスク管理に応用可能であるといえる。
【0100】
[本周波数解析手法の効果]
以上説明したように、本周波数解析手法は、非線形方程式の最適解を求めることにより、周波数分解能が解析窓長に依存せず、極めて高い信号解析精度を実現することができる。このような本周波数解析手法は、様々な工学系の分野において基本技術に関わるため、その利用価値が極めて大きい。例えば、本周波数解析手法は、高い周波数分解能を有するため、周波数解析を必要とする既存分野は勿論のこと、DFTに代表される調和形の周波数解析では限界があった応用分野に適用することにより、解析効果を大幅に改善することができる見通しがある。また、本周波数解析手法は、解析窓に依存しないことから、連続に変動する信号の一部を観測することにより、信号全体の特性を把握することができるため、予測問題への応用等も可能である。
【0101】
なお、本発明は、上述した実施の形態に限定されるものではない。例えば、上述した実施の形態では、信号解析装置によってソフトウェアによる周波数解析を行うものとして説明したが、本発明は、本周波数解析手法のアルゴリズムを実装したDSP(Digital Signal Processor)等、積和演算を行うことが可能であればハードウェアによっても実現することができる。
【0102】
このように、本発明は、その趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
Drawing
(In Japanese)【図1】
0
(In Japanese)【図2】
1
(In Japanese)【図3】
2
(In Japanese)【図4】
3
(In Japanese)【図5】
4
(In Japanese)【図6(A)】
5
(In Japanese)【図6(B)】
6
(In Japanese)【図7】
7
(In Japanese)【図8】
8
(In Japanese)【図9】
9
(In Japanese)【図10】
10