usingNeuroAnalyzerusingCairoMakieCairoMakie.activate!(type="png")usingDSPusingFFTWusingStatisticsusingRandomusingWaveletsusingContinuousWavelets# NEURO object with one channel from a signal vector# (used in the examples below)functionsignal2obj(s::Vector{Float64}, fs::Int64) obj =create_object(; data_type ="eeg") obj.header.recording[:sampling_rate] = fsadd_channel!(obj; data =reshape(s, 1, :, 1), label ="s",type="eeg", unit ="μV")create_time!(obj; fs = fs)return objend# Makie figure and axis settings matching NeuroAnalyzer plotsfig_size = (900, 450)ax_style = (titlesize =18, xlabelsize =18, ylabelsize =18, xticklabelsize =12, yticklabelsize =12, xticks =LinearTicks(10), xminorticksvisible =true, xminorticks =IntervalsBetween(10), xautolimitmargin = (0, 0), yautolimitmargin = (0.1, 0.1))# for small panels of multi-panel figuresax_style_small =merge(ax_style, (titlesize =14, xticks =LinearTicks(5), xminorticks =IntervalsBetween(5)))
What are wavelets?
Wavelets are short oscillatory waveforms that start and end at (or decay rapidly to) zero. Unlike the infinitely long sines and cosines of Fourier analysis, wavelets are localized in both time and frequency, which makes them suited to signals whose frequency content changes over time: bursts of oscillations in the EEG, spikes, transients, or sudden changes of state.
Core properties
Zero mean: a wavelet oscillates around zero (its integral is zero), so it responds to changes and oscillations, not to the constant (DC) level of the signal.
Finite energy: a wavelet is concentrated in a short time interval.
Adaptive windows: wavelets are stretched for low frequencies and compressed for high frequencies, so that the analysis window adapts to the frequency.
Fourier transform, STFT and wavelets
Time localization
Frequency resolution
Best suited for
Fourier transform
none (the whole signal)
fine, the same at all frequencies
stationary signals
Short-time FT (STFT)
fixed window length
the same at all frequencies (set by the window)
quasi-stationary signals
Wavelet transform
adaptive: short at high, long at low frequencies
adaptive: fine at low, coarse at high frequencies
non-stationary signals, transients, bursts
In an EEG recording:
Fourier: tells which frequencies are present, but not when.
STFT: shows how the spectrum changes over time, with the same time and frequency resolution at all frequencies; a window long enough to resolve slow rhythms blurs fast events, and vice versa.
Wavelets: use long windows for slow rhythms and short windows for fast activity, so slow oscillations and brief high-frequency events can be analyzed together.
where: \(a > 0\) is the scale: \(a > 1\) stretches the wavelet (lower frequency), \(a < 1\) compresses it (higher frequency), \(b\) is the translation: the wavelet is centered at \(t = b\), \(1 / \sqrt{a}\) keeps the energy of the wavelet the same at all scales.
# real Morlet mother wavelet, 1 Hz at scale a = 1psi(t) =cos(2π * t) *exp(-t^2/ (2*0.8^2))psi_ab(t, a, b) =psi((t - b) / a) /sqrt(a)t =-6:0.01:6fig, ax =lines(t, psi_ab.(t, 1, 0); label ="a = 1, b = 0", figure = (size = fig_size,), axis = (ax_style..., xlabel ="Time [s]", ylabel ="Amplitude"))lines!(ax, t, psi_ab.(t, 2, 0); color =:red, label ="a = 2, b = 0 (stretched)")lines!(ax, t, psi_ab.(t, 0.5, 3); color =:green, label ="a = 0.5, b = 3 (compressed, shifted)")axislegend(ax; position =:lb)fig
Scale and frequency
Scale and frequency are inversely related: doubling the scale halves the frequency of the wavelet. If the mother wavelet has a center frequency\(f_c\) (in cycles per sample), the frequency analyzed at scale \(a\) is:
\[
f = \frac{f_c}{a \, \Delta t}
\]
where \(\Delta t = 1 / f_s\) is the sampling interval. For example, a wavelet with a center frequency of 20 Hz at \(a = 1\) analyzes 10 Hz at \(a = 2\) and 40 Hz at \(a = 0.5\).
Large scales (stretched wavelets): low frequencies, slow changes (e.g. delta waves or baseline drift); good frequency resolution, poor time resolution.
Small scales (compressed wavelets): high frequencies, rapid changes (e.g. gamma activity or spikes); good time resolution, poor frequency resolution.
Dyadic scales
For the discrete wavelet transform, scales and translations are powers of two:
With suitable wavelets, the dyadic wavelets form an orthonormal basis: the signal is represented without redundancy and can be computed with fast filter banks (the fast wavelet transform).
Morlet Wavelet
The Morlet wavelet is the most widely used wavelet for time-frequency analysis of neurophysiological signals. It is a complex sine wave (the carrier) multiplied by a Gaussian (the envelope):
where: \(f\) is the frequency of the wavelet, \(\sigma\) is the standard deviation of the Gaussian (the width of the wavelet in time), \(n\) is the number of cycles, which sets the width of the Gaussian relative to the period of the wavelet.
The real part of the complex Morlet wavelet is a cosine tapered by the Gaussian, the imaginary part is a sine tapered by the Gaussian. Because the wavelet is complex, convolution with it gives both the amplitude and the phase of the signal at frequency \(f\).
generate_morlet(fs, f, t; ncyc, complex) returns a Morlet wavelet with frequency f sampled at fs from -t to t seconds:
With the same number of cycles, the 10 Hz wavelet is 5 times shorter than the 2 Hz wavelet: the time window adapts to the frequency.
The Morlet wavelet acts as a band-pass filter: its spectrum is a Gaussian centered at \(f\) with standard deviation \(\sigma_f = 1 / (2 \pi \sigma) = f / n\).
Number of Cycles
The number of cycles controls the trade-off between time and frequency precision:
Cycles
Wavelet length
Temporal precision
Frequency precision
Use
Few (3–5)
short
high
low (broad band)
transients, short bursts, event-related changes
Many (10–15)
long
low
high (narrow band)
sustained rhythms, separating close frequencies
The precision is usually described by the full width at half maximum (FWHM) of the Gaussian, in time and in frequency:
so that their product is constant: shortening the wavelet in time broadens it in frequency.
10 Hz Morlet wavelets with 3, 7 and 15 cycles, in time and in frequency:
t =-2:1/fs:2nfft =2^14f_axis =fftfreq(nfft, fs)fig =Figure(; size = (900, 600))ax1 =Axis(fig[1, 1]; ax_style_small..., title ="Time domain (real part)", xlabel ="Time [s]", ylabel ="Amplitude")ax2 =Axis(fig[2, 1]; ax_style_small..., title ="Frequency domain", xlabel ="Frequency [Hz]", ylabel ="Normalized amplitude")for (n, color) in ((3, :blue), (7, :red), (15, :green)) mw =generate_morlet(fs, 10, 2; ncyc = n, complex =true)lines!(ax1, t, real.(mw); color = color, label ="$n cycles")# amplitude spectrum, zero-padded for a smooth curve X =abs.(fft(vcat(mw, zeros(nfft -length(mw))))) X ./=maximum(X) idx =sortperm(f_axis)lines!(ax2, f_axis[idx], X[idx]; color = color)endhlines!(ax2, 0.5; color =:black, linestyle =:dot)xlims!(ax1, -1, 1)xlims!(ax2, 0, 25)axislegend(ax1; position =:rt)fig
FWHM measured from the wavelets and calculated from the formulas above:
for n in (3, 7, 15) mw =generate_morlet(fs, 10, 2; ncyc = n, complex =true)# time: width of the envelope at half maximum env =abs.(mw) idx =findall(env .>=0.5*maximum(env)) fwhm_t = t[idx[end]] - t[idx[1]]# frequency: width of the amplitude spectrum at half maximum X =abs.(fft(vcat(mw, zeros(nfft -length(mw))))) idx =findall(X .>=0.5*maximum(X)) fwhm_f =maximum(f_axis[idx]) -minimum(f_axis[idx]) σ = n / (2π *10)println(rpad("$n cycles:", 11),"FWHM time = ", round(1000* fwhm_t; digits =0), " ms ","(theory ", round(1000*2*sqrt(2*log(2)) * σ; digits =0), " ms), ","FWHM frequency = ", round(fwhm_f; digits =1), " Hz ","(theory ", round(2*sqrt(2*log(2)) / (2π * σ); digits =1), " Hz)")end
3 cycles: FWHM time = 109.0 ms (theory 112.0 ms), FWHM frequency = 7.8 Hz (theory 7.8 Hz)
7 cycles: FWHM time = 258.0 ms (theory 262.0 ms), FWHM frequency = 3.3 Hz (theory 3.4 Hz)
15 cycles: FWHM time = 555.0 ms (theory 562.0 ms), FWHM frequency = 1.6 Hz (theory 1.6 Hz)
A 3-cycle wavelet at 10 Hz averages the signal over ≈ 110 ms, but mixes frequencies over an ≈ 8 Hz wide band (e.g. theta and alpha); a 15-cycle wavelet separates frequencies 1.6 Hz apart, but averages over more than half a second.
Number of cycles across frequencies
With a constant number of cycles, the time window is long at low frequencies and short at high frequencies, and the frequency band width grows with frequency. A common compromise is to increase the number of cycles with frequency (e.g. 3 cycles at 2 Hz to 10 cycles at 40 Hz): low frequencies keep a reasonable time resolution and high frequencies a reasonable frequency resolution.
The number of cycles for each frequency can be spaced linearly or logarithmically between the lowest and the highest value:
Logarithmic spacing increases the number of cycles slowly at the low end and faster at the high end; it is the natural choice when the frequencies themselves are spaced logarithmically. In NeuroAnalyzer, Morlet wavelet functions accept ncyc as a single number or as a tuple, e.g. ncyc = (3, 10), in which case the number of cycles is spaced logarithmically from 3 to 10 across frequencies.
Morlet wavelets defined by FWHM
Instead of the number of cycles, the width of the wavelet can be defined directly by its FWHM in time (Cohen, 2019), which is easier to interpret and to report:
where \(h\) is the FWHM in seconds. generate_morlet_fwhm(fs, f, t; h) generates such a wavelet:
t =-1:1/fs:1mw =generate_morlet_fwhm(fs, 10, 1; h =0.25)fig, ax =lines(t, real.(mw); label ="real part", figure = (size = fig_size,), axis = (ax_style..., title ="Morlet wavelet, 10 Hz, FWHM = 250 ms", xlabel ="Time [s]", ylabel ="Amplitude"))lines!(ax, t, abs.(mw); color =:black, linestyle =:dash, label ="envelope")hlines!(ax, 0.5; color =:black, linestyle =:dot)vlines!(ax, [-0.125, 0.125]; color =:red, linestyle =:dot)axislegend(ax; position =:rt)fig
The FWHM should be at least about one cycle of the wavelet frequency (100 ms at 10 Hz); otherwise the wavelet is no longer an oscillation and its frequency band becomes very broad.
Wavelet Convolution
Convolving a signal with a complex Morlet wavelet gives a complex signal whose:
real part is the signal band-pass filtered around the wavelet frequency,
magnitude is the amplitude envelope of that frequency band,
angle is its phase.
A synthetic signal with a 10 Hz burst between 1 and 2 s, a 25 Hz burst between 2.5 and 3 s and noise:
fs =256t =collect(0:1/fs:4-1/fs)Random.seed!(1)x =0.3.*randn(length(t))burst10 =1.<= t .<2burst25 =2.5.<= t .<3x[burst10] .+=sin.(2π *10.* t[burst10])x[burst25] .+=sin.(2π *25.* t[burst25])
Convolution with a 10 Hz Morlet wavelet (5 cycles); the wavelet is normalized so that the envelope is in the units of the signal:
mw =generate_morlet(fs, 10, 1; ncyc =5, complex =true)mw ./=sum(abs.(mw)) /2# "same" convolution: remove half of the wavelet length at both endshalf = (length(mw) -1) ÷2c =conv(x, mw)[(half +1):(half +length(x))]fig, ax =lines(t, x; color =:gray, label ="signal", figure = (size = fig_size,), axis = (ax_style..., xlabel ="Time [s]", ylabel ="Amplitude"))lines!(ax, t, real.(c); color =:blue, label ="real part (10 Hz)")lines!(ax, t, abs.(c); color =:red, linewidth =2, label ="magnitude (10 Hz envelope)")axislegend(ax; position =:rt)fig
The envelope rises during the 10 Hz burst; the 25 Hz burst is outside the pass band of the wavelet. The envelope starts rising before 1 s and decays after 2 s: this temporal smoothing is set by the length of the wavelet (number of cycles).
Repeating the convolution for many frequencies gives a time-frequency representation (a scalogram when expressed by scale): power (squared magnitude) as a function of time and frequency.
Time-Frequency Analysis
The same signal analyzed with the short-time Fourier transform (1-second window), Morlet wavelet convolution and the continuous wavelet transform:
obj =signal2obj(x, fs)
plot_spectrogram(obj; ch ="s", method =:stft, flim = (0, 40), db =false, title ="STFT")
plot_spectrogram(obj; ch ="s", method =:mw, ncyc = (3, 10), flim = (0, 40), db =false, title ="Morlet wavelets, 3 to 10 cycles")
plot_spectrogram(obj; ch ="s", method =:cwt, flim = (0, 40), db =false, title ="Continuous wavelet transform")
STFT: the same frequency resolution (≈ 1 Hz) for both bursts, and the same temporal smearing at both frequencies; the onsets and offsets are blurred by the window, and the first and last half-windows (0–0.5 s, 3.5–4 s) are missing, because the time axis gives the window centers.
Morlet wavelets: the onsets and offsets are sharper, especially for the 25 Hz burst (short wavelets), but the frequency bands are broader, especially at 25 Hz.
CWT: similar to Morlet wavelet convolution, with fewer, logarithmically spaced frequencies (scales).
Many finely spaced scales at every time point; highly redundant.
Time-frequency analysis, visualization.
Discrete (DWT)
Dyadic scales; at each level the signal is split into approximation (low-pass) and detail (high-pass) and downsampled.
Compression, denoising, multiresolution analysis.
Stationary (SWT)
DWT without downsampling (the filters are upsampled instead); shift-invariant and redundant.
Denoising, feature extraction.
Wavelet packet transform
Both the approximation and the detail are decomposed at every level, giving uniform frequency bands.
Extracting specific frequency bands.
Multiresolution analysis
In the DWT, the signal is passed through a pair of filters: a low-pass filter gives the approximation (A) and a high-pass filter gives the detail (D). Only the approximation is decomposed further, so each level halves the frequency band:
At 256 Hz, D1 covers 64–128 Hz, D2 32–64 Hz, D3 16–32 Hz and A3 0–16 Hz. Because the bands are split in halves, EEG bands such as alpha (8–13 Hz) do not correspond to single DWT bands; the wavelet packet transform also splits the details and gives finer, uniform bands.
Common wavelets
Wavelet
Characteristics
Morlet
Complex sine wave with a Gaussian envelope; amplitude and phase; time-frequency analysis (CWT).
Mexican hat (Ricker)
Negative second derivative of a Gaussian; real-valued; detection of peaks and edges (CWT).
Paul
Complex; good time localization (CWT).
Haar
The simplest wavelet: a step function; compact support; detection of abrupt changes (DWT).
Daubechies (dbN)
Orthogonal, compact support, N vanishing moments; smoother with increasing N (DWT).
Symlets (symN)
Nearly symmetric modifications of Daubechies wavelets (DWT).
Coiflets
Vanishing moments for both the wavelet and the scaling function (DWT).
Biorthogonal
Separate analysis and synthesis filters; symmetric, linear-phase filters (DWT, e.g. JPEG 2000).
The choice depends on the application: Morlet wavelets for oscillations (amplitude and phase), Mexican hat or Haar wavelets for transients and edges, Daubechies or Symlets wavelets for denoising and compression.
Discrete Wavelet Decomposition
dwd() performs discrete wavelet decomposition. Two types are available, both without downsampling (all coefficient vectors have the length of the signal):
The decomposition is a wavelet packet tree: both the low-pass and the high-pass branches are decomposed. The default wavelet is Haar (wt = wavelet(WT.haar)); other wavelets from Wavelets.jl can be used, e.g. wt = wavelet(WT.db4). The number of levels is set with l (by default the maximum for the signal length).
The coefficients are returned as an array of channels × nodes × samples × epochs. The nodes (rows) are ordered level by level:
row 1: the signal (level 0),
rows 2–3: level 1 (low-pass, high-pass),
rows 4–7: level 2,
rows 8–15: level 3,
and the children of node \(k\) are nodes \(2k\) (low-pass) and \(2k + 1\) (high-pass). At level 3, each of the 8 nodes covers 1/8 of the frequency range (16 Hz at 256 Hz sampling rate), but the nodes are not in frequency order: after a high-pass step, the order of the following low- and high-pass bands is reversed. The 4 Hz component (0–16 Hz band) is in node 8 and the 40 Hz component (32–48 Hz band) in node 11:
node 8 standard deviation: 1.902
node 9 standard deviation: 0.218
node 10 standard deviation: 0.081
node 11 standard deviation: 0.459
node 12 standard deviation: 0.011
node 13 standard deviation: 0.059
node 14 standard deviation: 0.023
node 15 standard deviation: 0.131
Plotting the coefficients of a 2-level decomposition (the plot numbers the nodes without the signal: coefficient #1 is node 2, i.e. the level 1 low-pass node):
dc2 =dwd(obj; ch ="s", type=:acdwt, l =2)plot_dwc(dc2[1, :, :, 1]; t = obj.epoch_time)
Inverse discrete wavelet decomposition
idwd() reconstructs the signal from the nodes of the last level. Without changes, the reconstruction is exact:
maximum reconstruction error: 1.1102230246251565e-15
Selecting nodes with c reconstructs only the corresponding frequency bands; the other nodes are set to zero. Reconstructing the 4 Hz component from node 8:
x_4hz =idwd(d; type=:acdwt, c =8)fig, ax =lines(t, x; color =:gray, label ="signal", figure = (size = fig_size,), axis = (ax_style..., xlabel ="Time [s]", ylabel ="Amplitude"))lines!(ax, t, x_4hz; color =:red, linewidth =2, label ="reconstructed from node 8")xlims!(ax, 0, 1)axislegend(ax; position =:rt)fig
The Haar wavelet has a poor frequency selectivity, so neighboring bands leak into each other; smoother wavelets (e.g. WT.db4) separate the bands better.
Continuous Wavelet Decomposition
cwd() performs the continuous wavelet transform using ContinuousWavelets.jl. The default wavelet is wavelet(Morlet(2π), β = 2); it can be changed with wt, e.g. wt = wavelet(dog2) (Mexican hat) or wt = wavelet(paul2). The parameter β sets the spacing of the scales: smaller values give more scales per octave.
A signal whose frequency changes from 2 Hz to 16 Hz at 2 s:
plot_spectrogram(obj; ch ="s", method =:cwt, flim = (1, 30), db =false, title ="Continuous wavelet transform")
The scalogram shows the low frequency before 2 s and the high frequency after it; at 2 s the transient is localized sharply by the short, high-frequency wavelets and broadly by the long, low-frequency wavelets.
icwd() reconstructs the signal from the coefficients. With complex wavelets such as Morlet, the complex coefficients must be used (cwd(...; complex = true)); the real part alone does not contain the full information and the reconstruction is poor:
# middle of the signal (edges excluded)m =129:(length(x) -128)x_rec =icwd(ct[1, :, :, 1])println("relative error, complex coefficients: ",round(std(x_rec[m] .- x[m]) /std(x[m]); digits =3))ct_real =cwd(obj; ch ="s")x_rec =icwd(ct_real[1, :, :, 1])println("relative error, real coefficients: ",round(std(x_rec[m] .- x[m]) /std(x[m]); digits =3))
relative error, complex coefficients: 0.061
relative error, real coefficients: 0.359
The continuous wavelet frame is redundant and the inverse transform is approximate; the default reconstruction method (type = :pd, Penrose delta dual frame) is usually the most accurate. The accuracy depends on the wavelet: more scales per octave (smaller β, e.g. wt = wavelet(Morlet(2π), β = 1)) give a more accurate reconstruction; denoise_cwd() and artrem_cwd() use β = 1 by default. Frequencies above the highest wavelet frequency (≈ 87 Hz at 256 Hz sampling rate) are not represented and are lost in the reconstruction.
Wavelet Denoising
Denoising with discrete wavelet decomposition
DWT-based denoising (wavelet shrinkage) relies on the sparsity of the wavelet representation: the signal is concentrated in a few large coefficients, while white noise is spread over all coefficients with small amplitudes. Setting small coefficients to zero (thresholding) and reconstructing the signal removes most of the noise:
hard thresholding: coefficients below the threshold are set to zero, the others are kept,
soft thresholding: coefficients below the threshold are set to zero, the others are shrunk towards zero by the threshold.
Common thresholds include the universal threshold \(\sigma \sqrt{2 \ln n}\) (where \(\sigma\) is the noise standard deviation and \(n\) the number of samples), SURE (Stein’s unbiased risk estimate) and minimax thresholds.
l is the number of decomposition levels, and smooth the smoothing method: :regular thresholds all coefficients, :undersmooth does not threshold the lowest-frequency node. The denoising method is set with dnt (default: RelErrorShrink(SoftTH())); see WaveletsExt.jl for the available methods.
Denoising with continuous wavelet decomposition
denoise_cwd() removes a frequency band: the coefficients between nf - w and nf + w Hz are reconstructed by the inverse transform and subtracted from the signal. Because the inverse CWT is approximate and the wavelets overlap in frequency, a narrow band is removed only partially. A 10 Hz signal with 50 Hz noise, and the amplitude of the 50 Hz noise remaining for different band widths:
x_clean =sin.(2π *10.* t)x = x_clean .+0.5.*sin.(2π *50.* t)obj =signal2obj(x, fs)# amplitude of a sinusoid of frequency f (least squares, edges excluded)m =129:(length(t) -128)functionsine_amp(s, f) A = [sin.(2π * f .* t[m]) cos.(2π * f .* t[m])]returnhypot((A \ s[m])...)endfor w in (5, 10, 20) obj_dn =denoise_cwd(obj; ch ="s", nf =50, w = w) x_dn = obj_dn.data[1, :, 1]println("w = ", lpad(w, 2), " Hz: 50 Hz amplitude ",round(sine_amp(x_dn, 50); digits =2), " (before: 0.5), ","10 Hz amplitude ", round(sine_amp(x_dn, 10); digits =2)," (before: 1.0)")end
With a narrow band (w = 5) the noise is not removed (its amplitude even increases slightly); it is removed only with a wide band (here 30–70 Hz), which also removes all other activity in that band. For narrow-band noise such as power line interference, notch filters are more precise; see Power line noise.
Removing Artifacts
Because wavelet coefficients are localized in time and frequency, an artifact that occupies a limited time-frequency region can be removed by setting the coefficients in that region to zero and reconstructing the signal. artrem_cwd() does this for one channel and one epoch; the time segment tseg is given in seconds of epoch time and the frequency segment fseg in Hz.
To check the result, a synthetic eye blink (a slow, large deflection) is added to 10 seconds of a real EEG channel:
plot_spectrogram(obj; ch ="s", method =:cwt, flim = (0, 20), title ="EEG with a blink")
The blink occupies ≈ 3.5–4.5 s and frequencies below ≈ 5 Hz. The coefficients in this time-frequency region are reconstructed by the inverse transform and subtracted from the signal, so the rest of the signal is not changed:
x_clean = obj_clean.data[1, :, 1]fig, ax =lines(t, x_eeg .+ blink; color =:gray, label ="EEG with blink", figure = (size = fig_size,), axis = (ax_style..., xlabel ="Time [s]", ylabel ="Amplitude [μV]"))lines!(ax, t, x_eeg; color =:black, label ="original EEG")lines!(ax, t, x_clean; color =:red, label ="after artrem_cwd()")xlims!(ax, 2, 6)axislegend(ax; position =:rt)fig
win =3.4.<= t .<=4.6out = .!win .& (1.<= t .<=9)println("difference from the original EEG, 3.4–4.6 s: ",round(100*std(x_clean[win] .- x_eeg[win]) /std(blink[win]); digits =1), "% of the blink")println("difference outside the segment: ",round(100*std(x_clean[out] .- x_eeg[out]) /std(x_eeg[out]); digits =1), "% of the EEG")
difference from the original EEG, 3.4–4.6 s: 33.3% of the blink
difference outside the segment: 0.2% of the EEG
About two thirds of the blink is removed, and the signal outside the segment is practically unchanged. The EEG activity in the same time-frequency region (delta and theta waves) is removed with the blink, and the removal is incomplete because the inverse transform is approximate. The time and frequency segments should be as narrow as possible, and the result should always be inspected.
Practical Considerations
Frequency limits
The highest frequency must be below the Nyquist frequency (\(f_s / 2\)).
The lowest frequency is limited by the signal length: the wavelet must fit into the signal. A Morlet wavelet with \(n\) cycles at frequency \(f\) lasts several times \(n / f\) seconds; e.g. 5 cycles at 1 Hz require several seconds of data.
Frequencies closer to each other than the frequency resolution of the wavelets do not give independent information (e.g. 14.9 and 15 Hz); 20–40 logarithmically spaced frequencies are usually enough for EEG.
Edge effects
Near the beginning and the end of the signal, the wavelet extends beyond the data and the coefficients are distorted; the affected region (the cone of influence) is longer for low frequencies. Analyze long continuous segments or epochs with buffer periods before and after the time window of interest, and discard the edges.
Power and phase
The squared magnitude of the complex coefficients gives the power, and their angle the phase. Power is often expressed in dB or as a change relative to a baseline period, because power decreases with frequency (1/f).
Advantages and limitations
Advantages
Limitations
Time-frequency localization adapted to the frequency.
Time and frequency resolution cannot both be high (uncertainty principle).
Amplitude and phase from complex wavelets.
Results depend on the choice of the wavelet and its parameters (e.g. cycles).
Suited to non-stationary signals, bursts and transients.
CWT is redundant and computationally expensive; its inverse is approximate.
Sparse DWT representation for denoising and compression.
DWT bands are fixed by the dyadic scheme and may not match physiological bands.
Example workflow
Choose a wavelet, e.g. complex Morlet for oscillations, Haar or Daubechies for denoising.
Choose the frequencies and the number of cycles, e.g. logarithmically spaced frequencies and 3–10 cycles.
Compute the transform on long segments and discard the edges.
Extract features: power, phase, band power over time, or denoised and artifact-free signals.
References
Cohen MX. Analyzing Neural Time Series Data: Theory and Practice. MIT Press; 2014.
Cohen MX. A better way to define and describe Morlet wavelets for time-frequency analysis. NeuroImage. 2019;199:81–6.
Mallat S. A Wavelet Tour of Signal Processing. 3rd ed. Academic Press; 2009.
Torrence C, Compo GP. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society. 1998;79(1):61–78.
