Basic Concept: Wavelet Transform

The wavelet transform uses a group of short kernels (wavelets) of different frequencies convolved with the original signal. Thus, the signal only needs to be stable over the length of the wavelet (stationary).

Similar to Fourier transform, when the signal frequency matches the wavelet frequency, the dot product yields larger amplitude results.

Wavelet convolution can be utilized as a band-pass filter.

Convolution Kernels

Types of Wavelets

There are various types of convolution kernels (wavelets) in wavelet transforms. Different types affect temporal precision. The most commonly used type is the Gaussian kernel (Morlet Wavelet).

Types of Wavelets: MathWorks - Introduction to the Wavelet Families

Wavelet Families

Definition: A set of wavelets of the same type but different frequencies is referred to as a family.

Morlet Wavelet Family

The Morlet wavelet appears as a Gaussian envelope in the time domain with a sinusoidal function inside. In the frequency domain, it resembles a Gaussian function with its peak at the wavelet’s frequency, f.

Real-valued Morlet Wavelet

Mathematical formula:

asin⁡2πft∗ae−t2/(2s2)a\sin 2\pi ft * a e^{-t^2/(2s^2)}

Q: How is the Morlet wavelet constructed?

Construct a sine wave and multiply it by a Gaussian window to obtain the Morlet wavelet. It’s essential to consider the sampling frequency and the number of points in both the Gaussian window and the sine wave.

Q: What is the specific formula?

Sine wave formula: asin⁡2πfta\sin 2\pi ft

Gaussian formula: ae−t2/(2s2)a e^{-t^2/(2s^2)} where s=n2πfs=\frac{n}{2\pi f}

Wavelet formula: asin⁡2πft∗ae−t2/(2s2)a\sin 2\pi ft * a e^{-t^2/(2s^2)}

n represents wavelet cycles, balancing temporal precision and frequency precision, and f denotes the wavelet frequency.

Q: How should one choose frequency and cycles?

Minimum frequency: The theoretical lower limit is 1 times the epoch-related frequency, but typically, 4 times the epoch-related frequency is chosen. For instance, for an epoch of 1 second, at least 4Hz is selected as the wavelet frequency.

Maximum frequency: The theoretical upper limit is the Nyquist frequency, but it’s common to ensure four sample points in each cycle. For example, with a 500Hz sampling frequency, a maximum of 125Hz is chosen as the wavelet frequency.

Number of wavelet families: Generally, 20-30 different frequencies are selected. For EEG cognitive tasks, 4-60Hz is usually chosen.

Q: Are there any issues with real-valued wavelets?

Similar to Fourier transform, pure real-valued convolution kernels require the original signal and the kernel to have the same phase. Hence, complex-valued wavelets are introduced.

Code

Exercise

Suppose I have a signal that is 60 seconds long with a sampling frequency of 100Hz, and I want to perform wavelet convolution. How should I write the code?

Maltab

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
% wavelet parameters
fs = 100;
wTime = -1:(1/fs):0.99;
wFreq = 1:2:25;
n = 4;
wLength = length(wTime);

% creat test data
dataTime = -10:(1/fs):10;
data = sin(2*pi*4*dataTime);
data = [zeros(1,wLength/2) data zeros(1,wLength/2)];

% initialize matrix
wavelet_conv = zeros(length(wFreq),length(dataTime));

for i = 1:length(wFreq)
  f = wFreq(i);
  s = n/(2*pi*f);
  wavelet = cos(2*pi*f*wTime).*exp(-(wTime.^2)/(2*s^2));
  kernel = fliplr(wavelet);
  conv_result = conv(data,kernel,'same')/sum(kernel);
  wavelet_conv(i,:) = conv_result(wLength/2:wLength/2+length(dataTime)-1);
end

% plot
figure;
imagesc(dataTime, wFreq, wavelet_conv)
colormap("jet")
colorbar
xlabel('time')
ylabel('frequency')

Note that there is a side effect in the result, and there is a certain variance at 4Hz

Python

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
import numpy as np
import matplotlib.pyplot as plt

# wavelet parameters
fs = 100
wTime = np.arange(-1, 1, 1/fs)
wFreq = np.arange(1, 26, 2)
n = 4
wLength = len(wTime)

# create test data
dataTime = np.arange(-10, 10, 1/fs)
data = np.sin(2*np.pi*4*dataTime)
data = np.concatenate((np.zeros(wLength//2), data, np.zeros(wLength//2)))

# initialize matrix
wavelet_conv = np.zeros((len(wFreq), len(dataTime)))

for i in range(len(wFreq)):
    f = wFreq[i]
    s = n/(2*np.pi*f)
    wavelet = np.cos(2*np.pi*f*wTime) * np.exp(-(wTime**2)/(2*s**2))
    kernel = np.flip(wavelet)
    conv_result = np.convolve(data, kernel, 'same')/np.sum(kernel)
    wavelet_conv[i, :] = conv_result[wLength//2:wLength//2+len(dataTime)]

# plot
plt.figure()
plt.imshow(wavelet_conv, aspect='auto', extent=[min(dataTime), max(dataTime),  max(wFreq), min(wFreq) ], cmap='jet')
# 注意这里max(wFreq)和min(wFreq)需要对换
plt.colorbar()
plt.xlabel('time')
plt.ylabel('frequency')
plt.show()

Complex-valued Morlet Wavelet

Mathematical formula:

Ae−t2/2s2ei2πftwhereA=1(sπ)1/2Ae^{-t^2/2s^2}e^{i2\pi ft} \quad \text{where} A=\frac{1}{(s\sqrt{\pi})^{1/2}}

The left part is Gaussian, and the right part is a complex sinusoid.

Q: What is a complex wavelet?

The real and imaginary parts of the complex wavelet and their 3D representation are depicted below:

Real and imaginary parts of a complex wavelet

Q: How to interpret the results of convolution?

The result of complex wavelet convolution is a vector in the complex plane.

Convolution result in the complex plane

The magnitude of the complex vector represents how the signal’s power at the wavelet frequency changes over time (unaffected by phase).

The phase of the complex vector represents the phase relationship between the wavelet and the signal.

The real part of the complex vector is the result of the signal passing through a band-pass filter.

Real part, magnitude, and phase angle of the convolution result

Analyzing phase from the image

Q: How to choose parameters for the wavelet?

Number of cycles n: n determines the width of the wavelet. A larger n provides better frequency precision but poorer temporal precision.

Generally, for studying transient changes, 3-4 cycles are preferred, while for studying frequency changes, 7-10 cycles are chosen. Different cycle numbers can also be used at different frequencies.

Choosing the number of cycles

Y-axis scaling:

Using an exponential logarithmic axis highlights the low-frequency range.

Using a linear axis highlights the high-frequency range.

Y-axis scaling

Number of frequencies:

More frequencies yield higher resolution results. However, after a certain point, adding more frequencies may not provide significant benefits.

Number of frequencies

Resulting resolution with different frequencies

Length of wavelet: There is no maximum limit on the wavelet length. Its ends must approach zero, with its center at t=0. The sampling rate should match the signal’s.

Code

! practise

Suppose I have a 60s long signal with a sampling frequency of 100Hz and I want to perform complex-based wavelet convolution. How should the code be written?

Maltab

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
% wavelet parameters
fs = 100;
wTime = -1:(1/fs):0.99;
wFreq = 1:2:25;
n = 4;
wLength = length(wTime);

% creat test data
dataTime = -10:(1/fs):10;
data = sin(2*pi*4*dataTime);
data = [zeros(1,wLength/2) data zeros(1,wLength/2)];

% initialize matrix
wavelet_conv = zeros(length(wFreq),length(dataTime));

for i = 1:length(wFreq)
  f = wFreq(i);
  s = n/(2*pi*f);
  A = 1/((s*sqrt(pi))^0.5);
  wavelet = A.*exp(-(wTime.^2)/(2*s^2)).*exp(1i*2*pi*f.*wTime);
%   wavelet = exp(2*1i*pi*frex(fi).*time) .* exp(-time.^2./(2*(numcycles(fi)/(2*pi*frex(fi)))^2));
  kernel = fliplr(wavelet);
  conv_result = conv(data,kernel,'same');
  wavelet_conv(i,:) = conv_result(wLength/2:wLength/2+length(dataTime)-1);

end

% plot
figure;
imagesc(dataTime, wFreq, abs(wavelet_conv))
colormap("jet")
colorbar
xlabel('time')
ylabel('frequency')

Python

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import numpy as np
import matplotlib.pyplot as plt

# wavelet parameters
fs = 100
wTime = np.arange(-1, 1.01, 1/fs)
wFreq = np.arange(1, 26, 2)
n = 4
wLength = len(wTime)

# create test data
dataTime = np.arange(-10, 10.01, 1/fs)
data = np.sin(2*np.pi*4*dataTime)
data = np.concatenate((np.zeros(wLength//2), data, np.zeros(wLength//2)))

# initialize matrix
wavelet_conv = np.zeros((len(wFreq), len(dataTime)))

for i in range(len(wFreq)):
    f = wFreq[i]
    s = n/(2*np.pi*f)
    A = 1/np.sqrt(s*np.sqrt(np.pi))
    wavelet = A * np.exp(-(wTime**2)/(2*s**2)) * np.exp(1j*2*np.pi*f*wTime)
    kernel = np.flip(wavelet)
    conv_result = np.convolve(data, kernel, 'same')
    wavelet_conv[i, :] = conv_result[wLength//2:wLength//2+len(dataTime)]

# plot
plt.figure()
plt.imshow(np.abs(wavelet_conv), aspect='auto', extent=[min(dataTime), max(dataTime), max(wFreq), min(wFreq)], cmap='jet')
plt.colorbar()
plt.xlabel('time')
plt.ylabel('frequency')
plt.show()

Subsequent Analysis

Using real(wavelet_conv) = Result projected onto the real axis = Passed through a band-pass filter.

Using abs(wavelet_conv) = Magnitude of the vector = Power of the data at the wavelet frequency.

Using angle(wavelet_conv) = Angle of the vector = Phase angle of the data relative to the wavelet at a specific time point.

Efficiency tip: FFT → Multiplication → IFFT can replace convolution (conv) for improved computational efficiency.

Note: The content in this blog is class notes shared for educational purposes only. Some images and content are sourced from textbooks, teacher materials, and the internet. If there is any infringement, please contact aursus.blog@gmail.com for removal.