👩🏾‍🔬 🥛 👩 Wavelet-Analyse. Teil 2 👨‍🚒 😅 🚹

Einführung

Diese Veröffentlichung behandelt die Wavelet-Analyse von Zeitreihen. Die Hauptidee der Wavelet-Transformation entspricht den Besonderheiten vieler Zeitreihen, die die zeitliche Entwicklung ihrer Hauptmerkmale demonstrieren - Mittelwert, Varianz, Perioden, Amplituden und Phasen harmonischer Komponenten. Die überwiegende Mehrheit der in verschiedenen Wissensgebieten untersuchten Prozesse weist die oben genannten Merkmale auf.

Der Zweck dieser Veröffentlichung ist es, die Methode der kontinuierlichen Wavelet-Transformation von Zeitreihen unter Verwendung der PyWavelets-Bibliothek zu beschreiben.

Ein bisschen Geschichte

Geophysiker D. Morlaix in den späten 70er Jahren des 20. Jahrhunderts. Ich stieß auf das Problem, Signale von seismischen Sensoren zu analysieren, die über einen kurzen Zeitraum eine Hochfrequenzkomponente (seismische Aktivität) und über einen langen Zeitraum Niederfrequenzkomponenten (ruhiger Zustand der Erdkruste) enthielten. Mit der Fenster-Fourier-Transformation können Sie entweder die Hochfrequenzkomponente oder die Niederfrequenzkomponente analysieren, jedoch nicht beide Komponenten gleichzeitig.

Daher wurde ein Analyseverfahren vorgeschlagen, bei dem die Breite der Fensterfunktion für niedrige Frequenzen zunahm und für hohe Frequenzen abnahm. Die neue Fenstertransformation wurde als Ergebnis der Dehnung (Komprimierung) und des Zeitversatzes einer generierenden Funktion (sogenannte Skalierungsfunktion - Skalierungsfunktion, Skalette) erhalten. Diese Erzeugungsfunktion wurde von D. Morlaix als Wavelet bezeichnet.

Wavelet D. Morlaix

from pylab import* import scaleogram as scg axes = scg.plot_wav('cmor1-1.5', figsize=(14,3)) show()

Kontinuierliche Wavelet-Transformation (CWT)

Einstufige Wavelet-Transformation:

Coefs, Frequenzen = pywt.cwt (Daten, Skalen, Wavelet)

wo:

Daten: (array_like) -Eingabesignal;

skaliert (array_like) : - Die zu verwendende Wavelet-Skala. Sie können das Verhältnis f = scale2frequency (scale, Wavelet) / Sampling_period verwenden und bestimmen, welche physikalische Frequenz f ist. Hier f in Hertz bei Abtastperiode in Sekunden;

Wavelet (Wavelet-Objekt oder -Name) : - Zu verwendendes Wavelet;

Abtastperiode (float): - Abtastperiode für Ausgangsfrequenzen (optionaler Parameter). Die für Coefs berechneten Werte sind unabhängig von der Abtastperiode (d. H. Skalen werden nicht durch Periodenabtastung skaliert).

coefs (array_like) : - Kontinuierliches Wavelet - Transformation des Eingangssignals für eine bestimmte Skala und ein bestimmtes Wavelet;

Frequenzen (array_like) : - Wenn die Einheit der Abtastperiode die zweite ist und eingestellt ist, werden die Frequenzen in Hertz angezeigt. Andernfalls wird eine Abtastperiode von 1 angenommen.

Hinweis: Die Größe der Koeffizientenarrays hängt von der Länge des Eingabearrays und den Längen der angegebenen Skalen ab.

Holen Sie sich eine Liste der verfügbaren cwt-kompatiblen Wavelet-Namen:

 >>> 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 ',' morl ',' shan ']

Für die Wavelet-Analyse ist die Wahl des Mutter-Wavelets eine der Hauptaufgaben. Daher muss vor jeder Wavelet-Auswahl einfach das folgende Programm verwendet werden, mit dem Sie die dynamischen Eigenschaften des Wavelets bestimmen können:

Listing 1

 import numpy as np import pywt import matplotlib.pyplot as plt vav='cmor1.5-1.0' wav = pywt.ContinuousWavelet(vav) #  ,      print("       [{}, {}]".format( wav.lower_bound, wav.upper_bound)) width = wav.upper_bound - wav.lower_bound scales = [1, 2, 3, 4, 10, 15] max_len = int(np.max(scales)*width + 1) t = np.arange(max_len) fig, axes = plt.subplots(len(scales), 2, figsize=(12, 6)) for n, scale in enumerate(scales): #       cwt int_psi, x = pywt.integrate_wavelet(wav, precision=10) step = x[1] - x[0] j = np.floor( np.arange(scale * width + 1) / (scale * step)) if np.max(j) >= np.size(int_psi): j = np.delete(j, np.where((j >= np.size(int_psi)))[0]) j = j.astype(np.int) # normalize int_psi     int_psi /= np.abs(int_psi).max() #     filt = int_psi[j][::-1] # CWT          #          . nt = len(filt) t = np.linspace(-nt//2, nt//2, nt) axes[n, 0].plot(t, filt.real, t, filt.imag) axes[n, 0].set_xlim([-max_len//2, max_len//2]) axes[n, 0].set_ylim([-1, 1]) axes[n, 0].text(50, 0.35, 'scale = {}'.format(scale)) f = np.linspace(-np.pi, np.pi, max_len) filt_fft = np.fft.fftshift(np.fft.fft(filt, n=max_len)) filt_fft /= np.abs(filt_fft).max() axes[n, 1].plot(f, np.abs(filt_fft)**2) axes[n, 1].set_xlim([-np.pi, np.pi]) axes[n, 1].set_ylim([0, 1]) axes[n, 1].set_xticks([-np.pi, 0, np.pi]) axes[n, 1].set_xticklabels([r'$-\pi$', '0', r'$\pi$']) axes[n, 1].grid(True, axis='x') axes[n, 1].text(np.pi/2, 0.5, 'scale = {}'.format(scale)) axes[n, 0].set_xlabel('time (samples)') axes[n, 1].set_xlabel('frequency (radians)') axes[0, 0].legend(['real', 'imaginary'], loc='upper left') axes[0, 1].legend(['Power'], loc='upper left') axes[0, 0].set_title('filter: wavelet - %s'%vav) axes[0, 1].set_title(r'|FFT(filter)|$^2$') plt.show()

Als nächstes werden wir Wavelet-Funktionen und ihre Eigenschaften mit dem obigen Programm betrachten:

Mexikanischer Hut Wavelet Mexh:

v a r p h i l e f t (t r i g h t) = f r a c 2 s q r t 3 s q r t [4] p i c d o t e x p^{- f r a c t^{2} 2} c d o t l e f t (1 - t^{2} r i g h t)

$\ varphi \ left (t \ right) = \ frac {2} {\ sqrt {3} \ sqrt [4] {\ pi}} \ cdot exp ^ {- \ frac {t ^ {2}} {2} } \ cdot \ left (1-t ^ {2} \ right)$

Wavelet "Morl" Morlaix:

v a r p h i l e f t (t r i g h t) = e x p^{- f r a c t^{2} 2} c d o t c o s (5 t)

$\ varphi \ left (t \ right) = exp ^ {- \ frac {t ^ {2}} {2}} \ cdot cos (5t)$

Komplexes Morlet-Wavelet "cmorB-C"

v a r p h i l e f t (t r i g h t) = f r a c 1 s q r t p i B e x p^{- f r a c t^{2} B} c d o t e x p^{j 2 p i C t}

$\ varphi \ left (t \ right) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

wobei B die Bandbreite ist; C ist die Mittenfrequenz.

Nachfolgend (ohne besondere Angabe) werden die Werte von B, C mit einem Gleitkomma gesetzt.

Gaußsche Wavelets "gausP"

v a r p h i l e f t (t r i g h t) = C c d o t e x p^{- t^{2}}

$\ varphi \ left (t \ right) = C \ cdot exp ^ {- t ^ {2}}$

wobei: P - eine ganze Zahl von 1 bis 8 in den Namen des Wavelets eingefügt wird, zum Beispiel "gaus8"; C-eingebaute Normalisierungskonstante.

Integrierte Gaußsche Wavelets “cgauP”

v a r p h i l e f t (t r i g h t) = C c d o t e x p^{- t^{2}} c d o t e x p^{- j t}

$\ varphi \ left (t \ right) = C \ cdot exp ^ {- t ^ {2}} \ cdot exp ^ {- jt}$

wobei: P - eine ganze Zahl von 1 bis 8 in den Namen des Wavelets eingefügt wird, zum Beispiel "gaus8" und der Reihenfolge der Ableitung der Wavelet-Funktion entspricht; C-eingebaute Normalisierungskonstante.

Shannon Wavelets "shanB-C"

v a r p h i l e f t (t r i g h t) = s q r t B c d o t f r a c s i n (p i B t) p i B T c d o t e x p^{j 2 p i C t}

$\ varphi \ left (t \ right) = \ sqrt {B} \ cdot \ frac {sin (\ pi B t)} {\ pi B T} \ cdot exp ^ {j 2 pi C t}$

wobei: B die Bandbreite ist; C ist die Mittenfrequenz.

CWT in PyWavelets gilt für diskrete Daten - Windungen mit Abtastwerten des Wavelet-Integrals. Wenn der Maßstab zu klein ist, erhalten wir ein diskretes Filter mit unzureichender Abtastung, was zu einer Glättung führt, wie in der Grafik für das Wavelet 'cmor1.5-1.0' gezeigt.

In der linken Spalte zeigen die Grafiken diskrete Filter, die bei der Faltung für verschiedene Skalen verwendet werden. Die rechte Spalte enthält die entsprechenden Fourier-Leistungsspektren jedes Filters. Mit den Skalen 1 und 2 für den Graphen 'cmor1.5-1.0' ist ersichtlich, dass eine Glättung aufgrund einer Verletzung der Nyquist-Einschränkung auftritt. Das angegebene Phänomen kann bei der Auswahl von C und B in anderen Wavelets auftreten. Wenn Sie also mit Wavelets arbeiten, müssen Sie das Programm Listing 1 verwenden.

Um eine bestimmte Skala mit einer bestimmten Signalfrequenz in Beziehung zu setzen, muss die Signalabtastperiode bekannt sein. Mit der Funktion pywt.scale2frequency () können Sie die Liste der Skalen in die entsprechenden Frequenzen konvertieren. Die richtige Auswahl der Skalen hängt vom ausgewählten Wavelet ab. Daher sollte pywt.scale2frequency () verwendet werden, um eine Vorstellung vom entsprechenden Frequenzbereich des Signals zu erhalten.

 >>> import pywt >>> dt = 0.01 # 100 Hz sampling >>> frequencies = pywt.scale2frequency('cmor1.5-1.0', [1, 2, 3, 4]) / dt >>> frequencies

Array ([100., 50., 33.33333333, 25.])

Wavelet - Zeitreihenanalyse mit CWT PyWavelets

Der el-Nino-Datensatz ist ein Zeitreihendatensatz zur Verfolgung von El Nino und enthält vierteljährliche Messungen der Meeresoberflächentemperatur von 1871 bis 1997. Um die Leistungsfähigkeit eines Skalendiagramms zu verstehen, visualisieren wir es für einen el-Nino-Datensatz zusammen mit den ursprünglichen Zeitreihendaten und seiner Fourier-Transformation.

Für die Wavelet-Analyse von Zeitreihen müssen die folgenden Schritte ausgeführt werden:

1. Wählen Sie das Mutter-Wavelet aus: Wählen Sie das komplexe Morlet-Wavelet "cmorB-C" aus:

v a r p h i l e f t (t r i g h t) = f r a c 1 s q r t p i B e x p^{- f r a c t^{2} B} c d o t e x p^{j 2 p i C t}

$\ varphi \ left (t \ right) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

Bandbreite - B und Mittenfrequenz - Aus der wir im nächsten Schritt bestimmen werden.

2. Bestimmen Sie die Mittenfrequenz, indem Sie für unsere Zeitreihen dt = 0,25 nehmen:

 import pywt dt = 0.25 # 4 Hz sampling scale=range(1,4) frequencies = pywt.scale2frequency('cmor1.0-0.5', scale) / dt print(frequencies)

[2. 1. 0,66666667]

3. Wir finden die Fourier-Transformation des Mutter-Wavelets cmor1.0-0.5 (wir verwenden Listing 1):
Kontinuierliches Wavelet wird über den gesamten Bereich ausgewertet [-8.0, 8.0]

4. Wählen Sie die Daten für die Zeitreihe aus:

 pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|")

Die Daten sind vierteljährliche Messungen der Meeresoberflächentemperatur von 1871 bis 1997. Für diese Daten gilt: t0 = 1871 dt = 0,25

5. Wir erstellen eine Zeitreihe mit einem Signal und seinem gleitenden Durchschnitt auf einem Diagramm:

Listing 2:

 from numpy import * from scipy import * import pandas as pd from pylab import * import pywt def get_ave_values(xvalues, yvalues, n = 5): signal_length = len(xvalues) if signal_length % n == 0: padding_length = 0 else: padding_length = n - signal_length//n % n xarr = array(xvalues) yarr = array(yvalues) xarr.resize(signal_length//n, n) yarr.resize(signal_length//n, n) xarr_reshaped = xarr.reshape((-1,n)) yarr_reshaped = yarr.reshape((-1,n)) x_ave = xarr_reshaped[:,0] y_ave = nanmean(yarr_reshaped, axis=1) return x_ave, y_ave def plot_signal_plus_average(time, signal, average_over = 5): fig, ax = subplots(figsize=(15, 3)) time_ave, signal_ave = get_ave_values(time, signal, average_over) ax.plot(time, signal, label='') ax.plot(time_ave, signal_ave, label = '   (n={})'.format(5)) ax.set_xlim([time[0], time[-1]]) ax.set_ylabel(' ', fontsize=18) ax.set_title(' +   ', fontsize=18) ax.set_xlabel('', fontsize=18) ax.legend() show() df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") N = df_nino.shape[0] t0=1871 dt=0.25 time = arange(0, N) * dt + t0 signal = df_nino.values.squeeze() scales = arange(1, 128) plot_signal_plus_average(time, signal)

6. Wir führen die Fourier-Transformation und die Spektrum-Modalität aus den Zeitreihen durch:

Listing 3:

 from numpy import * from scipy import * import pandas as pd from pylab import * import pywt def get_ave_values(xvalues, yvalues, n = 5): signal_length = len(xvalues) if signal_length % n == 0: padding_length = 0 else: padding_length = n - signal_length//n % n xarr = array(xvalues) yarr = array(yvalues) xarr.resize(signal_length//n, n) yarr.resize(signal_length//n, n) xarr_reshaped = xarr.reshape((-1,n)) yarr_reshaped = yarr.reshape((-1,n)) x_ave = xarr_reshaped[:,0] y_ave = nanmean(yarr_reshaped, axis=1) return x_ave, y_ave def get_fft_values(y_values, T, N, f_s): f_values = linspace(0.0, 1.0/(2.0*T), N//2) fft_values_ = fft(y_values) fft_values = 2.0/N * abs(fft_values_[0:N//2]) return f_values, fft_values def plot_fft_plus_power(time, signal): dt = time[1] - time[0] N = len(signal) fs = 1/dt fig, ax = subplots(figsize=(15, 3)) variance = std(signal)**2 f_values, fft_values = get_fft_values(signal, dt, N, fs) fft_power = variance * abs(fft_values) ** 2 # FFT power spectrum ax.plot(f_values, fft_values, 'r-', label='FFT ') ax.plot(f_values, fft_power, 'k--', linewidth=1, label=' ') ax.set_xlabel('[ / ]', fontsize=18) ax.set_ylabel('', fontsize=18) ax.legend() show() df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") N = df_nino.shape[0] t0=1871 dt=0.25 time = arange(0, N) * dt + t0 signal = df_nino.values.squeeze() scales = arange(1, 128) plot_fft_plus_power(time, signal)

7. Wir bestimmen die Skalen: Skalen = Bereich (1, 128); Stufen = [2 ** - 4, 2 ** - 3, 2 ** - 2, 2 ** - 1, 2 ** 0, 2 ** 1, 2 ** 2, 2 ** 3].

8. Erstellen Sie ein Skalendiagramm:

Listing 4:

 from numpy import * import pandas as pd from pylab import * import pywt def plot_wavelet(time, signal, scales, waveletname = 'cmor1.0-0.4', cmap = plt.cm.seismic, title = '-( ) ', ylabel = ' ()', xlabel = ''): dt = time[1] - time[0] [coefficients, frequencies] = pywt.cwt(signal, scales, waveletname, dt) power = (abs(coefficients)) ** 2 period = 1. / frequencies levels = [2**-4 , 2**-3 , 2**-2 , 2**-1 , 2**0 , 2**1 , 2**2 , 2**3] contourlevels = log2(levels) fig, ax = subplots(figsize=(15, 10)) im = ax.contourf(time, log2(period), log2(power), contourlevels, extend='both',cmap=cmap) ax.set_title(title, fontsize=20) ax.set_ylabel(ylabel, fontsize=18) ax.set_xlabel(xlabel, fontsize=18) yticks = 2**arange(np.ceil(log2(period.min())), ceil(log2(period.max()))) ax.set_yticks(log2(yticks)) ax.set_yticklabels(yticks) ax.invert_yaxis() ylim = ax.get_ylim() ax.set_ylim(ylim[0], -1) cbar_ax = fig.add_axes([0.95, 0.5, 0.03, 0.25]) fig.colorbar(im, cax=cbar_ax, orientation="vertical") show() df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") N = df_nino.shape[0] t0=1871 dt=0.25 time = arange(0, N) * dt + t0 signal = df_nino.values.squeeze() scales = arange(1, 128) plot_wavelet(time, signal, scales)

Das Skalendiagramm zeigt, dass der größte Teil der Spektrumleistung über einen Zeitraum von 2-8 Jahren konzentriert ist, was einer Frequenz von 0,125 - 0,5 Hz entspricht. Im Fourier-Spektrum konzentriert sich die Mode auch um diese Frequenzwerte. Die Wavelet-Transformation liefert uns jedoch auch temporäre Informationen, die Fourier-Transformation jedoch nicht.

Zum Beispiel sehen wir im Skalendiagramm, dass es vor 1920 viele Schwankungen gab, während es zwischen 1960 und 1990 nicht so viele gab. Wir können auch sehen, dass es im Laufe der Zeit eine Verschiebung von kürzeren zu längeren Zeiträumen gibt.

Verwendung des Skalogrammmoduls

Das Skalogramm ist ein praktisches Werkzeug zur Analyse von 1D-Daten mit kontinuierlicher Wavelet-Transformation basierend auf der PyWavelets-Bibliothek. Das Modul bietet ein zuverlässiges Werkzeug für die schnelle Datenanalyse.

Das Skalogramm hat folgende Funktionen:

einfache Anrufersignatur für Anfänger
lesbare Achsen und reine Matplotlib-Integration
viele Optionen zum Zoomen, Spektrumsfilter, Farbbalkenimplementierung usw.
Unterstützung für Frequenz und Frequenzeinheiten gemäß Kennzeichnung
Geschwindigkeit, verwendet den N * log (N) -Algorithmus für Transformationen
Portabilität: Getestet mit Python2.7 und Python3.7
detaillierte Fehlermeldungen und Dokumentation mit Beispielen

In den Verwendungsbeispielen wird der Datenzugriff jedoch nicht korrekt aufgezeichnet:

 nino3 = "http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat"

Die folgende Warnung wird angezeigt:

nino3 = pd.read_table (nino3)
was zu einer Warnung führt:
Warnung (vom Warnmodul):
Datei “C: /Users/User/Desktop/wavelet/pr.py”, Zeile 7
nino3 = pd.read_table (nino3)
FutureWarning: read_table ist veraltet. Verwenden Sie stattdessen read_csv und übergeben Sie sep = '\ t'.

Der Datenzugriff sollte folgendermaßen geschrieben werden:

 url="http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat" nino3 = pd.read_csv(url, sep = "|")

Um die Verwendung des Skalogramms zu zeigen und die Ergebnisse mit den Ergebnissen des vorherigen Beispiels zu vergleichen, verwenden wir dieselben Daten - für vierteljährliche Messungen der Meeresoberflächentemperatur von 1871 bis 1997. Für diese Daten gilt: t0 = 1871 dt = 0,25

Listing 4:

 import pandas as pd import pywt from numpy import* import scaleogram as scg from pylab import* url="sst_nino3.dat" nino3 = pd.read_csv(url, sep = "|") data = nino3.values.squeeze() N = data.size t0 = 1871; dt = 0.25 time = t0 + arange(len(data))*dt wavelet = 'cmor1-0.5' ax = scg.cws(time, data, scales=arange(1, 128), wavelet=wavelet, figsize=(14, 3), cmap="jet", cbar=None, ylabel=' ()', xlabel=" []", yscale="log", title='-  \n( )') ticks = ax.set_yticks([2,4,8, 16,32]) ticks = ax.set_yticklabels([2,4,8, 16,32]) show()

Wenn wir das Skalendiagramm mit dem erhaltenen Scan vergleichen, sind sie mit Ausnahme der Farbpalette identisch und zeigen daher in der Zeitreihe die gleiche Dynamik.
Listing 4 ist einfacher als Listing 3 und bietet eine bessere Leistung.

Wenn das Frequenzspektrum und das Leistungsspektrum nach Fourier keine zusätzlichen Informationen über die Dynamik der Zeitreihen zulassen, ist die Verwendung eines Skalogramms vorzuziehen.

Schlussfolgerungen

Ein Beispiel für die Wavelet-Analyse einer Zeitreihe mit Hilfe der PyWavelets-Bibliothek wird gegeben.

Wavelet-Analyse. Teil 2

Einführung

Kontinuierliche Wavelet-Transformation (CWT)

Wavelet - Zeitreihenanalyse mit CWT PyWavelets

Verwendung des Skalogrammmoduls

Schlussfolgerungen

More articles: