```python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
```

## Convolution

Convolution is often employed to obtain the degree of similarities between analyzed signal and prototype functions.
It's major advantage is to use known atoms to investigate unknown phenomena.
It is an important concept for Fourier transform, Wavelet transform etc.

The convolution of two functions $f(t)$ and $g(t)$ is:
\begin{equation}
h(t) = f(t)* g(t)= \int_{-\infty}^{\infty}f(\tau)g(t-\tau)d\tau
\end{equation}

For discrete time series, the convolution is a sum instead of an integral.
\begin{equation}
h(n) = \sum_{m=-\infty}^{\infty}f(m)g(n-m)
\end{equation}

The most common fast convolution algorithms use fast Fourier transform (FFT) algorithms via the circular convolution theorem.


```python
def conv(f, g, same=True):
    """
    discrete convolution:
        i: [0, len(f)+len(g)-1)
        j: [max(0, i-len(g)+1), min(i, len(f)-1)]
    """
    len_f = len(f)
    len_g = len(g)
    len_h = len_f + len_g - 1
    h = np.zeros(len_h)
    for i in range(len_h):
        for j in range(max(0, i-len_g+1), min(i, len_f-1)+1):
            h[i] += f[j]*g[i-j]
    if same:
        l = (len(g)-1)//2
        h = h[l:len(h)-(len(g)-l-1)]
    return h
```


```python
def gaussian(t, sigma, mu):
    return np.e ** (-0.5*((t-mu)/sigma)**2)

def Plot(x,y,z):
    # plot the morlet atom
    plt.figure(figsize=(10,5))
    ax = plt.gca()
    ax.grid(color='#b7b7b7', linestyle='-', linewidth=0.5, alpha=0.5)
    plt.plot(t,np.real(x), color='#333333', linestyle=':', alpha=0.5)
    plt.plot(t,np.real(y), color='#333333', linestyle=':', alpha=0.5)
    plt.plot(t,np.real(z), color='#121212', linestyle='-')
    plt.plot(t,np.imag(z), color='#121212', linestyle='-.')

# sampling rate
Fs = 100
# sampling interval
Ts = 1.0/Fs
# time vector
t = np.arange(0, 10, Ts)
# frequency of the atom
w = 2*np.pi*1
x=np.e**(-1j*w*t)
# construct the morlet atom
y=gaussian(t, 1, 5)
morlet=x*y    
Plot(x,y,morlet)

```

### Morlet wavelet (standard version)
\begin{equation}
\psi(t) = \pi^{-1/4}e^{jwt}e^{-t^2/2}
\end{equation}


### Daughter wavelet atoms
\begin{equation}
\psi_{a,b}(t) = \frac1{\sqrt{a}}\psi\left(\frac{t - b}{a}\right)
\end{equation}


```python
def morlet(N, s=1.0, w=6.0):
    """
    Complex Morlet wavelet.
    
    :param N: (int) Length of the wavelet
    :param s: (float) Scaling factor. Default is 1.0
    :param w: (float) Omega0. Default is 5.0
    :return: morlet: (ndarray) Morlet wavelet
    """
    t = np.linspace(-1/s*2*np.pi, 1/s*2*np.pi, N)
    return np.pi**(-0.25)*np.exp(1j*w*t)*np.exp(-0.5*(t**2))

def fourier(x):
    """
    Fourier transform.
    
    :param x: (ndarray) input signal
    :return: X: (ndarray) one side frequency of input signal
    """
    n = len(x)         # length pf the signal
    X = np.fft.fft(x)  # fast fourier transform
    X /= n             # normalization
    X = X[range(n//2)] # one side frequency range
    return abs(X)

def Plot(morlet_atoms, scales, morlet_freqs, freq):
    alphas = [0.4, 0.7, 1]
    plt.figure(figsize=(10,5))
    plt.subplot(2,1,1)
    ax = plt.gca()
    ax.grid(color='#b7b7b7', linestyle='-', linewidth=0.5, alpha=0.5)
    for i in range(len(morlet_atoms)):
        plt.plot(np.real(morlet_atoms[i]), color="#333333", alpha=alphas[i], label=str(scales[i]))
    plt.legend(ncol=4, loc='upper right', fontsize='small')
    plt.subplot(2,1,2)
    ax = plt.gca()
    ax.grid(color='#b7b7b7', linestyle='-', linewidth=0.5, alpha=0.5)
    for i in range(len(morlet_atoms)):
        plt.plot(freq, morlet_freqs[i], color="#333333", alpha=alphas[i])
    
fs = 1
dt=1/fs
n = 256
T = n * dt
freq = np.arange(n)/T    # two side frequencies
freq = freq[range(n//2)] # one side frequencies

scales = [0.5, 1, 2]
morlet_atoms = [1/np.sqrt(s) * morlet(n, s) for s in scales]
morlet_freqs = [fourier(x) for x in morlet_atoms]
Plot(morlet_atoms, scales, morlet_freqs, freq)
```

### Continuous wavelet transform
\begin{equation}
W_f(a,b) = \frac1{\sqrt{a}}\int_{-\infty}^{+\infty}f(t)\psi^*\left(\frac{t - b}{a}\right)dt
\end{equation}
It is possible to calculate the wavelet transform using convolution in time space, it is considerably faster to do the calculations in Fourier space.


```python
t = np.linspace(-1, 1, 256, endpoint=False)
data  = 2*np.sin(2 * np.pi * 1 * t) + np.sin(2 * np.pi * 5 * t)*(t>0)
scales = np.arange(1, 21)

wavelet_spectrum = np.zeros([len(scales), len(data)], dtype=np.complex)
for ind, scale in enumerate(scales):
    wavelets = morlet(min(10 * scale, len(data)), scale)
    wavelet_spectrum[ind, :] = conv(dasta, wavelets)

wavelet_power_spectrum = abs(wavelet_spectrum)**2

fig = plt.figure(figsize=(10, 5))
fig.subplots_adjust(left=0.09, bottom=0.09, right=0.95, top=0.95,
                    hspace=0.05, wspace=0.05)
plt.subplot(2,1,1)
ax = plt.gca()
ax.grid(color='#b7b7b7', linestyle='-', linewidth=0.5, alpha=0.5)
plt.setp(ax.get_xticklabels(), visible=False)
plt.plot(data, '#333333')
plt.subplot(2,1,2)
plt.imshow(wavelet_spectrum.real, cmap='PRGn', aspect='auto',
           vmax=abs(wavelet_spectrum.real).max(), vmin=-abs(wavelet_spectrum.real).max())
plt.show()

```

65e13a5f5d6b0fe2dadd550dbc41f27137cd4b71
ipynb
Jupyter Notebook
004_wavelet_transform.ipynb
ValleyZw/hill
51d4bc76d7c49b3dccf88659fade806906a9b4e2
MIT
2019-07-22T07:19:39.000Z
2019-07-22T07:19:39.000Z
004_wavelet_transform.ipynb
ValleyZw/hill
51d4bc76d7c49b3dccf88659fade806906a9b4e2
MIT
004_wavelet_transform.ipynb
ValleyZw/hill
51d4bc76d7c49b3dccf88659fade806906a9b4e2
MIT
2020-05-19T17:52:16.000Z
2020-05-19T17:52:16.000Z
Qwen/Qwen-72B
1. YES
2. YES