April 26, 2019 at 18:02

Wavelet analysis. The basics

  • Tutorial

Introduction


The English word wavelet (from the French "ondelette") literally translates as "short (small) wave." In various translations of foreign articles into Russian, there are also terms: “burst”, “burst function”, “low-wave function”, “wave”, etc.

Wavelet transform (VP) is widely used for signal analysis. In addition, it finds great application in the field of data compression. The VP of a one-dimensional signal is its representation in the form of a generalized series or Fourier integral over a system of basis functions.

$ \ psi _ {ab} (t) = \ frac {1} {\ sqrt {a}} \ psi \ left (\ frac {tb} {a} \ right) $, (1)

constructed from the mother (source) wavelet$ \ psi (t) $possessing certain properties due to time shift operations (b) and temporal scale changes (a).

Factor$ 1 / \ sqrt {a} $ensures the independence of the norm of functions (1) from the scaling number (a). For the given values ​​of the parameters a and b, the function$ \ psi_ {ab} (t) $ and there is a wavelet generated by the mother wavelet $ \ psi (t) $.

An example is the Mexican hat wavelet in the time and frequency domains:

Wavelet listing for time domain
from numpy import*
import matplotlib.pyplot as plt
x= arange(-4,30,0.01)
def w(a,b,t):    
    f =(1/a**0.5)*exp(-0.5*((t-b)/a)**2)* (((t-b)/a)**2-1)
    return f
plt.title("Вейвлет «Мексиканская шляпа»:\n$1/\sqrt{a}*exp(-0,5*t^{2}/a^{2})*(t^{2}-1)$")
y=[w(1,12,t) for t in x]
plt.plot(x,y,label="$\psi(t)$ a=1,b=12") 
y=[w(2,12,t) for t in x]
plt.plot(x,y,label="$\psi_{ab}(t)$ a=2 b=12")   
y=[w(4,12,t) for t in x]
plt.plot(x,y,label="$\psi_{ab}(t)$ a=4 b=12")   
plt.legend(loc='best')
plt.grid(True)
plt.show()




Listing for the wavelet spectrum
from numpy import*
from pylab import *
from scipy import *
import os
def w(a,b,t):    
    f =(1/a**0.5)*exp(-0.5*((t-b)/a)**2)* (((t-b)/a)**2-1)
    return f
x= arange(-4,30,0.2)
def plotSpectrum(y,Fs):
 n = len(y) 
 k = arange(n)
 T = n/Fs
 frq = k/T
 frq = frq[range(int(n/2))] 
 Y = fft(y)/n 
 Y = Y[range(int(n/2))]
 return Y,frq 
 xlabel('f (Hz)')
 ylabel('|Y(f)|')
Fs=1024.0
y=[w(1,12,t) for t in x]
Y,frq=plotSpectrum(y,Fs)
plot(frq,abs(Y),label="$\psi(\omega)$ a=1,b=12")
y=[w(2,12,t) for t in x]
Y,frq=plotSpectrum(y,Fs)
plot(frq,abs(Y),label="$\psi_{ab}(\omega)$ a=2 b=12")
y=[w(4,12,t) for t in x]
Y,frq=plotSpectrum(y,Fs)
plot(frq,abs(Y),label="$\psi_{ab}(\omega)$ a=4 b=12")
plt.title("Вейвлет «Мексиканская шляпа» частотная область $\omega$")
legend(loc='best')
grid(True)
show()




Conclusion:

1. There is indeed a relationship between the concept of Fourier harmonics and the scale of the wavelet. The main thing in this relationship is the inverse proportion of the natural frequency and scale. In addition, reducing the scale, we increase the bandwidth of the wavelet spectrum.

2. By changing the scale (increasing a leads to a narrowing of the Fourier spectrum of the function$ \ psi_ {ab} (t) $), wavelets are able to detect the difference in characteristics on different scales (frequencies), and due to the shift, analyze the signal properties at different points on the entire studied interval. Therefore, when analyzing non-stationary signals, due to the locality of wavelets, they get a significant advantage over the Fourier transform, which gives only global information about the frequencies (scales) of the analyzed signal, since the system of functions used (complex exponent or sines and cosines) is defined on infinite interval.

3. The above listings written in the freely distributed high-level Python language allow you to select functions for wavelets that meet the specified requirements. However, it is additionally necessary to take into account all the main signs of wavelets.

The main signs of the wavelet


Limitedness. The square of the norm of the function must be finite:
$ \ left \ | \ psi \ right \ | ^ {2} = \ int _ {- \ infty} ^ {\ infty} \ left | \ psi (t) \ right | ^ {2} dt <\ infty $. (2)

Localization. VP, in contrast to the Fourier transform, uses a localized initial function both in time and in frequency. To do this, it is enough that the following conditions are met:

$ \ left | \ psi (t) \ right | \ leq C (1+ \ left | t \ right |) ^ {- 1- \ varepsilon} $ and
$ \ left | S _ {\ psi} (\ omega) \ right | \ leq C (1+ \ left | \ omega \ right |) ^ {- 1- \ varepsilon} $ at $ \ varepsilon> 0 $, (3)

For example, a delta function$ \ delta (t) $and harmonic function do not satisfy the necessary condition for simultaneous localization in the time and frequency domains.

Zero average. The graph of the original function should oscillate (be alternating) around zero on the time axis and have zero area:

$ \ int _ {- \ infty} ^ {\ infty} \ psi (t) dt = 0 $. (4)

From this condition, it becomes clear the choice of the name “wavelet” - a small wave.

Equal to zero function area$ \ psi (t) $, i.e. zero moment, leads to the fact that the Fourier transform$ S _ {\ psi} (\ omega) $ this function is equal to zero for $ \ omega = 0 $and has the form of a band-pass filter. For various values ​​(a), this will be a set of bandpass filters.

Often for applications it is necessary that not only zero, but all the first n moments be equal to zero:

$ \ int _ {- \ infty} ^ {\ infty} t ^ {^ {n}} \ psi (t) dt = 0 $. (5)

Waves of the nth order make it possible to analyze the finer (high-frequency) structure of the signal, suppressing its slowly changing components.

Self-made. A characteristic feature of VP is its self-similarity. All wavelets of a specific family$ \ psi_ {ab} (t) $ have the same number of oscillations as the mother wavelet $ \ psi (t) $, since obtained from it by means of scale transformations (a) and shift (b).

Continuous Wavelet Transform


Continuous (integral) wavelet transform (NVP or СWT - continuous wavelet transform). We construct the basis$ \ psi_ {ab} (t) $ using continuous scale transformations (a) and transfers (b) of the mother wavelet $ \ psi (t) $with arbitrary values ​​of the basis parameters a and b in the formula (1).

Then, by definition, the direct (analysis) and reverse (synthesis) NVP (i.e., PNVP and ONVP) of the signal S (t) are written as follows:
$ W_ {s} (a, b) = (S (t), \ psi _ {ab} (t)) = \ frac {1} {\ sqrt {a}} \ int _ {- \ infty} ^ {\ infty} S (t) \ psi \ left (\ frac {tb} {a} \ right) dt $, (6)

$ S (t) = \ frac {1} {C _ {\ psi}} \ int _ {- \ infty} ^ {\ infty} \ int _ {- \ infty} ^ {\ infty} W_ {s} (a, b ) \ psi _ {ab} (t) \ frac {dadb} {a ^ {2}} $, (7)

where$ C _ {\ psi} $ - normalizing coefficient,
$ C _ {\ psi} = \ int _ {- \ infty} ^ {\ infty} \ left | \ psi (\ omega) \ right | ^ {2} \ left | \ omega \ right | ^ {- 1} d \ omega <\ infty $
where: (•, •) is the scalar product of the corresponding factors,
$ \ mathbf {\ psi (\ omega)} $- Fourier transform of a wavelet $ \ psi (t) $. For orthonormal wavelets$ C _ {\ psi} = 1 $.

It follows from (6) that the wavelet spectrum$ W_ {s} (b, a) $(wavelet spectrum, or time-scale-spectrum), unlike the Fourier spectrum (single spectrum), is a function of two arguments: the first argument a (time scale) is similar to the oscillation period, i.e. is inverse to frequency, and the second b –– is similar to the signal offset along the time axis.

It should be noted that$ W_ {s} (b, a_ {0}) $ characterizes the time dependence (at $ a = a_ {0}) $, while the dependencies $ W_ {s} (a, b_ {0}) $ it is possible to match the frequency dependence (for $ b = b_ {0} $)

If the studied signal S (t) is a single pulse of duration$ \ tau_ {u} $concentrated in the neighborhood $ t = t_ {0} $, then its wavelet spectrum will have the greatest value in the vicinity of the point with coordinates $ a = \ tau_ {u}, b = t_ {0} $.

Presentation methods$ W_ {s} (b, a) $may be different. Spectrum$ W_ {s} (b, a) $is a surface in three-dimensional space. However, often instead of the image of the surface, its projection onto the ab plane is presented with iso levels (or figures of various colors), which make it possible to trace the change in the intensity of the EP amplitudes at different scales (a) and in time (b).

In addition, they depict lines of local extrema of these surfaces, the so-called skeleton (sceleton), which reveals the structure of the analyzed signal.

Continuous wavelet transform when determining the wavelet spectrum based on the mother wavelet - “Mexican hat”.



Other mother wavelets used for NVP are also used:



Continuous VP has found wide application in signal processing. In particular, wavelet analysis (VA) gives unique opportunities to recognize local and “subtle” features of signals (functions), which is important in many areas of radio engineering, communications, radio electronics, geophysics and other branches of science and technology. Let us consider this possibility with some simple examples.

Harmonic oscillation.

The signal has the form:$ s (t) = Asin (\ omega t- \ phi) $

Where:$ A = 1 V, \ omega = \ frac {2 \ pi} {T} = \ frac {2 \ pi} {50}, \ psi = 0 $

Wavelet-forming function:$ MHAT (t): = \ frac {d ^ {2}} {dt ^ {2}} exp (-t ^ {2} / 2) $,

Wavelets:$ \ psi (a, b, t) = \ frac {1} {\ sqrt {a}} MHAT \ left (\ frac {tb} {a} \ right) $,

Wavelet spectrum: N: = 256, a: = 1..30, b: = 0..50,$ W (a, b): = \ int _ {- N} ^ {N} \ psi (a, b, t) s (t) dt, N_ {ab}: = W (a, b). $.

Program listing
from scipy.integrate import quad
from numpy import*
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
N=256
T=50
def S(t):
    return sin(2*pi*t/T)
plt.figure()
plt.title(' Гармоническое колебание', size=12)
y=[S(t) for t in arange(0,100,1)]
x=[t for t in arange(0,100,1)]
plt.plot(x,y)
plt.grid()
def w(a,b):    
    f = lambda t :(1/a**0.5)*exp(-0.5*((t-b)/a)**2)* (((t-b)/a)**2-1)*S(t)
    r= quad(f, -N, N)
    return round(r[0],3)
x = arange(1,50,1)
y = arange(1,50,1)
z = array([w(i,j) for j in y for i in x])
X, Y = meshgrid(x, y)
Z = z.reshape(49,49)
fig = plt.figure("Вейвлет- спектр: гармонического колебания")
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet)
ax.set_xlabel(' Масштаб:a')
ax.set_ylabel('Задержка: b')
ax.set_zlabel('Амплитуда ВП: $ N_{ab}$')
plt.figure("2D-график для z = w (a,b)")
plt.title('Плоскость ab с цветовыми областями ВП', size=12)
plt.contourf(X, Y, Z,100)
plt.show()




Signal graph.



Graph of a two-parameter spectrum$ N_ {ab} = W (ab) $displayed as a surface in three-dimensional space.



It should be noted that the section W (a, b) for the time scale$ a = a_ {0} $characterizes the initial oscillation s (t). Moreover, its amplitude is maximum at$ a_ {0}: 1 / \ omega $. Dependencies$ W (a_ {0}, b_ {0}) $ can match the current spectrum of oscillations at $ b = b_ {0} $.

The sum of two harmonic oscillations.

The signal has the form:$ s (t): = A_ {1} sin (\ omega _ {1} t) + A_ {2} sin (\ omega _ {2} t) $

Where: $ A_ {1} = A_ {2} = 1V, \ omega_ {1} = \ frac {2 \ pi} {T_ {1}}, \ omega_ {2} = \ frac {2 \ pi} {T_ {2 }}, T_ {1} = 50s, T_ {2} = 10s $.

$ MHAT (t): = \ frac {d ^ {2}} {dt ^ {2}} \ left [e ^ {- t ^ {2/2}} \ right] $, N: = 256,
$ \ psi (a, b, t): = MHAT \ left (\ frac {tb} {a} \ right), W (a, b): = \ int _ {- N} ^ {N} \ psi (a , b, t) f (t) dt $, a: = 1 ... 30, b: = 0 ... 50, $ N_ {ab}: = W (a, b) $.

Program listing
from scipy.integrate import quad
from numpy import*
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
N=256
def S(t):
    return sin(2*pi*t/10)+sin(2*pi*t/50)
plt.figure(' Сумма двух гармонических колебаний')
plt.title(' Сумма двух гармонических колебаний', size=12)
y=[S(t) for t in arange(0,250,1)]
x=[t for t in arange(0,250,1)]
plt.plot(x,y)
plt.grid()
def w(a,b):    
    f = lambda t :(1/a**0.5)*exp(-0.5*((t-b)/a)**2)* (((t-b)/a)**2-1)*S(t)
    r= quad(f, -N, N)
    return round(r[0],3)
x = arange(1,50,1)
y = arange(1,50,1)
z = array([w(i,j) for j in y for i in x])
X, Y = meshgrid(x, y)
Z = z.reshape(49, 49)
fig = plt.figure("Вейвлет-спектр:2-х гармонических колебаний")
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet)
ax.set_xlabel(' Масштаб:a')
ax.set_ylabel('Задержка: b')
ax.set_zlabel('Амплитуда ВП: $ N_{ab}$')
plt.figure("2D-график для z = w (a,b)")
plt.title('Плоскость ab с цветовыми областями ВП', size=12)
plt.contourf(X, Y, Z, 100)
plt.figure()
q=[w(2,i) for i in y]
p=[i for i in y]
plt.plot(p,q,label='w(2,b)')
q=[w(15,i) for i in y]
plt.plot(p,q,label='w(15,b)')
q=[w(30,i) for i in y]
plt.plot(p,q,label='w(30,b)')
plt.legend(loc='best')
plt.grid(True)
plt.figure()
q=[w(i,13) for i in x]
p=[i for i in x]
plt.plot(p,q,label='w(a,13)')
q=[w(i,17) for i in x]
plt.plot(p,q,label='w(a,17)')
plt.legend(loc='best')
plt.grid(True)
plt.show()





Signal graph.



The plot of the two-parameter spectrum W (a, b) is displayed as a surface in three-dimensional space.



The plane of parameters a, b on which the results of the EP are highlighted in colored areas.



The graph shows the “cross sections” of the wavelet spectrum for two values ​​of the parameter a. With a relatively small parameter of the time scale a, for$ a_ {1} = 2 (a_ {1}: 1 / \ omega_ {2}) $, the cross section of the spectrum carries information only about the high-frequency component of the signal, filtering out (suppressing) its low-frequency component.

As a increases, the extension of the basis function$ \ psi (\ frac {tb} {a}) $, therefore, narrowing its spectrum, and narrowing the bandwidth of the frequency "window". As a result, when$ a_ {2} = 15 (a_ {2}: 1 / \ omega_ {1}) $the cross section of the spectrum is only a low-frequency component of the signal.

With a further increase in a, the window band still decreases and the level of this low-frequency component decreases to a constant component (for a> 25), which is equal to zero for the analyzed signal.



The graph shows the cross sections of the wavelet spectrum W (a, b) characterizing the
current signal spectrum for$ b_ {1} = 13 $ and $ b_ {2} = 17 $. The spectrum of the considered signal, in contrast to the harmonic, contains a high-frequency component in the region of small values ​​of the time scale a (a: 1..3), which corresponds to the second component of the signal$ A_ {2} sin (\ omega_ {2} t) $.

Rectangular momentum.

$ U: = 5, t_ {0}: = 20, \ tau: = 60 $
$s(t):=\begin{vmatrix}U, if t_{0}\leq t\leq t_{0}+\tau \\ 0, otherwise \end{vmatrix}$
$MHAT (t):=\frac{d^{2}}{dt^{2}}exp\left ( \frac{-t^{2}}{2} \right )$
$N:=128,\psi (a,b,t):=MHAT\left ( \frac{t-b}{a} \right ), W(a,b):=\int_{-N}^{N}\psi (a,b,t)f(t)dt$
$a:=1..50, b:=0..100, N_{ba}:=W(a,b)$

Program listing
from scipy.integrate import quad
from numpy import*
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
N=256
def S(t):
    U=5;t0=20;tau=60
    if  t0<=t<=t0+tau:
        return U
    else:
        return 0
plt.figure()
plt.title('Прямоугольный импульс', size=12)
y=[S(t) for t in arange(0,120,1)]
x=[t for t in arange(0,120,1)]
plt.plot(x,y)
plt.grid()
def w(a,b):    
    f = lambda t :(1/a**0.5)*exp(-0.5*((t-b)/a)**2)* (((t-b)/a)**2-1)*S(t)
    r= quad(f, -N, N)
    return round(r[0],3)
x = arange(1,100,1)
y = arange(1,100,1)
z = array([w(i,j) for j in y for i in x])
X, Y = meshgrid(x, y)
Z = z.reshape(99, 99)
fig = plt.figure("3D-график вейвлет спектрограммы")
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet)
ax.set_xlabel(' Масштаб:a')
ax.set_ylabel('Задержка: b')
ax.set_zlabel('Амплитуда ВП: $ N_{ab}$')
plt.figure("2D-график для z = w (a,b)")
plt.title('Плоскость ab с цветовыми областями ВП', size=12)
plt.contourf(X, Y, Z,100)
plt.show()











The wavelet spectra are shown in graphs, the wavelet spectrum well conveys the subtle features of the signal - its jumps on the samples b = 20 and b = 80 (for a: 1..10).

conclusions


This publication is educational in nature, it provides basic information about wavelet analysis in general, and simple examples in the freely distributed high-level programming language Python show the features of continuous wavelet analysis with extensive graphical 3D and 2D visualization.

PS The author does not detract from the unconditional advantages of wavelet analysis using Matlab both in terms of the number of wavelets and speed. But in Python there is still room for further development: scipy.signal and PyWavelets.
Tags:

Also popular now: