Перед вами обучающий материал по основам цифровой обработки сигналов с использованием средств языка программирования Python. Предполагается, что читатель имеет базовые знания из области высшей математики, а также владеет языком Python и хотя бы поверхностно знает различные python-библиотеки - numpy/scipy, matplotlib и другие.
Для пользователей MATLAB / GNU Octave освоение материала с точки зрения программного кода не составит труда, поскольку основные функции и их атрибуты во многом идентичны и схожи с методами из python-библиотек.
Вейвлет-преобразование — интегральное преобразование, представляющее собой свертку вейвлет-функции с сигналом. Вейвлет-преобразование похоже на оконное преобразование Фурье, так как переводит сигнал из временного представления в частотно-временное.
Wavelet дословно переводится как "короткая (маленькая) волна". Вейвлет это волнообразная функция с быстро затухающей амплитудой и нулевым средним.
Вейвлеты преобразует функцию в двумерную поверхность:
$$ y(t) \to T(t, f) \tag{15.1}$$где $T(t, f)$ - вклад частоты $f$ в момент времени $t$ в сигнале.
Чтобы произвольная функция $\Phi(t)$ могла считаться вейвлет-функцией, она должна удовлетворять двум условиям:
то есть не иметь компоненты нулевой частоты.
то есть быть локализованной во времени, не бесконечной.
Рассмотрим и визуализируем несколько популярных вейвлет-функций.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
В python есть библиотека pyWavelets, позволяющая легко вычислять вейвлет-преобразования. Библиотека содержит реализации множества вейвлет-функций из различных семейств, как дискретных, так и непрерывных.
import pywt
# All wavelet families available in pywt
pywt.families(short=False)
['Haar', 'Daubechies', 'Symlets', 'Coiflets', 'Biorthogonal', 'Reverse biorthogonal', 'Discrete Meyer (FIR Approximation)', 'Gaussian', 'Mexican hat wavelet', 'Morlet wavelet', 'Complex Gaussian wavelets', 'Shannon wavelets', 'Frequency B-Spline wavelets', 'Complex Morlet wavelets']
# Haar wavelet
_, haar, _ = pywt.Wavelet('haar').wavefun()
# Gaussian wavelet
gauss, _ = pywt.ContinuousWavelet('gaus8').wavefun()
# Morlet wavelet
morlet, _ = pywt.ContinuousWavelet('morl').wavefun()
# Mexican hat
mex_hat, _ = pywt.ContinuousWavelet('mexh').wavefun()
# Plot wavelet functions
wavelets = [haar, gauss, morlet, mex_hat]
titles = ['Haar', 'Gaussian', 'Morlet', 'Mexican hat']
fig = plt.figure(figsize=(16, 10), dpi=120)
for i, sig in enumerate(wavelets):
plt.subplot(2, 2, i+1)
plt.title(titles[i], fontsize=18)
plt.plot(sig)
plt.ylabel('Amplitude', fontsize=16)
plt.xlabel('Time', fontsize=16)
plt.grid()
plt.tight_layout()
Так как морлет является комплексной функцией и имеет как вещественную, так и мнимую часть, для наглядности рассмотрим его в 3D.
from scipy import signal
# График интерактивный
%matplotlib notebook
fs = 500
n_cycles = 5
freq = 5
scaling = 1.0
wavelet_len = int(n_cycles * fs / freq)
morlet = signal.morlet(wavelet_len, n_cycles, scaling)
# Plot real and imaginary parts in a 3D plot
fig = plt.figure(dpi=120)
ax = fig.add_subplot(111, projection='3d')
ax.plot(np.linspace(0, scaling, morlet.size), morlet.real, morlet.imag)
plt.title('Morlet in 3D')
ax.set(xlabel='Scaling', ylabel='Real Amplitude', zlabel='Imag Amplitude')
plt.show()
Вейвлеты, используемые в вейвлет-преобразовании, получаются с помощью сдвига и масштабирования базовой вейвлет-функции (материнского вейвлета) наподобие рассмотренных выше.
$a$ - коэффициент масштабирования.
Как видно из формулы, шкала вейвлета обратно пропорционально его частоте. Растянутые вейвлеты (с низкой частотой) помогают улавливать медленные изменения в сигнале; сжатые вейвлеты (с высокой частотой), наоборот - быстрые изменения. Ниже показан один вейвлет с разным масштабом (частотой).
%matplotlib inline
w = pywt.Wavelet('sym3')
(phi, psi, x) = w.wavefun(level=3)
(phi2, psi2, x2) = w.wavefun(level=5)
phi = np.append(phi, np.zeros(120))
f, (ax1, ax2) = plt.subplots(1, 2, sharex=True, figsize=(16, 6), dpi=120)
ax1.plot(phi)
ax1.grid()
ax1.set_xlabel('Time', fontsize=16)
ax1.set_ylabel('Amplitude', fontsize=16)
ax2.plot(phi2)
ax2.grid()
ax2.set_xlabel('Time', fontsize=16)
ax2.set_ylabel('Amplitude', fontsize=16)
plt.show()
Видно, что с уменьшением частоты вейвлет становится более растянутым по времени, и наоборот.
С изменением параметра сдвига $b$ базовый вейвлет сдвигается по оси $X$, позволяя захватить сигнал на всей его протяженности. Ниже приведена иллюстрация сдвига вейвлета вдоль сигнала.
from ipywidgets import *
import numpy as np
import matplotlib.pyplot as plt
# График интерактивный
%matplotlib notebook
ns = 1024
time = np.arange(ns)
morlet, smh = pywt.ContinuousWavelet('morl').wavefun()
x = np.linspace(0, 5 * np.pi, 256)
fig = plt.figure(dpi=120)
plt.title('Wavelet moving along the signal')
plt.grid()
line, = plt.plot(morlet)
line2, = plt.plot(np.cos(4 * x))
def update(w = 1.0):
line.set_xdata(range(round(20 * w), round(20 * w) + len(morlet)))
fig.canvas.draw_idle()
interact(update)
plt.show()
interactive(children=(FloatSlider(value=1.0, description='w', max=3.0, min=-1.0), Output()), _dom_classes=('wi…
Таким образом, вейвлеты, получаемые преобразование материнского вейвлета, выглядят следующим образом: $$\Phi_{a, b}(t) = \Phi(\frac{t - b}{a}) \tag{15.10}$$ Меняя параметры сдвига и масштаба, мы можем проанализировать, какие частоты присутствуют в сигнале в определенный момент времени.
Вейвлет-преобразование задается как свертка сигнала и материнского вейвлета $\Phi_{a, b}$ с различными параметрами масштаба $a$ и сдвига $b$. По сути мы в каждый момент времени $t$ оценивает кросс-корреляцию сигнала и вейвлета определенной частоты; полученное значение будет максимальным по модулю там, где частоты вейвлета и сигнала совпадают.
Есть два вида вейвлет-преобразований: непрерывное и дискретное. Они различаются тем, что для непрерывного преобразования масштаб $a$ и сдвиг $b$ являются непрерывными, следовательно, есть бесконечное количество вейвлетов с разным масштабом и сдвигом. Дискретное преобразование, наоборот, использует конечное множество вейвлетов, полученных из конечного множества значений сдвига и масштаба.
где $\hat{\Phi}_{a, b}(t)$ - функция комплексной переменной.
Вейвлет-преобразлвание $T(a, b)$ показывает, насколько материнский вейвлет $\Phi_{a, b}(t)$ с частотой $\frac{t}{a}$ в определенный момент времени $t - b$ похож на исходный сигнал. Когда частота вейвлета близка к частоте сигнала, они совпадают.
Амплитуда у этой осциллирующей функции будет самой большой там, где частоты вейвлета и сигнала полностью совпадают. Вещественная часть - вклад определенной частоты в определенный момент времени.
Рассмотрим, как работает непрерывное вейвлет-преобразование, на примере сложного сигнала, состоящего из наложения колебаний различной частоты. Для этого визуализируем скейлограмму вейвлет-преобразования.
Скейлограмма показывает результат вейвлет-преобразования над сигналом. По оси $X$ на скейлограмме откладывается время, а по оси $Y$ - шкала (то есть масштаб) вейвлета.
# Установка библиотеки scaleogram для работы со скейлограммами
# Раскомментируйте при первом запуске
#! git clone http://github.com/alsauve/scaleogram
#%cd scaleogram
#! python ./setup.py install --user
#%cd ..
import librosa
from IPython.display import Audio
# Example file
audio, sr = librosa.load(librosa.util.example('trumpet'))
print('Trumpet solo')
Audio(audio, rate=sr)
Trumpet solo
%matplotlib inline
import scaleogram as scg
# Range of scales to perform the transform
scales = scg.periods2scales(np.arange(1, 60))
# Plot the original signal
fig1, ax1 = plt.subplots(1, 1, figsize=(15, 6), dpi=120)
lines = ax1.plot(audio)
#ax1.set_xlim(0, sr)
ax1.grid()
ax1.set_title("Original signal")
fig1.tight_layout()
# Compute and plot the scaleogram
ax2 = scg.cws(audio, scales=scales, figsize=(16, 6))
Как видно из графика, вейвлет-преобразование позволяет получить информацию о том, какие частоты присутствуют в сигнале в каждый момент времени.
Хотя вейвлет-преобразование и позволяет получать информацию о частотах и времени их появления в сигнале одновременно, его возможности ограничены фундаментальным принципом неопределенности Гейзенберга. Согласно этому принципу, чем точнее измеряется одна характеристика, тем менее точно мы можем измерить вторую.
Однако, в отличие от оконного преобразования Фурье, у которого разрешение по времени и частоте остается постоянным (при фиксированном размере окна), у вейвлет-преобразований разрешение зависит от шкалы (то есть частоты) вейвлета.
Ниже приведены графики uncertainty boxes для сигнала во временной области и трех видов преобразований: Фурье, оконного преобразования Фурье и вейвлет-преобразования. Грани прямоугольников показывают, каким количеством информации об одной характеристике мы жертвуем, чтобы получить информацию о другой характеристике. Видно ключевое отличие вейвлет-преобразования: разрешение не является фиксированным и зависит от частоты.
from IPython.display import Image
Image(filename='../img/Comparisonoftransformations.jpeg')
Во временной области у нас есть вся информация о времени, однако полностью отсутствует информация о частоте; в частотной области - наоборот.
В оконном преобразовании Фурье разрешение по времени и по частоте остается постоянным, а в случае с вейвлет-преобразованием разрешение меняется: на низких частотах у нас высокое разрешение по частоте, и низкое по времени, на высоких частотах - наоборот.
Это логично, так как низкие частоты длятся продолжительное количество времени, поэтому точное разрешение по времени не так важно, а вот небольшой сдвиг по частоте (2 Гц вместо 1) может значительно изменить ситуацию. И наборот, высокие частоты локализованы на маленьком отрезке времени, поэтому нам важно высокое разрешение по времени, а конкретное значение частоты (500 или 501 Гц) становится менее существенно.
Так как изображения являются 2D-данными, мы можем отдельно посчитать одномерное вейвлет-преобразование для строк и для столбцов. Вейвлет-преобразования позволяют разделить информацию, представленную на картинке, на две части: аппроксимации и детали (подсигналы).
Сигнал проходит через 2 фильтра: фильтр верхних частот и нижних частот. Полученные компоненты высокой частоты это детали, а низкой - аппроксимации.
Сначала фильтры высоких и низких частот применяются к изображению построчно, выходные сигналы из этих фильтров даунсемплятся. Затем к полученным сигналам по столбцам применяются фильтры верхних и. нижних частот и снова производится даунсемплинг.
Image(filename='../img/DWT-decomposition.png')
Тогда мы получим 4 набора коэффициентов: аппроксимация, горизонтальные детали изображения, вертикальные детали изображения и диагональные детали изображения.
Рассмотрим преобразование изображений на примере.
import cv2
image = cv2.imread('../img/kitten.jpeg')
plt.imshow(image)
plt.title('Original image')
plt.show()
# Wavelet transform of image, and plot approximation and details
titles = ['Approximation', ' Horizontal detail',
'Vertical detail', 'Diagonal detail']
coeffs2 = pywt.dwt2(image[:, :, 0], 'bior1.3')
LL, (LH, HL, HH) = coeffs2
fig = plt.figure(figsize=(16, 13), dpi=120)
for i, a in enumerate([LL, LH, HL, HH]):
ax = fig.add_subplot(2, 2, i + 1)
ax.imshow(a, interpolation="nearest", cmap=plt.cm.gray)
ax.set_title(titles[i], fontsize=18)
fig.tight_layout()
plt.show()
Можно делать несколько последовательных вейвлет-преобразований подряд. Тогда на каждом уровне последующем уровне раскладываться будет только аппроксимация с предыдущего уровня. Рассмотрим на примере трехуровневое вейвлет-преобразование.
from pywt._doc_utils import wavedec2_keys, draw_2d_wp_basis
x = pywt.data.camera().astype(np.float32)
shape = x.shape
max_lev = 3 # how many levels of decomposition to draw
label_levels = 3 # how many levels to explicitly label on the plots
fig, axes = plt.subplots(2, 4, figsize=[16, 8], dpi=120)
for level in range(0, max_lev + 1):
if level == 0:
# show the original image before decomposition
axes[0, 0].set_axis_off()
axes[1, 0].imshow(x, cmap=plt.cm.gray)
axes[1, 0].set_title('Image', fontsize=16)
axes[1, 0].set_axis_off()
continue
# plot subband boundaries of a standard DWT basis
draw_2d_wp_basis(shape, wavedec2_keys(level), ax=axes[0, level],
label_levels=label_levels)
axes[0, level].set_title('{} level\ndecomposition'.format(level), fontsize=16)
# compute the 2D DWT
c = pywt.wavedec2(x, 'db2', mode='periodization', level=level)
# normalize each coefficient array independently for better visibility
c[0] /= np.abs(c[0]).max()
for detail_level in range(level):
c[detail_level + 1] = [d/np.abs(d).max() for d in c[detail_level + 1]]
# show the normalized coefficients
arr, slices = pywt.coeffs_to_array(c)
axes[1, level].imshow(arr, cmap=plt.cm.gray)
axes[1, level].set_title('Coefficients\n({} level)'.format(level), fontsize=16)
axes[1, level].set_axis_off()
plt.tight_layout()
plt.show()
Источник: pyWavelets Documentation
Так как каждая функция, удовлетворяющая условиям (14.2), (14.3) является вейвлетом, можно создавать свои функции для вейвлет-преобразований.
В библиотеке pyWavelets для создания своей вейвлет-функции нужно создать объект класса Wavelet
, определить для него название name
и используемый банк фильтров filter_bank
. Фильтры нужно задавать в следующем порядке:
# Construct custom wavelet
const = 10
c = np.sqrt(2)/2
dec_lo, dec_hi, rec_lo, rec_hi = [-c*const, c], [c, c], [c, -c*const], [c, c]
filter_bank = [dec_lo, dec_hi, rec_lo, rec_hi]
custom_wavelet = pywt.Wavelet(name="custom", filter_bank=filter_bank)
# Apply wavelet filters
coeffs = pywt.dwt2(image[:, :, 0], custom_wavelet)
cA, (cH, cV, cD) = coeffs
# Plot result
fig = plt.figure(figsize=(16, 13), dpi=120)
for i, a in enumerate([cA, cH, cV, cD]):
ax = fig.add_subplot(2, 2, i + 1)
ax.imshow(a, interpolation="nearest", cmap=plt.cm.gray)
ax.set_xticks([])
ax.set_yticks([])
fig.tight_layout()
plt.show()