Вейвлет – анализ. Часть 2

Введение

В данной публикации рассматривается вейвлет – анализ временных рядов. Основная идея вейвлет-преобразования отвечает специфике многих временных рядов, демонстрирующих эволюцию во времени своих основных характеристик – среднего значения, дисперсии, периодов, амплитуд и фаз гармонических компонент. Подавляющее большинство процессов, изучаемых в различных областях знаний, имеют вышеперечисленные особенности.

Целью настоящей публикации является описание методики непрерывного вейвлет- преобразования временных рядов средствами библиотеки PyWavelets..

Немного истории

Инженер-геофизик Д. Морле в конце 70-х годов XX в. столкнулся с проблемой анализа сигналов от сейсмодатчиков, которые содержали высокочастотную компоненту (сейсмическая активность) в течение короткого промежутка времени и низкочастотные составляющие (спокойное состояние земной коры) – в течение длительного периода. Оконное преобразование Фурье позволяет анализировать либо высокочастотную составляющую, либо низкочастотную составляющую, но не обе составляющие сразу.

Поэтому, был предложен метод анализа, в котором ширина оконной функции для низких частот увеличивалась, а для высоких частот – уменьшалась. Новое оконное преобразование получалось в результате растяжения (сжатия) и смещения по времени одной порождающей (так называемой скейлинг-функции – scaling function, scalet) функции. Эта порождающая функция была названа вейвлетом Д. Морле.

Вейвлет Д. Морле

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


Непрерывное вейвлет-преобразование (CWT)

Одноуровневое вейвлет – преобразование:
coefs, frequencies =pywt.cwt(data, scales, wavelet)

где:

data : (array_like) -Входной сигнал;

scales (array_like):- Масштаб вейвлета для использования. Можно использовать соотношение f =scale2frequency(scale, wavelet)/sampling_period и определить какая физическая частота f. Здесь, f в hertz когда sampling_period в секундах;

wavelet (Wavelet object or name): — Вейвлет для использования;

sampling_period (float):- Период дискретизации для выходных частот (необязательный параметр). Значения, вычисленные для coefs, не зависят от выбора sampling_period (т. е. scales не масштабируется выборкой периода.);

coefs (array_like):- Непрерывное вейвлет — преобразование входного сигнала для заданных масштабов и вейвлет;

frequencies (array_like ):- Если единица периода выборки секунды и задана, то частоты выводятся в hertz. В противном случае предполагается период выборки, равный 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’, ‘morl’, ‘shan’]

Для вейвлет — анализа выбор материнского вейвлета является одной из ключевых задач. Поэтому, перед каждым выбором вейвлета просто необходимо пользоваться следующей программой, позволяющей определить динамические свойства вейвлета:

Листинг 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()

Далее будем рассматривать вейвлет-функции и их свойства с использованием приведенной программы:

Мексиканская шляпа вейвлет «mexh»:

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

Вейвлет «morl»Морле:

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

Комплексный вейвлет Морле «cmorB-C»

$varphileft ( t right )= frac{1}{sqrt{pi}B}exp^{-frac{t^{2}}{B}}cdot exp^{j2pi C t}$

где B — пропускная способность; C — центральная частота.
Здесь и далее (без специального указания) величины B,C задаются с плавающей точкой.

Гауссовы вейвлеты «gausP»

$varphileft ( t right )=Ccdot exp^{-t^{2}}$

где: P — целое число от 1 до 8 вставляется в имя вейвлета, например «gaus8»; C- встроенная константа нормализации.

Комплексные гауссовы вейвлеты «cgauP»

$varphileft ( t right )=Ccdot exp^{-t^{2}}cdot exp^{-jt}$

где: P — целое число от 1 до 8 вставляется в имя вейвлета, например «сgaus8» и соответствуют порядку производной от вейвлет функции; C- встроенная константа нормализации.

Вейвлеты Шеннона «shanB-C»

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

где: B — ширина полосы; C — центральная частота.

CWT в PyWavelets применяется к дискретным данным — свертки с образцами интеграла вейвлета. Если scale слишком мало, то мы получаем дискретный фильтр с неадекватной выборкой, приводящий к сглаживанию, как показано на графике для вейвлета ‘cmor1.5-1.0’.

В левой колонке на графиках показаны дискретные фильтры, используемые в свертке при различных шкалах. Правый столбец — соответствующие спектры мощности Фурье каждого фильтра. При шкалах 1 и 2 для графика ‘cmor1.5-1.0’ видно, что сглаживание происходит из-за нарушения ограничение Найквиста. Указанное явление может возникнуть и у других вейвлетов при выборе С и B, поэтому при работе с вейвлетами необходимо пользоваться программой – листинг 1.

Чтобы связать заданный масштаб с определенной частотой сигнала, период дискретизации сигнала должен быть известен. При помощи функции pywt.scale2frequency() можно преобразовать список масштабов в соответствующие частоты. Правильный выбор масштабов зависит от выбранного вейвлета, поэтому pywt.scale2frequency() следует использовать для получения представления о соответствующем частотном диапазоне сигнала.

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

Вейвлет — анализ временных рядов средствами CWT PyWavelets

Набор данных el-Nino представляет собой набор данных временных рядов, используемый для отслеживания El Nino и содержит ежеквартальные измерения температуры морской поверхности с 1871 по 1997 год. Чтобы понять силу масштабограммы, давайте визуализируем ее для набора данных el-Nino вместе с исходными данными временных рядов и его преобразованием Фурье.

Для вейвлет-анализа временных рядов необходимо выполнить следующие действия по пунктам:
1. Выбрать материнский вейвлет: Выбираем комплексный вейвлет Морле «cmorB-C»:

$varphileft ( t right )= frac{1}{sqrt{pi}B}exp^{-frac{t^{2}}{B}}cdot exp^{j2pi C t}$

Пропускную способность – B и центральную частоту – С которого будем определять на следующем этапе.
2. Определить центральную частоту, приняв dt=0.25 для нашего временного ряда:

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. Находим преобразование Фурье материнского вейвлета 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. Строим временной ряд с сигналом и его скользящее среднее на одном графике:

Листинг 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. Проводим преобразование Фурье и модности спектра от временного ряда:

Листинг 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. Определяем масштабы: scales = arange(1, 128); levels = [2**-4, 2**-3, 2**-2, 2**-1, 2**0, 2**1, 2**2, 2**3].
8. Строим масштабограмму:

Листинг 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)

На масштабограмме видно, что большая часть мощности спектра сконцентрирована за период 2-8 лет, это соответствует частоте 0,125 – 0,5 Гц. В спектре Фурье модность так же концентрируется вокруг этих значений частоты. Однако вейвлет — преобразование также дает нам временную информацию, а преобразование Фурье — нет.

Например, на масштабной диаграмме мы видим, что до 1920 года было много колебаний, в то время как между 1960 и 1990 годами их было не так много. Мы также можем видеть, что с течением времени происходит сдвиг от более коротких периодов к более длинным.

Использование модуля scaleogram

Scaleogram-удобный инструмент для анализа 1D данных с непрерывным вейвлет — преобразованием построен на базе библиотеки PyWavelets. Модуль призван обеспечить надежный инструмент для быстрого анализа данных:
Scaleogram имеет следующие особенности:
• простая подпись вызова для начинающих
• читаемые оси и чистая интеграция matplotlib
• много вариантов для изменения масштаба, фильтра спектра, внедрения colorbar, etc…
• поддержка единиц периодичности и частоты в соответствии с маркировкой
• скорость, использует алгоритм N * log(N) для преобразований
• портативность: испытано с python2.7 и python3.7
• подробные сообщения об ошибках и документация с примерами
Однако, в примерах использования не правильно записано обращение к данным:

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

При этом появляется следующее предупреждение:

nino3 = pd.read_table(nino3)
что приводит к предупреждению:
Warning (from warnings module):
File «C:/Users/User/Desktop/wavelet/pr.py», line 7
nino3 = pd.read_table(nino3)
FutureWarning: read_table is deprecated, use read_csv instead, passing 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

Листинг 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()


Если сравнивать масштабограмму с полученной сканограммой, то, за исключением цветовой палитры, они идентичны, а следовательно показывают одинаковую динамику временного ряда.
Листинг 4 проще листинга 3 и имеет большее быстродействие.

Когда частотный спектр и спектр мощности по Фурье не позволяют получить дополнительную информацию о динамике временного ряда, тогда использование скалограммы является предпочтительным.

Выводы:

Приведен пример вейвлет-анализа временного ряда средствами библиотеки PyWavelets

 
Источник

Читайте также