|
|
| ミランコヴィッチメニューへ戻る 序章 本書は時系列のウェーブレット解析を扱っています。 ウェーブレット変換の基本的な考え方は、多くの時系列の特異性に対応しており、それらの主要な特性である平均、分散、期間、振幅、高調波成分の位相の経時的な進化を示しています。 様々な知識分野で研究されているプロセスの大部分は、上記のような特徴を持っています。 本書の目的は、PyWaveletsライブラリを用いた時系列の連続ウェーブレット変換の手法を説明することである。 ちょっとした歴史 70年代後半、地球物理学者のD.モーレットは、短期間に高周波成分（地震活動）を、長期間に低周波成分（地殻の静穏状態）を含む地震センサーからの信号を解析するという問題に直面しました。 フーリエ変換ウィンドウでは、高周波成分と低周波成分のどちらかを分析することができますが、両方の成分を一度に分析することはできません。 そこで、低周波では窓関数の幅を大きくし、高周波では幅を小さくする解析方法を提案した。 両方の生成関数（いわゆるスケーリング関数、scalet）を伸張（圧縮）し、時間的にシフトさせることにより、新しいウィンドウ変換を得た。この生成関数はモーレットウェーブレットと呼ばれていました。 ウェーブレットD.モーレット   連続ウェーブレット変換(CWT) 1 レベルウェーブレット変換。 coefs, frequencies =pywt.cwt(data, scales, wavelet) data : (array_like):-入力信号。 scales (array_like):- 利用するウェーブレットのスケール．比 f =scale2frequency(scale, wavelet)/sampling_period を使用して、f がどの物理周波数であるかを判断することができます。ここで、fはsampling_periodを秒単位で指定した場合のヘルツ数です。 wavelet (Wavelet オブジェクトまたは名前): - 使用するウェーブレット． sampling_period (float): - 出力周波数のサンプリング周期 (オプションパラメータ) coefsのために計算された値は、サンプリング周期の選択に依存しない（すなわち、サンプリング周期によってスケーリングされない）。 coefs (array_like ): - 指定されたスケール，ウェーブレットに対する入力信号の連続ウェーブレット変換． frequencies (array_like ): - サンプリング周期の単位が秒で指定されている場合は，ヘルツ単位で出力されます．それ以外の場合は、サンプリング周期を1とする。 注意：係数配列のサイズは，入力配列の長さと与えられたスケールの長さに依存します． 利用可能なcwt互換のウェーブレット名のリストを取得します。 >>>> import pywt >>>> wavlist = pywt.wavelist(kind='continuous') >>> wavlistを使用しています。 ['cgau1', 'cgau2', 'cgau3', 'cgau4', 'cgau5', 'cgau6', 'cgau7', 'cgau8', 'cmor', 'fbsp', 'gaus1', 'gaus2', 'gaus3', 'gaus4', 'gaus5', 'gaus6', 'gaus7', 'gaus8', 'mexh', 'moll', 'shan'] ウェーブレット解析では、親ウェーブレットの選択が重要な作業の一つです。したがって、各ウェーブレットの選択の前に、以下のプログラムを使用してウェーブレットの動的特性を決定することが単純に必要である。 リスト1. 次に、上記のプログラムを用いて、ウェーブレット関数とその性質について検討する。 メキシコの帽子のウェーブレット "mexh"   モーレのウェーブレット  複素モーレットウェーブレット "cmorB-C"  B - 帯域幅、C - 中心周波数。 以下、（特に表示せずに）B,Cの値は浮動小数点で与えます。  ガウスウェーブレット "gausP"  de: P - ウェーブレットの名前に挿入される 1 から 8 までの整数，例えば "gaus8"; C - 組み込みの正規化定数．  複素ガウスウェーブレット "cgauP"  ここで，Pはウェーブレットの名前に挿入される1から8までの整数，例えば "cgaus8 "であり，ウェーブレット関数の微分の次数に対応します．  シャノンウェーブレット "shanB-C"  ここで: B - 帯域幅; C - 中心周波数.  PyWaveletsのCWTは離散データ、ウェーブレット積分のサンプルを用いた畳み込みに適用されます。スケールが小さすぎると、'cmor1.5-1.0' ウェーブレットのグラフに示されているように、サンプリングが不十分で平滑化につながる離散フィルタが得られます。 グラフの左列は、異なるスケールでの畳み込みに使用される離散フィルタを示しています。右の列は各フィルタの対応するフーリエパワースペクトルです。 cmor1.5-1.0' グラフの縮尺1と2では、平滑化がナイキスト制約の違反によるものであることがわかります。上記の現象は、CとBを選択したときに他のウェーブレットでも発生する可能性があるので、ウェーブレットを扱うときには、プログラム - リスト1 - を使用する必要があります。 所定のスケールを特定の信号周波数に関連付けるためには、信号のサンプリング周期を知る必要があります。pywt.scale2frequency()を使用すると、スケールのリストを適切な周波数に変換することができます。正しいスケールの選択は選択されたウェーブレットに依存するので、 pywt.scale2frequency()を使用して信号の周波数範囲を把握してください。 >>>> import pywt >>> 100Hz サンプリング >>>> frequencies = pywt.scale2frequency('cmor1.5-1.0', [1, 2, 3, 4]) / dt >>> 周波数 array([100., 50., 33.33333, 25. ]) ウェーブレット - CWT PyWaveletsを用いた時系列解析 エルニーニョデータセットは、エルニーニョの追跡に使用される時系列データセットで、1871年から1997年までの四半期ごとの海面水温の測定値が含まれています。 スケールグラムの威力を理解するために、エルニーニョのデータセットについて、元の時系列データとそのフーリエ変換を用いて、スケールグラムを可視化してみましょう。 時系列のウェーブレット解析を行うには、以下の手順をポイントバイポイントで行う必要があります。 1. 親ウェーブレットの選択：複素モーレットウェーブレット "cmorB-C "を選択します。  スループット-Ｂおよび中心周波数-Ｃは、次のステップで決定されます。 2.時系列にdt=0.25を取り、中心周波数を決定します。 インポートpywt dt = 0.25 # 4 Hz スケール=範囲(1,4) 周波数 = pywt.scale2frequency('cmor1.0-0.5', scale) / dt プリント(周波数) [2. 1. 0.66666667] 3. 母ウェーブレット cmor1.0-0.5のフーリエ変換を求めます(リスト1を使用)。 連続ウェーブレットは、全範囲[-8.0, 8.0]にわたって評価されます。  4. 時系列のデータを選択します。 pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") データは1871年から1997年までの四半期ごとの海面水温測定値です。このデータの場合：t0=1871 dt=0.25 5. シグナルとその移動平均線を持つ時系列を1つのグラフにプロットします。 リスト2  6. 時系列からフーリエ変換を行い、スペクトルを修正します。 リスト3  7. スケールの定義: scales = arange(1, 128); levels = [2**-4, 2**-3, 2**-2, 2**-1, 2**0, 2**1, 2**2, 2**3]。 8. スケーリング図を作ってみましょう。 リスト4.  散布図は、スペクトルパワーの大部分が2-8年の期間に集中していることを示しており、これは0.125-0.5Hzの周波数に対応しています。 フーリエスペクトルでは、モディシティもまた、これらの周波数値の周りに集中しています。しかし、ウェーブレット変換は時間的な情報も与えてくれますが、フーリエ変換はそうではありません。 例えば、1920年以前は変動が多かったのに対し、1960年から1990年まではあまり変動がなかったことが縮尺図でわかります。また、時間の経過とともに、短い期間から長い期間へと移行していることもわかります。 スケールグラムモジュールの使用 Scaleogramは、PyWaveletsライブラリをベースに構築された連続ウェーブレット変換を用いた1次元データ解析のための便利なツールです。このモジュールは、データを高速に分析するための信頼性の高いツールを提供するために設計されています。 スケールグラムには以下のような特徴があります。 初心者のためのシンプルなコールサイン 読みやすい軸とクリーンなmatplotlibの統合 スケーリング、スペクトルフィルタ、カラーバーの実装など、多くのオプションがあります。 マーキングによる周期性と周波数単位のサポート 速度は、変換に N * log(N) アルゴリズムを使用します。 移植性: python2.7 と python3.7 でテストしました。 詳細なエラーメッセージと例を含むドキュメント しかし、使用例ではデータ参照が正しく書かれていません。 nino3 = "http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat" 以下のような警告が表示されます。 nino3 = pd.read_table(nino3) と警告が出てしまいます。 警告（警告モジュールから）。 ファイル "C:/Users/User/Desktop/wavelet/pr.py", 7行目 nino3 = pd.read_table(nino3) FutureWarning: read_table is deprecated, use read_csv instead, pass sep='\t' データ参照は次のように記述する必要があります。 url="http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat" nino3 = pd.read_csv(url, sep = "|") スカログラムの使い方を示し、前の例と結果を比較するために、1871年から1997年までの四半期ごとの海面水温の測定結果と同じデータを使用します。これらのデータの場合：t0=1871 dt=0.25 リスト5  カラーパレットを除いて、スキャッタグラムと結果のスキャングラムを比較すると、それらは同一であり、したがって、同じ時系列のダイナミクスを示します。 リスト5はリスト3よりもシンプルで高速です。 周波数スペクトルやフーリエパワースペクトルでは時系列のダイナミクスに関する追加情報が得られない場合は、スカログラムの使用が好ましい。 結論 PyWaveletsライブラリを用いた時系列のウェーブレット解析の例を示します。 Tags. ウェーブレット解析 スペクトル解析 フーリエ級数 スケーリング じかんちえん ハブ。 パイソン 数学 Windows開発 ポピュラーサイエンス 物理学 |