JP2714880B2 - 物理的な波形を非調和周波数成分に分析する方法 - Google Patents
物理的な波形を非調和周波数成分に分析する方法Info
- Publication number
- JP2714880B2 JP2714880B2 JP51787194A JP51787194A JP2714880B2 JP 2714880 B2 JP2714880 B2 JP 2714880B2 JP 51787194 A JP51787194 A JP 51787194A JP 51787194 A JP51787194 A JP 51787194A JP 2714880 B2 JP2714880 B2 JP 2714880B2
- Authority
- JP
- Japan
- Prior art keywords
- waveform
- residual
- sine
- waveform data
- cosine
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Description
【発明の詳細な説明】 技術分野 本発明は、コンピュータを用いた波形解析、波形合
成、雑音抑圧、信号検出、伝送帯域圧縮の分野に関す
る。
成、雑音抑圧、信号検出、伝送帯域圧縮の分野に関す
る。
背景技術 FFT（高速フーリエ変換）としてよく知られているフ
ーリエ解析は調和解析の一つであり、分析される周波数
は与えられた波形の長さをＴとすると、n/T（ｎ＝1,2,
…）なる高調波の関係になる。ピリオドグラムは周期解
析法の一つとして知られているが、与えられる波形の長
さは分析される周期よりも十分長いことが必要であり、
また周期波形が比較的顕著なものでないと分析の精度が
得られない。相関法の場合も同様である。時系列の分析
で用いられる自己回帰モデル、自己回帰移動平均モデル
（ARMA）は、デジタルフィルタにおける低域フィルタあ
るいは帯域フィルタに相当するものであり、与えられる
波形の長さは上記と同様に十分なものであることが求め
られる。また、非調和解析の一つであるプロニー法は雑
音があると正しい分析が行なえない。
ーリエ解析は調和解析の一つであり、分析される周波数
は与えられた波形の長さをＴとすると、n/T（ｎ＝1,2,
…）なる高調波の関係になる。ピリオドグラムは周期解
析法の一つとして知られているが、与えられる波形の長
さは分析される周期よりも十分長いことが必要であり、
また周期波形が比較的顕著なものでないと分析の精度が
得られない。相関法の場合も同様である。時系列の分析
で用いられる自己回帰モデル、自己回帰移動平均モデル
（ARMA）は、デジタルフィルタにおける低域フィルタあ
るいは帯域フィルタに相当するものであり、与えられる
波形の長さは上記と同様に十分なものであることが求め
られる。また、非調和解析の一つであるプロニー法は雑
音があると正しい分析が行なえない。
本発明の目的は、比較的短い波形から、必ずしも一つ
の周期とは調和しない正弦波の周期と振幅を高い精度で
検出する方法およびその合成処理方法を提供することに
ある。
の周期とは調和しない正弦波の周期と振幅を高い精度で
検出する方法およびその合成処理方法を提供することに
ある。
発明の開示 上記課題を解決するために、本発明は、物理的な波形
を非調和周波数成分に分析する方法であって、（ａ）所
与の区間の物理的波形を対応するデジタル波形データに
変換し、（ｂ）デジタル波形データに第１の所定周期を
有する第１の正弦波形と第１の余弦波形をそれぞれ乗算
して第１及び第２の積を求め、（ｃ）前記第１の所定周
期の整数倍であり前記所与の区間に等しいかそれよりも
小さい第１の所定区間に亘って第１及び第２の積をそれ
ぞれ加算して第１及び第２の和値を求め、（ｄ）前記第
１及び第２の和値に基づきそれぞれ前記第１の正弦波形
及び第１の余弦波形の振幅を求め、（ｅ）このようにし
て求めたそれぞれの振幅と前記第１の所定周期を有する
正弦波形及び余弦波形を前記デジタル波形データから差
し引いて残差波形データを求め、（ｆ）それぞれ第２乃
至第（ｎ）の所定周期を有する第２乃至第（ｎ）の正弦
波形及び余弦波形を用いて前記デジタル波形データに対
し前記ステップ（ｂ）乃至（ｅ）を順次実行することに
より（ｎ−１）個の残差波形データを求め、（ｇ）前記
（ｎ）個の残差波形データのパワーを調べることにより
前記パワーが最小となる前記正弦波形及び余弦波形を第
１の非調和周波数成分とし、（ｈ）最小のパワーの残差
波形データを与える前記正弦波形及び余弦波形を前記デ
ジタル波形データから差し引いて第１の残差波形を求
め、（ｉ）前記第１の残差波形に対し前記ステップ
（ｂ）乃至（ｈ）を順次実行することにより第２の非調
和周波数成分を求め、（ｊ）前記第１の残差波形から前
記第２の非調和周波数成分を差し引いて第２の残差波形
を求め、（ｋ）前記第２の残差波形に対し前記ステップ
（ｂ）乃至（ｈ）を順次実行することにより第３の非調
和周波数成分を求める、という手順を繰り返すことによ
り、第（ｎ）の残差波形から第（ｎ＋１）の非調和周波
数成分を求めることよりなる物理的な波形を非調和周波
数成分に分析する方法を提供する。
を非調和周波数成分に分析する方法であって、（ａ）所
与の区間の物理的波形を対応するデジタル波形データに
変換し、（ｂ）デジタル波形データに第１の所定周期を
有する第１の正弦波形と第１の余弦波形をそれぞれ乗算
して第１及び第２の積を求め、（ｃ）前記第１の所定周
期の整数倍であり前記所与の区間に等しいかそれよりも
小さい第１の所定区間に亘って第１及び第２の積をそれ
ぞれ加算して第１及び第２の和値を求め、（ｄ）前記第
１及び第２の和値に基づきそれぞれ前記第１の正弦波形
及び第１の余弦波形の振幅を求め、（ｅ）このようにし
て求めたそれぞれの振幅と前記第１の所定周期を有する
正弦波形及び余弦波形を前記デジタル波形データから差
し引いて残差波形データを求め、（ｆ）それぞれ第２乃
至第（ｎ）の所定周期を有する第２乃至第（ｎ）の正弦
波形及び余弦波形を用いて前記デジタル波形データに対
し前記ステップ（ｂ）乃至（ｅ）を順次実行することに
より（ｎ−１）個の残差波形データを求め、（ｇ）前記
（ｎ）個の残差波形データのパワーを調べることにより
前記パワーが最小となる前記正弦波形及び余弦波形を第
１の非調和周波数成分とし、（ｈ）最小のパワーの残差
波形データを与える前記正弦波形及び余弦波形を前記デ
ジタル波形データから差し引いて第１の残差波形を求
め、（ｉ）前記第１の残差波形に対し前記ステップ
（ｂ）乃至（ｈ）を順次実行することにより第２の非調
和周波数成分を求め、（ｊ）前記第１の残差波形から前
記第２の非調和周波数成分を差し引いて第２の残差波形
を求め、（ｋ）前記第２の残差波形に対し前記ステップ
（ｂ）乃至（ｈ）を順次実行することにより第３の非調
和周波数成分を求める、という手順を繰り返すことによ
り、第（ｎ）の残差波形から第（ｎ＋１）の非調和周波
数成分を求めることよりなる物理的な波形を非調和周波
数成分に分析する方法を提供する。
以下、本発明の作用を数式を用いて説明する。
波形データＷ（ｍ）（１≦ｍ≦Ｍ）に所定周期Ｔ、振幅
１の一定振幅正弦波形Ｓ（ｍ）および一定振幅余弦波形
Ｃ（ｍ）、すなわち Ｓ（ｍ）＝sin（２πm/T） Ｃ（ｍ）＝cos（２πm/T） …（１） をそれぞれ掛けて所定区間、たとえばｍ＝Ｍ−Ｌからｍ
＝Ｍまで積分した値をそれぞれＸ（Ｔ）およびＹ（Ｔ）
とすれば、 で与えられる。ここで、積分は離散的な数値系列で示さ
れるＷ（ｍ）,S（ｍ）,C（ｍ）の積について、和（サム
メーション）の形で表わされている。
１の一定振幅正弦波形Ｓ（ｍ）および一定振幅余弦波形
Ｃ（ｍ）、すなわち Ｓ（ｍ）＝sin（２πm/T） Ｃ（ｍ）＝cos（２πm/T） …（１） をそれぞれ掛けて所定区間、たとえばｍ＝Ｍ−Ｌからｍ
＝Ｍまで積分した値をそれぞれＸ（Ｔ）およびＹ（Ｔ）
とすれば、 で与えられる。ここで、積分は離散的な数値系列で示さ
れるＷ（ｍ）,S（ｍ）,C（ｍ）の積について、和（サム
メーション）の形で表わされている。
所定周期Ｔの正弦波形の振幅Ａ（Ｔ）および余弦波形
の振幅Ｂ（Ｔ）は、それぞれ、 Ａ（Ｔ）＝Ｘ（Ｔ）／（L/2） Ｂ（Ｔ）＝Ｙ（Ｔ）／（L/2） …（３） で与えられ、残差波形データＲ（ｍ）は、 Ｒ（ｍ）＝Ｗ（ｍ）−Ａ（Ｔ）Ｓ（ｍ）−Ｂ（Ｔ）Ｃ
（ｍ） …（４） と表される。これを二乗して一定区間ＪからＭにわたっ
て加算した残差量Ｑ（Ｔ）は、 で与えられる。
の振幅Ｂ（Ｔ）は、それぞれ、 Ａ（Ｔ）＝Ｘ（Ｔ）／（L/2） Ｂ（Ｔ）＝Ｙ（Ｔ）／（L/2） …（３） で与えられ、残差波形データＲ（ｍ）は、 Ｒ（ｍ）＝Ｗ（ｍ）−Ａ（Ｔ）Ｓ（ｍ）−Ｂ（Ｔ）Ｃ
（ｍ） …（４） と表される。これを二乗して一定区間ＪからＭにわたっ
て加算した残差量Ｑ（Ｔ）は、 で与えられる。
ここでＷ（ｍ）を振幅Ｖ、周期（Ｔ＋ｄ）、位相Ｐを
もった正弦波形として、Ｌ＝nT（ｎ＝1,2,3,…）とすれ
ば、Ａ（Ｔ）およびＢ（Ｔ）は、それぞれ、 Ａ（Ｔ）＝Ｖ｛cos（Ｐ）cos（M1）＋sin（Ｐ）sin
（M2）｝ Ｂ（Ｔ）＝Ｖ｛sin（Ｐ）cos（M3）−cos（Ｐ）sin
（M4）｝ …（６） Mj＝２πdEj/T,0＜Ej＜１ （ｊ＝1,2,3,4） となり、ｄ＝０すなわち前記所定周期とＷ（ｍ）の正弦
波形の周期が一致したときＱ（Ｔ）は極小（最小）とな
る。これより第１の正弦波形の振幅A1および余弦波形の
振幅B1は、 A1＝Vcos（Ｐ） B1＝Vsin（Ｐ） …（７） で与えられ、第１の正弦波形と余弦波形の和は、上記の
振幅Ｖと位相Ｐをもった正弦波形と等価になる。これを
数式で表わすならば、次のようになる。
もった正弦波形として、Ｌ＝nT（ｎ＝1,2,3,…）とすれ
ば、Ａ（Ｔ）およびＢ（Ｔ）は、それぞれ、 Ａ（Ｔ）＝Ｖ｛cos（Ｐ）cos（M1）＋sin（Ｐ）sin
（M2）｝ Ｂ（Ｔ）＝Ｖ｛sin（Ｐ）cos（M3）−cos（Ｐ）sin
（M4）｝ …（６） Mj＝２πdEj/T,0＜Ej＜１ （ｊ＝1,2,3,4） となり、ｄ＝０すなわち前記所定周期とＷ（ｍ）の正弦
波形の周期が一致したときＱ（Ｔ）は極小（最小）とな
る。これより第１の正弦波形の振幅A1および余弦波形の
振幅B1は、 A1＝Vcos（Ｐ） B1＝Vsin（Ｐ） …（７） で与えられ、第１の正弦波形と余弦波形の和は、上記の
振幅Ｖと位相Ｐをもった正弦波形と等価になる。これを
数式で表わすならば、次のようになる。
Vsin（２πm/T＋Ｐ） ＝A1sin（２πm/T）＋B1cos（２πm/T） Ｖ＝{A2+B2}1/2 Ｐ＝tan1(B1/A1) …（８） なお、式５において二乗して加算するところを、絶対
値をとって加算することもできる。更に、式８から明ら
かなように、波形データから正弦波形および余弦波形を
除くことと、その合成でなる正弦波形を除くことは全く
同じである。
値をとって加算することもできる。更に、式８から明ら
かなように、波形データから正弦波形および余弦波形を
除くことと、その合成でなる正弦波形を除くことは全く
同じである。
Ｗ（ｍ）が複数の異なる周期の正弦波形から成る場合
は、Ｑ（Ｔ）が最小もしくは極小となるＴを求め、その
ときのＡ（Ｔ）とＢ（Ｔ）を第１の正弦波形の振幅A1と
余弦波形の振幅B1とし、Ｗ（ｍ）からこれら振幅をもっ
た正弦波形と余弦波形を差し引いて第１の残差波形を求
めて、この第１の残差波形についても同様にして分析を
行ない、新たな正弦波形および余弦波形を求めるという
手順をくり返す。
は、Ｑ（Ｔ）が最小もしくは極小となるＴを求め、その
ときのＡ（Ｔ）とＢ（Ｔ）を第１の正弦波形の振幅A1と
余弦波形の振幅B1とし、Ｗ（ｍ）からこれら振幅をもっ
た正弦波形と余弦波形を差し引いて第１の残差波形を求
めて、この第１の残差波形についても同様にして分析を
行ない、新たな正弦波形および余弦波形を求めるという
手順をくり返す。
分析が十分になされたかどうかを判断するには、残差
波形の収束を次のようにして調べることができる。
波形の収束を次のようにして調べることができる。
Ｗ（ｍ）から第１ないし第ｎの正弦波形を除いた第ｎ
の残差波形と、更に第（ｎ＋１）の正弦波形を除いた第
（ｎ＋１）の残差波形について、その残差量を比較して
収束を見る。第ｎおよび第（ｎ＋１）の残差波形の二乗
値をそれぞれ一定区間ｍ＝Ｊからｍ＝Ｍまで積分（加
算）した値をＩ（ｎ）およびＩ（ｎ＋１）とする。一般
には、ｎ＝1,2,3,…となるに従ってＩ（ｎ）は小さくな
るが、 Ｉ（ｎ＋１）＝Ｉ（ｎ） …（９） となったとき、あるいは近似的に式９が成り立つとき、
そのｎが収束点になる。
の残差波形と、更に第（ｎ＋１）の正弦波形を除いた第
（ｎ＋１）の残差波形について、その残差量を比較して
収束を見る。第ｎおよび第（ｎ＋１）の残差波形の二乗
値をそれぞれ一定区間ｍ＝Ｊからｍ＝Ｍまで積分（加
算）した値をＩ（ｎ）およびＩ（ｎ＋１）とする。一般
には、ｎ＝1,2,3,…となるに従ってＩ（ｎ）は小さくな
るが、 Ｉ（ｎ＋１）＝Ｉ（ｎ） …（９） となったとき、あるいは近似的に式９が成り立つとき、
そのｎが収束点になる。
上記のようにしてＷ（ｍ）の正弦波形成分および余弦
波形成分として振幅、周期がそれぞれA1,A2,…,AkとB1,
B2,…,BkとT1,T2,…,Tkなる値で求められたとすれば、
周期Tnなる振動の振幅Vnは、 Vn＝(An 2+Bn 2)1/2 …（10） （ｎ＝1,2,…,k） で与えられる。なお、周波数は1/Tnで与えられ、位相Pn
は Pn＝tan1(Bn/An) …（11） で与えられる。
波形成分として振幅、周期がそれぞれA1,A2,…,AkとB1,
B2,…,BkとT1,T2,…,Tkなる値で求められたとすれば、
周期Tnなる振動の振幅Vnは、 Vn＝(An 2+Bn 2)1/2 …（10） （ｎ＝1,2,…,k） で与えられる。なお、周波数は1/Tnで与えられ、位相Pn
は Pn＝tan1(Bn/An) …（11） で与えられる。
図１は、上記の分析手順を示すブロック図である。図
１において、D1は波形記憶部、D2,D3,D4およびD7は演算
部、D5は計算値記憶部、D6は比較判別および演算部、D8
はパラメータ記憶部を表わす。入力端子ａから被測定波
形がD1に送られ、D2で周期Ｔなる正弦波形および余弦波
形の振幅Ａ（Ｔ）,B（Ｔ）と、それぞれの振幅をもった
正弦波形および余弦波形が求められ、Ａ（Ｔ）,B（Ｔ）
およびＴはD5に送られ、波形はD3へ送られて被測定波形
から除かれ、残差波形が求められ、D4で残差波形のパワ
ーに相当する残差量が求められ、D5に送られる。Ｔを変
更して、上記D2からの手順が繰り返される。こうして与
えられた総てのＴについての振幅と残差量がD5に記憶さ
れ、残差量の最小となるものがD6で判別され、その周期
T1と一対の振幅A1,B1がD8に記憶され、これに応ずると
ころの第１の正弦波形と余弦波形がD7に送られ、D1に記
憶された波形からこれらが除去され、第１の残差波形が
求められる。D1に記憶されていた波形は第１の残差波形
に更新され、上記被測定波形についてなされたと同様な
処理が行われ、第２の正弦波形と余弦波形を特定するパ
ラメータである周期T2と一対の振幅A2,B2がD8に記憶さ
れ、第２の残差波形がD7で求められる。D1の波形は第２
の残差波形に更新され、以下同様にして与えられたＴの
範囲で、指定された数、例えばＮ回の繰り返し処理によ
って、Ｎ組のパラメータ（Tk,Ak,Bk）が求められ、これ
らはD8に記憶される。D8に記憶されたパラメータを用い
れば、被測定波形を合成することができ、更に合成され
た波形を外挿して予測波形を求めることができる。ま
た、周期（あるいは周波数）をある範囲に制限して上記
合成を行なえば、理想的な帯域通過波形を得ることもで
きる。
１において、D1は波形記憶部、D2,D3,D4およびD7は演算
部、D5は計算値記憶部、D6は比較判別および演算部、D8
はパラメータ記憶部を表わす。入力端子ａから被測定波
形がD1に送られ、D2で周期Ｔなる正弦波形および余弦波
形の振幅Ａ（Ｔ）,B（Ｔ）と、それぞれの振幅をもった
正弦波形および余弦波形が求められ、Ａ（Ｔ）,B（Ｔ）
およびＴはD5に送られ、波形はD3へ送られて被測定波形
から除かれ、残差波形が求められ、D4で残差波形のパワ
ーに相当する残差量が求められ、D5に送られる。Ｔを変
更して、上記D2からの手順が繰り返される。こうして与
えられた総てのＴについての振幅と残差量がD5に記憶さ
れ、残差量の最小となるものがD6で判別され、その周期
T1と一対の振幅A1,B1がD8に記憶され、これに応ずると
ころの第１の正弦波形と余弦波形がD7に送られ、D1に記
憶された波形からこれらが除去され、第１の残差波形が
求められる。D1に記憶されていた波形は第１の残差波形
に更新され、上記被測定波形についてなされたと同様な
処理が行われ、第２の正弦波形と余弦波形を特定するパ
ラメータである周期T2と一対の振幅A2,B2がD8に記憶さ
れ、第２の残差波形がD7で求められる。D1の波形は第２
の残差波形に更新され、以下同様にして与えられたＴの
範囲で、指定された数、例えばＮ回の繰り返し処理によ
って、Ｎ組のパラメータ（Tk,Ak,Bk）が求められ、これ
らはD8に記憶される。D8に記憶されたパラメータを用い
れば、被測定波形を合成することができ、更に合成され
た波形を外挿して予測波形を求めることができる。ま
た、周期（あるいは周波数）をある範囲に制限して上記
合成を行なえば、理想的な帯域通過波形を得ることもで
きる。
図面の簡単な説明 図１は、波形データの非調和分析および合成処理方法
の分析手順を示すブロック図であり、ａは入力端子、D1
は波形記憶部、D2,D3,D4,D7はそれぞれ演算部、D5は計
算値記憶部、D6は比較判別および演算部、D8はパラメー
タ記憶部を示す。
の分析手順を示すブロック図であり、ａは入力端子、D1
は波形記憶部、D2,D3,D4,D7はそれぞれ演算部、D5は計
算値記憶部、D6は比較判別および演算部、D8はパラメー
タ記憶部を示す。
図２は、波形データと本方法により求まるスペクトル
であり、（１）は波形データ、（２）はスペクトルを示
す。
であり、（１）は波形データ、（２）はスペクトルを示
す。
図３は、図２の波形データについてフーリエ解析によ
り求まるスペクトルである。
り求まるスペクトルである。
図４は、波形データ（実線）と予測波形（点線）を示
す図であり、（１）は本方式による場合、（２）はフー
リエ解析による場合を示す。
す図であり、（１）は本方式による場合、（２）はフー
リエ解析による場合を示す。
図５は、波形データの分析、合成、予測の実験例を示
す図であり、（ａ）は被測定波形、（ｂ）は合成波形と
予測波形、（ｃ）は残差波形、（ｄ）は被測定波形のス
ペクトル、（ｃ）は本方法で求まったスペクトル、
（ｆ）はFFTで求まったスペクトルを示す図である。
す図であり、（ａ）は被測定波形、（ｂ）は合成波形と
予測波形、（ｃ）は残差波形、（ｄ）は被測定波形のス
ペクトル、（ｃ）は本方法で求まったスペクトル、
（ｆ）はFFTで求まったスペクトルを示す図である。
図６は、波形データにランダム雑音が加わった場合の
分析、合成、予測の実験例を示す図であり、（ａ）は被
測定波形、（ｂ）は合成波形と予測波形、（ｃ）は残差
波形、（ｄ）は信号波形成分のスペクトル、（ｅ）は本
方法で求まったスペクトル、（ｆ）はFFTで求まったス
ペクトルをそれぞれ示す図である。
分析、合成、予測の実験例を示す図であり、（ａ）は被
測定波形、（ｂ）は合成波形と予測波形、（ｃ）は残差
波形、（ｄ）は信号波形成分のスペクトル、（ｅ）は本
方法で求まったスペクトル、（ｆ）はFFTで求まったス
ペクトルをそれぞれ示す図である。
発明を実施するための最良の形態 500個の数値系列からなる波形データＷ（ｍ）の非調
和分析の実施例を説明する。
和分析の実施例を説明する。
実施例（１） 所定区間を所定周期の１周期分とした場
合に関し、Ｔ＝500におけるＸ（500）はＷ（１）Ｓ
（１）,W（２）Ｓ（２），…,W（500）Ｓ（500）の和で
与えられ、Ｙ（500）はＷ（１）Ｃ（１）,W（２）Ｃ
（２），…,W（500）Ｃ（500）の和で与えられる。この
ときの振幅は式３でＬ＝500として与えられる。同様に
して、Ｔ＝499に対してＸ（499）がｍ＝２から500まで
のＷ（ｍ）Ｓ（ｍ）の積和で、Ｙ（499）がＷ（ｍ）Ｃ
（ｍ）の積和で与えられ、式３でＬ＝499として振幅が
与えられる。以下同様である。ただし各周期Ｔにおける
Ｓ（ｍ）,C（ｍ）は式１においてそのＴの値を代入して
求める。
合に関し、Ｔ＝500におけるＸ（500）はＷ（１）Ｓ
（１）,W（２）Ｓ（２），…,W（500）Ｓ（500）の和で
与えられ、Ｙ（500）はＷ（１）Ｃ（１）,W（２）Ｃ
（２），…,W（500）Ｃ（500）の和で与えられる。この
ときの振幅は式３でＬ＝500として与えられる。同様に
して、Ｔ＝499に対してＸ（499）がｍ＝２から500まで
のＷ（ｍ）Ｓ（ｍ）の積和で、Ｙ（499）がＷ（ｍ）Ｃ
（ｍ）の積和で与えられ、式３でＬ＝499として振幅が
与えられる。以下同様である。ただし各周期Ｔにおける
Ｓ（ｍ）,C（ｍ）は式１においてそのＴの値を代入して
求める。
Ｑ（Ｔ）が極小となる周期を求める場合は、ある周期
における値がその前後の周期の値よりも小さくなるもの
を選択することによっても求められるが、次のようにす
ることもできる。
における値がその前後の周期の値よりも小さくなるもの
を選択することによっても求められるが、次のようにす
ることもできる。
所定周期Ｔを500から１までとした場合これらを分割
して例えば第１の周期帯域をＴ＝500から400、第２の周
期帯域をＴ＝401から300、第３の周期帯域をＴ＝301か
ら200、第４の周期帯域をＴ＝201から100、第５の周期
帯域をＴ＝101から１として、各周期帯域毎にＱ（Ｔ）
が最小となる周期を求め、その中から各帯域において帯
域の端に当たるものを除けばよい。一定区間は例えばｍ
＝１から500までとする。
して例えば第１の周期帯域をＴ＝500から400、第２の周
期帯域をＴ＝401から300、第３の周期帯域をＴ＝301か
ら200、第４の周期帯域をＴ＝201から100、第５の周期
帯域をＴ＝101から１として、各周期帯域毎にＱ（Ｔ）
が最小となる周期を求め、その中から各帯域において帯
域の端に当たるものを除けばよい。一定区間は例えばｍ
＝１から500までとする。
実施例（２） 所定区間を所定周期の２周期分とした場
合に関し、Ｔ＝250におけるＸ（250）はＷ（１）Ｓ
（１）,W（２）Ｓ（２），…,W（500）Ｓ（500）の和で
与えられ、Ｙ（250）はＷ（１）Ｃ（１）,W（２）Ｃ
（２），…,W（500）Ｃ（500）の和で与えられる。振幅
は式３でＬ＝500として与えられる。同様にしてＴ＝249
に対してＸ（249）がｍ＝３から500までのＷ（ｍ）Ｓ
（ｍ）の積和で、Ｙ（249）が同じくＷ（ｍ）Ｃ（ｍ）
の積和で与えられ、振幅は式３でＬ＝498として与えら
れる。以下同様である。ただし各周期ＴにおけるＳ
（ｍ）,C（ｍ）は式１においてそのＴの値を代入して求
める。
合に関し、Ｔ＝250におけるＸ（250）はＷ（１）Ｓ
（１）,W（２）Ｓ（２），…,W（500）Ｓ（500）の和で
与えられ、Ｙ（250）はＷ（１）Ｃ（１）,W（２）Ｃ
（２），…,W（500）Ｃ（500）の和で与えられる。振幅
は式３でＬ＝500として与えられる。同様にしてＴ＝249
に対してＸ（249）がｍ＝３から500までのＷ（ｍ）Ｓ
（ｍ）の積和で、Ｙ（249）が同じくＷ（ｍ）Ｃ（ｍ）
の積和で与えられ、振幅は式３でＬ＝498として与えら
れる。以下同様である。ただし各周期ＴにおけるＳ
（ｍ）,C（ｍ）は式１においてそのＴの値を代入して求
める。
Ｑ（Ｔ）の極小となる周期を求める場合は上記（１）
で説明された方法を用いることができる。
で説明された方法を用いることができる。
実施例（３） 所定区間を所定周期の３周期分とした場
合に関し、Ｔ＝166におけるＸ（166）はＷ（３）Ｓ
（３），…,W（500）Ｓ（500）の和で与えられＹ（16
6）はＷ（３）Ｃ（３），…,W（500）Ｃ（500）の和で
与えられる。振幅は式３でＬ＝498として与えられる。
同様にしてＴ＝165に対してＸ（165）がｍ＝６から500
までのＷ（ｍ）Ｓ（ｍ）の積和で、Ｙ（165）が同じく
Ｗ（ｍ）Ｃ（ｍ）の積和で与えられ、振幅は式３でＬ＝
495として与えられる。以下同様である。ただし各周期
ＴにおけるＳ（ｍ）,C（ｍ）は式１においてそのＴの値
を代入して求める。
合に関し、Ｔ＝166におけるＸ（166）はＷ（３）Ｓ
（３），…,W（500）Ｓ（500）の和で与えられＹ（16
6）はＷ（３）Ｃ（３），…,W（500）Ｃ（500）の和で
与えられる。振幅は式３でＬ＝498として与えられる。
同様にしてＴ＝165に対してＸ（165）がｍ＝６から500
までのＷ（ｍ）Ｓ（ｍ）の積和で、Ｙ（165）が同じく
Ｗ（ｍ）Ｃ（ｍ）の積和で与えられ、振幅は式３でＬ＝
495として与えられる。以下同様である。ただし各周期
ＴにおけるＳ（ｍ）,C（ｍ）は式１においてそのＴの値
を代入して求める。
Ｑ（Ｔ）の極小となる周期を求める場合は上記（１）
で説明された方法を用いることができる。
で説明された方法を用いることができる。
Ｑ（Ｔ）が極小となる正弦波形を除いて残差波形を求
めるごとに、収束状態を判断するために残差波形の二乗
値を所定区間にわたって積分（加算）し、式９を用いて
収束点を求めるが、前述のように帯域を分割して行なう
場合、実施例（１）について述べるならば、積分（加
算）の区間は第１の周期帯域ではｍ＝１から500まで、
第２の周期帯域ではｍ＝100から500まで、第３の周期帯
域ではｍ＝200から500まで、第４の周期帯域ではｍ＝30
0から500まで、第５の周期帯域ではｍ＝400から500まで
となる。実施例（２）について述べるならば、第１から
第３の周期帯域は実施例（１）と同じで、第４の周期帯
域ではｍ＝100から500まで、第５の周期帯域ではｍ＝30
0から500までとなる。同様に実施例（３）の場合、第５
の周期帯域ではｍ＝200から500までとなり、他の周期帯
域は実施例（１）に準ずる。
めるごとに、収束状態を判断するために残差波形の二乗
値を所定区間にわたって積分（加算）し、式９を用いて
収束点を求めるが、前述のように帯域を分割して行なう
場合、実施例（１）について述べるならば、積分（加
算）の区間は第１の周期帯域ではｍ＝１から500まで、
第２の周期帯域ではｍ＝100から500まで、第３の周期帯
域ではｍ＝200から500まで、第４の周期帯域ではｍ＝30
0から500まで、第５の周期帯域ではｍ＝400から500まで
となる。実施例（２）について述べるならば、第１から
第３の周期帯域は実施例（１）と同じで、第４の周期帯
域ではｍ＝100から500まで、第５の周期帯域ではｍ＝30
0から500までとなる。同様に実施例（３）の場合、第５
の周期帯域ではｍ＝200から500までとなり、他の周期帯
域は実施例（１）に準ずる。
図２は、周期Ｔの正弦波の１周期半が波形データであ
る場合の波形（１）と本方法によって分析して求めた場
合のスペクトル（２）を示したものである。
る場合の波形（１）と本方法によって分析して求めた場
合のスペクトル（２）を示したものである。
図３は、図２の波形（１）を調和解析（フーリエ解
析）によって分析した場合に得られるスペクトルを示し
たものである。
析）によって分析した場合に得られるスペクトルを示し
たものである。
図４は、それぞれ実線で示された波形データを分析
し、合成波形を求めて外挿して得た予測波形（点線）を
示したものであり、（１）は本方法による場合、（２）
が調和解析（フーリエ解析）による場合を示している。
し、合成波形を求めて外挿して得た予測波形（点線）を
示したものであり、（１）は本方法による場合、（２）
が調和解析（フーリエ解析）による場合を示している。
本発明による波形データの非調和分析および合成処理
方法によれば、必ずしも一つの振動とは調和しない独立
な振動波形の和で成るような変動波形、典型的なものと
して概周期波形について、それを構成する正弦波形を特
定することができる。特に波形データ長が短い場合で
も、１周期以上の長さがあれば部分的波形からでもスプ
リアスを生じることなくそれを検出することができる。
従って、この分析を物理的変動波形データに適用すれ
ば、その変動の中にある実際の周期的変動成分を検出で
き、経済分野での時系列データに適用すれば従来の方法
よりも精度の高い分析を行なうことができ、更に合成波
形から予測データを得ることもできる。また、信号伝送
（記録）においては、本方法によって求まったパラメー
タ（Ak,Bk,Tk）を伝送（記録）することによって、帯域
圧縮を行なうことができる。この場合、信号波形を所定
のデータ長ごとに区切って分析を行ない、各区間ごとに
上記パラメータを求めて、これを逐次伝送（記録）す
る。受信側（再生側）においては、上記パラメータを用
いて各正弦波を得、これらを加算して信号波形を復元す
る。
方法によれば、必ずしも一つの振動とは調和しない独立
な振動波形の和で成るような変動波形、典型的なものと
して概周期波形について、それを構成する正弦波形を特
定することができる。特に波形データ長が短い場合で
も、１周期以上の長さがあれば部分的波形からでもスプ
リアスを生じることなくそれを検出することができる。
従って、この分析を物理的変動波形データに適用すれ
ば、その変動の中にある実際の周期的変動成分を検出で
き、経済分野での時系列データに適用すれば従来の方法
よりも精度の高い分析を行なうことができ、更に合成波
形から予測データを得ることもできる。また、信号伝送
（記録）においては、本方法によって求まったパラメー
タ（Ak,Bk,Tk）を伝送（記録）することによって、帯域
圧縮を行なうことができる。この場合、信号波形を所定
のデータ長ごとに区切って分析を行ない、各区間ごとに
上記パラメータを求めて、これを逐次伝送（記録）す
る。受信側（再生側）においては、上記パラメータを用
いて各正弦波を得、これらを加算して信号波形を復元す
る。
なお、本方法により、信号波形とランダム雑音が分離
できることを、その分析、合成、予測の実験例によって
示す。
できることを、その分析、合成、予測の実験例によって
示す。
図５は、本分析法による波形のスペクトル推定と、合
成波形による波形予測の実施例を示したものである。同
図において（ａ）は被測定波形、（ｄ）はそのスペクト
ルを示している。被測定波形の区間〔0,512〕を分析し
て推定したスペクトルが（ｅ）、分析で得た正弦波の合
成と、512以後を予測した波形が（ｂ）、区間〔0,512〕
における被測定波形と合成波形の差が（ｃ）である。ス
ペクトルの縦軸は相対振幅、横軸は相対的な周波数1/T
を示す。また、比較のためにFFTによるスペクトル推定
結果を（ｆ）に示す。
成波形による波形予測の実施例を示したものである。同
図において（ａ）は被測定波形、（ｄ）はそのスペクト
ルを示している。被測定波形の区間〔0,512〕を分析し
て推定したスペクトルが（ｅ）、分析で得た正弦波の合
成と、512以後を予測した波形が（ｂ）、区間〔0,512〕
における被測定波形と合成波形の差が（ｃ）である。ス
ペクトルの縦軸は相対振幅、横軸は相対的な周波数1/T
を示す。また、比較のためにFFTによるスペクトル推定
結果を（ｆ）に示す。
図６は、図５の被測定波形にランダム雑音を加えた場
合について図５の場合と同様な分析処理を行なった結果
を示したもので、雑音の影響および雑音抑圧効果が分
る。（ａ）はランダム雑音を加えられた被測定波形、
（ｄ）は被測定波形（ノイズなし）のスペクトルを示し
ている。以下の説明は図５に準じる。なお、図５および
図６の分析は、ｎ＝4,N＝30とし、Ｔは10から128まで１
ずつ変化させた。
合について図５の場合と同様な分析処理を行なった結果
を示したもので、雑音の影響および雑音抑圧効果が分
る。（ａ）はランダム雑音を加えられた被測定波形、
（ｄ）は被測定波形（ノイズなし）のスペクトルを示し
ている。以下の説明は図５に準じる。なお、図５および
図６の分析は、ｎ＝4,N＝30とし、Ｔは10から128まで１
ずつ変化させた。
本発明による波形データの非調和分析および合成処理
方法の効果を最も端的に示すものは、波形の周期性検出
の精度と局所性に関するものである。本分析法によれば
最低１波長の部分波形によって周期検出ができるため
に、例えば指数関数的に減衰する複数の異なる周期をも
った正弦波形から成る減衰波形について、分析点を移動
しつつ各正弦波の振幅を求めれば、これら振幅の変化か
ら個々の正弦波の減衰量を精度よく求めることができ
る。これによって減衰振動の複雑な変化を各周波数成分
の振幅変化として観測することができる。これは楽器音
の分析や地震波形の分析に応用できる。
方法の効果を最も端的に示すものは、波形の周期性検出
の精度と局所性に関するものである。本分析法によれば
最低１波長の部分波形によって周期検出ができるため
に、例えば指数関数的に減衰する複数の異なる周期をも
った正弦波形から成る減衰波形について、分析点を移動
しつつ各正弦波の振幅を求めれば、これら振幅の変化か
ら個々の正弦波の減衰量を精度よく求めることができ
る。これによって減衰振動の複雑な変化を各周波数成分
の振幅変化として観測することができる。これは楽器音
の分析や地震波形の分析に応用できる。
なお、本発明による波形データの非調和分析および合
成処理方法における分析においては、分析周期（あるい
は周波数）を帯域に分けて低い周波数帯域から順次分析
を進めることができる。この場合は帯域分割数をｓとす
ると演算時間は約1/sになる。また演算部（例えば図１
のD2〜D3）をｓ組並列に配置して、ｓ個の帯域をそれぞ
れ同時に演算するならば、その部分の演算時間は約1/s2
になる。これらはいずれも本発明に変更を加えるもので
はない。
成処理方法における分析においては、分析周期（あるい
は周波数）を帯域に分けて低い周波数帯域から順次分析
を進めることができる。この場合は帯域分割数をｓとす
ると演算時間は約1/sになる。また演算部（例えば図１
のD2〜D3）をｓ組並列に配置して、ｓ個の帯域をそれぞ
れ同時に演算するならば、その部分の演算時間は約1/s2
になる。これらはいずれも本発明に変更を加えるもので
はない。
Claims (3)
- 【請求項１】物理的な波形を非調和周波数成分に分析す
る方法であって、 （ａ）所与の区間の物理的波形を対応するデジタル波形
データに変換し、 （ｂ）デジタル波形データに第１の所定周期を有する第
１の正弦波形と第１の余弦波形をそれぞれ乗算して第１
及び第２の積を求め、 （ｃ）前記第１の所定周期の整数倍であり前記所与の区
間に等しいかそれよりも小さい第１の所定区間に亘って
第１及び第２の積をそれぞれ加算して第１及び第２の和
値を求め、 （ｄ）前記第１及び第２の和値に基づきそれぞれ前記第
１の正弦波形及び第１の余弦波形の振幅を求め、 （ｅ）このようにして求めたそれぞれの振幅と前記第１
の所定周期を有する正弦波形及び余弦波形を前記デジタ
ル波形データから差し引いて残差波形データを求め、 （ｆ）それぞれ第２乃至第（ｎ）の所定周期を有する第
２乃至第（ｎ）の正弦波形及び余弦波形を用いて前記デ
ジタル波形データに対し前記ステップ（ｂ）乃至（ｅ）
を順次実行することにより（ｎ−１）個の残差波形デー
タを求め、 （ｇ）前記（ｎ）個の残差波形データのパワーを調べる
ことにより前記パワーが最小となる前記正弦波形及び余
弦波形を第１の非調和周波数成分とし、 （ｈ）最小のパワーの残差波形データを与える前記正弦
波形及び余弦波形を前記デジタル波形データから差し引
いて第１の残差波形を求め、 （ｉ）前記第１の残差波形に対し前記ステップ（ｂ）乃
至（ｈ）を順次実行することにより第２の非調和周波数
成分を求め、 （ｊ）前記第１の残差波形から前記第２の非調和周波数
成分を差し引いて第２の残差波形を求め、 （ｋ）前記第２の残差波形に対し前記ステップ（ｂ）乃
至（ｈ）を順次実行することにより第３の非調和周波数
成分を求める、という手順を繰り返すことにより、第
（ｎ）の残差波形から第（ｎ＋１）の非調和周波数成分
を求めることよりなる物理的な波形を非調和周波数成分
に分析する方法。 - 【請求項２】（ｎ）個の残差波形データのパワーを調べ
る前記ステップ（ｇ）は（ｎ）個の残差波形データの二
乗値を所定の区間に亘って加算することよりなる請求項
１に記載の方法。 - 【請求項３】（ｎ）個の残差波形データのパワーを調べ
る前記ステップ（ｇ）は（ｎ）個の残差波形データの絶
対値を所定の区間に亘り加算することよりなる請求項１
に記載の方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP51787194A JP2714880B2 (ja) | 1993-02-02 | 1994-02-02 | 物理的な波形を非調和周波数成分に分析する方法 |
Applications Claiming Priority (5)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP5120993 | 1993-02-02 | ||
| JP5-51209 | 1993-02-02 | ||
| JP5-65865 | 1993-02-16 | ||
| JP6586593 | 1993-02-16 | ||
| JP51787194A JP2714880B2 (ja) | 1993-02-02 | 1994-02-02 | 物理的な波形を非調和周波数成分に分析する方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| JP2714880B2 true JP2714880B2 (ja) | 1998-02-16 |
Family
ID=27294244
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP51787194A Expired - Fee Related JP2714880B2 (ja) | 1993-02-02 | 1994-02-02 | 物理的な波形を非調和周波数成分に分析する方法 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP2714880B2 (ja) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2020178919A1 (ja) * | 2019-03-01 | 2020-09-10 | 三菱電機株式会社 | 信号解析装置、信号解析システムおよび信号解析方法 |
-
1994
- 1994-02-02 JP JP51787194A patent/JP2714880B2/ja not_active Expired - Fee Related
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2020178919A1 (ja) * | 2019-03-01 | 2020-09-10 | 三菱電機株式会社 | 信号解析装置、信号解析システムおよび信号解析方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US5029509A (en) | Musical synthesizer combining deterministic and stochastic waveforms | |
| US5749047A (en) | Method and apparatus for recognizing a test signal and determining signal transfer characteristics therefrom | |
| JP5140676B2 (ja) | 計算による音楽テンポの推定 | |
| EP0153787B1 (en) | System of analyzing human speech | |
| JP3704336B2 (ja) | 波形検出装置とそれを使用した状態監視システム | |
| EP1974225B1 (en) | Signal analyser | |
| RU2010139018A (ru) | Устройство и способ преобразования звукового сигнала в параметрическое представление, устройство и способ модификации параметрического представления, устройство и способ синтеза параметрического представления звукового сигнала | |
| EP3121608B1 (en) | Method of modeling characteristics of a non linear system. | |
| JPH0844394A (ja) | 励起パラメータの評価 | |
| KR100911685B1 (ko) | 코히어런트하지 않게 샘플링된 데이타의 파워 스펙트럼을측정하기 위한 저누설 방법 | |
| Lopatinskaia et al. | Monitoring varying speed machinery vibrations—I. The use of non-stationary recursive filters | |
| JP2714880B2 (ja) | 物理的な波形を非調和周波数成分に分析する方法 | |
| US5862516A (en) | Method of non-harmonic analysis and synthesis of wave data | |
| JPH09243679A (ja) | 任意区間波形を用いた非調和的周波数分析法 | |
| CN109975867B (zh) | 一种频率域信号混叠的地震数据拓频方法 | |
| JPH05281996A (ja) | ピッチ抽出装置 | |
| CN107124678B (zh) | 一种音频谐波失真的测量系统 | |
| JP3599546B2 (ja) | 波形検出装置およびその装置を利用した状態監視システム | |
| Wu et al. | Time domain averaging based on fractional delay filter | |
| JP2001222289A (ja) | 音響信号分析方法及び装置並びに音声信号処理方法及び装置 | |
| Ribeiro et al. | Non-stationary analysis and noise filtering using a technique extended from the original prony method | |
| Maula et al. | Spectrum identification of peking as a part of traditional instrument of gamelan | |
| JPH1020886A (ja) | 波形データに存在する調和波形成分の検出方式 | |
| JP3646702B2 (ja) | 周波数分析装置 | |
| KR0128851B1 (ko) | 극성이 다른 가변길이 듀얼 임펄스의 스펙트럼 하모닉스 매칭에 의한 피치 검출 방법 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| LAPS | Cancellation because of no payment of annual fees |