😐 👴🏽 💯 Analisis wavelet. Bagian 2 🔇 🧑🏾‍🤝‍🧑🏽 👇🏼

Pendahuluan

Publikasi ini membahas analisis wavelet dari deret waktu. Gagasan utama transformasi wavelet sesuai dengan spesifikasi banyak deret waktu yang menunjukkan evolusi dalam waktu karakteristik utama mereka - nilai rata-rata, varian, periode, amplitudo, dan fase komponen harmonik. Sebagian besar proses yang dipelajari dalam berbagai bidang pengetahuan memiliki fitur di atas.

Tujuan dari publikasi ini adalah untuk menggambarkan metode transformasi wavelet berkelanjutan dari deret waktu menggunakan pustaka PyWavelets ..

Sedikit sejarah

Ahli Geofisika D. Morlaix pada akhir 70-an abad XX. Saya mengalami masalah menganalisis sinyal dari sensor seismik yang mengandung komponen frekuensi tinggi (aktivitas seismik) untuk periode waktu yang singkat dan komponen frekuensi rendah (keadaan tenang kerak bumi) untuk jangka waktu yang lama. Transformasi Window Fourier memungkinkan Anda untuk menganalisis komponen frekuensi tinggi atau komponen frekuensi rendah, tetapi tidak kedua komponen sekaligus.

Oleh karena itu, metode analisis diusulkan di mana lebar fungsi jendela meningkat untuk frekuensi rendah dan menurun untuk frekuensi tinggi. Transformasi jendela baru diperoleh sebagai hasil dari peregangan (kompresi) dan penggantian waktu dari satu fungsi pembangkit (fungsi penskalaan - fungsi penskalaan, scalet). Fungsi pembangkit ini disebut wavelet oleh D. Morlaix.

Wavelet D. Morlaix

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

Transformasi Wavelet Berkelanjutan (CWT)

Transformasi wavelet tingkat-tunggal:

coefs, frequency = pywt.cwt (data, skala, wavelet)

dimana:

data: (seperti array_) Sinyal input ;

skala (seperti array_) : - Skala wavelet untuk digunakan. Anda dapat menggunakan rasio f = scale2frequency (scale, wavelet) / sampling_ period dan menentukan berapa frekuensi fisik f. Di sini, f di hertz ketika sampling_ periode dalam detik;

wavelet (objek atau nama wavelet ) : - Wavelet untuk digunakan;

sampling_ period (float): - Periode sampling untuk frekuensi output (parameter opsional). Nilai yang dihitung untuk koefisien tidak bergantung pada sampling_ periode (yaitu, skala tidak diskalakan dengan periode sampling.);

coefs (seperti array_) : - Wavelet kontinu - transformasi sinyal input untuk skala dan wavelet yang diberikan;

frekuensi (seperti array_) : - Jika unit periode sampling adalah kedua dan diatur, maka frekuensi ditampilkan dalam hertz. Kalau tidak, periode pengambilan sampel 1 diasumsikan.

Catatan: Ukuran array koefisien tergantung pada panjang array input dan panjang skala yang diberikan.

Dapatkan daftar nama wavelet yang kompatibel dengan cwt yang tersedia:

 >>> import pywt >>> wavlist = pywt.wavelist(kind='continuous') >>> wavlist

['cgau1', 'cgau2', 'cgau3', 'cgau4', 'cgau5', 'cgau6', 'cgau7', 'cgau8', 'cgau8', 'cmor', 'fbsp', 'gaus1', 'gaus2', ' gaus3 ',' gaus4 ',' gaus5 ',' gaus6 ',' gaus7 ',' gaus8 ',' mexh ',' morl ',' morl ',' shan ']

Untuk analisis wavelet, pilihan wavelet ibu adalah salah satu tugas utama. Oleh karena itu, sebelum setiap pemilihan wavelet, Anda hanya perlu menggunakan program berikut yang memungkinkan Anda untuk menentukan properti dinamis wavelet:

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()

Selanjutnya, kami akan mempertimbangkan fungsi wavelet dan propertinya menggunakan program di atas:

Topi Meksiko wavelet mexh:

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

Wavelet "morl" Morlaix:

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

Wavelet Morlet Kompleks "cmorB-C"

$\ varphi \ kiri (t \ kanan) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi Ct}$

di mana B adalah bandwidth; C adalah frekuensi pusat.

Selanjutnya (tanpa indikasi khusus) nilai-nilai B, C ditetapkan dengan titik mengambang.

Wavelet Gaussian “gausP”

$\ varphi \ kiri (t \ kanan) = C \ cdot exp ^ {- t ^ {2}}$

di mana: P - bilangan bulat dari 1 hingga 8 dimasukkan ke dalam nama wavelet, misalnya, "gaus8"; C-built-in normalisasi konstan.

Wavelet Gaussian Terpadu “cgauP”

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

di mana: P - bilangan bulat dari 1 hingga 8 dimasukkan ke dalam nama wavelet, misalnya, "gaus8" dan sesuai dengan urutan turunan dari fungsi wavelet; C-built-in normalisasi konstan.

Gelombang air Shannon "shanB-C"

$\ varphi \ kiri (t \ kanan) = \ sqrt {B} \ cdot \ frac {sin (\ pi Bt)} {\ pi B T} \ cdot exp ^ {j 2 pi Ct}$

di mana: B adalah bandwidth; C adalah frekuensi pusat.

CWT dalam PyWavelets berlaku untuk data diskrit - konvolusi dengan sampel integral wavelet. Jika skalanya terlalu kecil, maka kita mendapatkan filter diskrit dengan pengambilan sampel yang tidak memadai, yang mengarah ke perataan, seperti yang ditunjukkan pada grafik untuk wavelet 'cmor1.5-1.0'.

Di kolom kiri, grafik menunjukkan filter diskrit yang digunakan dalam konvolusi untuk berbagai skala. Kolom kanan adalah spektrum daya Fourier yang sesuai dari setiap filter. Dengan skala 1 dan 2 untuk grafik 'cmor1.5-1.0' dapat dilihat bahwa smoothing terjadi karena pelanggaran terhadap batasan Nyquist. Fenomena yang ditunjukkan dapat terjadi pada wavelet lain ketika memilih C dan B, oleh karena itu, ketika bekerja dengan wavelet, Anda harus menggunakan program - Listing 1.

Untuk menghubungkan skala tertentu dengan frekuensi sinyal tertentu, periode pengambilan sampel sinyal harus diketahui. Menggunakan fungsi pywt.scale2frequency () , Anda bisa mengonversi daftar skala ke frekuensi yang sesuai. Pilihan skala yang tepat tergantung pada wavelet yang dipilih, jadi pywt.scale2frequency () harus digunakan untuk mendapatkan gambaran rentang frekuensi sinyal yang sesuai.

 >>> 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.])

Analisis Wavelet - Time Series Menggunakan CWT PyWavelets

Dataset el-Nino adalah dataset deret waktu yang digunakan untuk melacak El Nino dan berisi pengukuran suhu permukaan laut triwulanan dari tahun 1871 hingga 1997. Untuk memahami kekuatan grafik skala, mari kita visualisasikan untuk dataset el-Nino bersama dengan data deret waktu asli dan transformasi Fourier-nya.

Untuk analisis wavelet dari deret waktu, langkah-langkah berikut harus dilakukan:

1. Pilih wavelet ibu: Pilih wavelet Morlet kompleks "cmorB-C":

$\ varphi \ kiri (t \ kanan) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi Ct}$

Bandwidth - B dan frekuensi tengah - Dari mana kita akan menentukan pada langkah selanjutnya.

2. Tentukan frekuensi pusat dengan mengambil dt = 0,25 untuk deret waktu kami:

 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. Kami menemukan transformasi Fourier dari mother wavelet cmor1.0-0.5 (kami menggunakan daftar 1):
Wavelet berkelanjutan akan dievaluasi pada seluruh rentang [-8.0, 8.0]

4. Pilih data untuk deret waktu:

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

Data tersebut adalah pengukuran suhu permukaan laut triwulanan dari tahun 1871 hingga 1997. Untuk data ini: t0 = 1871 dt = 0.25

5. Kami membangun deret waktu dengan sinyal dan moving average pada satu grafik:

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. Kami melakukan transformasi Fourier dan modalitas spektrum dari seri waktu:

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. Kami menentukan skala: skala = arange (1, 128); level = [2 ** - 4, 2 ** - 3, 2 ** - 2, 2 ** - 1, 2 ** 0, 2 ** 1, 2 ** 2, 2 ** 3].

8. Buat grafik skala:

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)

Grafik skala menunjukkan bahwa sebagian besar kekuatan spektrum terkonsentrasi selama 2-8 tahun, ini sesuai dengan frekuensi 0,125 - 0,5 Hz. Dalam spektrum Fourier, fashion juga berkonsentrasi pada nilai-nilai frekuensi ini. Namun, transformasi wavelet juga memberi kita informasi sementara, tetapi transformasi Fourier tidak.

Sebagai contoh, pada diagram skala, kita melihat bahwa sebelum 1920 ada banyak fluktuasi, sedangkan antara 1960 dan 1990 tidak ada banyak. Kita juga dapat melihat bahwa dari waktu ke waktu ada pergeseran dari periode yang lebih pendek ke periode yang lebih lama.

Menggunakan modul scaleogram

Scaleogram adalah alat yang mudah digunakan untuk menganalisis data 1D dengan transformasi wavelet kontinu berdasarkan pustaka PyWavelets. Modul ini dirancang untuk menyediakan alat yang andal untuk analisis data cepat.

Scaleogram memiliki beberapa fitur berikut:

tanda tangan panggilan pemula sederhana
sumbu dapat dibaca dan integrasi matplotlib murni
banyak pilihan untuk zoom, filter spektrum, implementasi colorbar, dll ...
dukungan untuk unit frekuensi dan frekuensi sesuai dengan tanda
kecepatan, menggunakan algoritma N * log (N) untuk transformasi
portabilitas: diuji dengan python2.7 dan python3.7
rincian pesan kesalahan dan dokumentasi dengan contoh-contoh

Namun, dalam contoh penggunaan, akses data tidak direkam dengan benar:

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

Peringatan berikut muncul:

nino3 = pd.read_table (nino3)
mengarah ke peringatan:
Peringatan (dari modul peringatan):
File “C: / Pengguna / Pengguna / Desktop /wavelet /pr.py”, baris 7
nino3 = pd.read_table (nino3)
FutureWarning: read_table sudah usang, gunakan read_csv sebagai gantinya, lewat sep = '\ t'

Akses data harus ditulis seperti ini:

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

Untuk menunjukkan penggunaan skalogram dan membandingkan hasilnya dengan hasil dari contoh sebelumnya, kami menggunakan data yang sama - pada pengukuran triwulanan suhu permukaan laut dari tahun 1871 hingga 1997. Untuk data ini: 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()

Jika kita membandingkan grafik skala dengan hasil pemindaian, maka, dengan pengecualian palet warna, mereka identik, dan karenanya menunjukkan dinamika yang sama dalam deret waktu.
Listing 4 lebih sederhana daripada listing 3 dan memiliki kinerja yang lebih baik.

Ketika spektrum frekuensi dan spektrum daya menurut Fourier tidak memungkinkan untuk memperoleh informasi tambahan tentang dinamika deret waktu, maka penggunaan skalogram lebih disukai.

Kesimpulan

Contoh analisis wavelet dari deret waktu dengan menggunakan perpustakaan PyWavelets diberikan.

Analisis wavelet. Bagian 2

Pendahuluan

Transformasi Wavelet Berkelanjutan (CWT)

Analisis Wavelet - Time Series Menggunakan CWT PyWavelets

Menggunakan modul scaleogram

Kesimpulan

More articles: