using NeuroAnalyzer
using CairoMakie
CairoMakie.activate!(type = "png")
using DSP
using FFTW
using Statistics
# NEURO object with one channel from a signal vector (used in the examples below)
function signal2obj(s::Vector{Float64}, fs::Int64)
obj = create_object(; data_type = "eeg")
obj.header.recording[:sampling_rate] = fs
add_channel!(obj; data = reshape(s, 1, :, 1), label = "s", type = "eeg", unit = "μV")
create_time!(obj; fs = fs)
return obj
endThe Hilbert transform is used to compute the analytic signal of a real-valued signal. The analytic signal gives the instantaneous amplitude, phase and frequency of the signal at every time point, which makes it a basic tool for the analysis of brain oscillations (envelopes, phase synchronization, cross-frequency coupling).
The Hilbert transform \(H\{x(t)\}\) of a real signal shifts the phase of every frequency component by −90° (−π/2) without changing its amplitude: a cosine becomes a sine and a sine becomes a negative cosine.
The analytic signal \(z(t)\) is a complex signal whose real part is the original signal and whose imaginary part is its Hilbert transform:
\[ z(t) = x(t) + i H\{x(t)\} \]
For a cosine, the analytic signal is a complex sine wave (Euler’s formula):
\[ x(t) = A \cos(2\pi f t + \theta) \Rightarrow z(t) = A \cos(2\pi f t + \theta) + i A \sin(2\pi f t + \theta) = A e^{i (2\pi f t + \theta)} \]
so the magnitude of \(z(t)\) is the amplitude \(A\) and its angle is the phase \(2\pi f t + \theta\) at every time point.
In the frequency domain, the analytic signal contains only positive frequencies: the negative frequencies (the complex conjugates of the positive ones for a real signal) are removed and the positive ones doubled. This is how it is computed (e.g. by DSP.hilbert()):
fs = 100
t = (0:(200 - 1)) ./ fs
x = cos.(2π * 5 .* t)
X = fft(x)
n = length(x)
X[2:(n ÷ 2)] .*= 2 # positive frequencies
X[(n ÷ 2 + 2):end] .= 0 # negative frequencies
z_fft = ifft(X)
maximum(abs.(z_fft .- hilbert(x)))3.3306690738754696e-16The real part of the analytic signal is the original cosine, the imaginary part is a sine (shifted by −90°):
z = hilbert(x)
fig, ax = lines(t, real.(z); label = "real part: x(t)",
axis = (xlabel = "Time [s]", ylabel = "Amplitude"))
lines!(ax, t, imag.(z); color = :red, label = "imaginary part: H{x(t)}")
xlims!(ax, 0, 0.6)
axislegend(ax; position = :rb)
fig
and in the complex plane the analytic signal rotates on a circle of radius \(A\), once per cycle:
lines(t, real.(z), imag.(z);
axis = (type = Axis3, xlabel = "Time [s]", ylabel = "Real", zlabel = "Imaginary"))
From the analytic signal \(z(t)\):
In NeuroAnalyzer, htransform() (for signal vectors and arrays) and transform(obj; h = true) (for NEURO objects) return the analytic signal (c), amplitudes (a), powers (p) and phases (ph); frqinst(obj) returns the instantaneous frequency in Hz.
The instantaneous amplitude is the envelope of the signal. An amplitude-modulated signal: a 200 Hz carrier whose amplitude varies at 20 Hz:
fs = 1000
t = collect(0:(1 / fs):0.1)
s = @. (1 + cos(2π * 20 * t)) * sin(2π * 200 * t)
a = htransform(s).a
fig, ax = lines(t, s; color = :black, label = "signal",
axis = (xlabel = "Time [s]", ylabel = "Amplitude"))
lines!(ax, t, a; color = :red, label = "envelope")
lines!(ax, t, -a; color = :red)
axislegend(ax; position = :rt)
fig
The envelope is accurate when the amplitude changes more slowly than the oscillation itself (the spectra of the envelope and of the carrier do not overlap). The envelope is symmetric about zero, so a DC offset must be removed (e.g. remove_dc()) before the Hilbert transform; otherwise the offset is included in the envelope.
The instantaneous phase of a 10 Hz cosine increases linearly from −π to π in each cycle and then jumps back to −π (it is wrapped into \([-\pi, \pi]\)):
fs = 1000
t = collect(0:(1 / fs):1)
s = generate_cosine(10, t, 1)
ph = htransform(s).ph
plot_phase(ph, t; unit = :rad, type = :line, title = "Instantaneous Phase", xlabel = "Time [s]")
Phase unwrapping adds the appropriate multiples of \(2\pi\) at each jump, giving a continuous phase:
plot_phase(DSP.unwrap(ph), t; unit = :rad, type = :line, title = "Unwrapped Instantaneous Phase", xlabel = "Time [s]")
The unwrapped phase of a 10 Hz oscillation increases by \(2\pi \times 10\) radians per second. Note the deviation at the edges of the signal: the Hilbert transform (computed with the FFT) assumes a periodic signal and is distorted at the edges (edge effects); the first and last samples (or cycles) should be discarded.
The instantaneous frequency is the rate of change of the unwrapped phase:
\[ f(t) = \frac{1}{2\pi} \frac{\mathrm{d}\varphi(t)}{\mathrm{d}t} \approx \frac{f_s}{2\pi} \left(\varphi[n + 1] - \varphi[n]\right) \]
For a chirp whose frequency increases linearly from 2 to 20 Hz, the instantaneous frequency follows the sweep:
fs = 500
t = collect(0:(1 / fs):(4 - 1 / fs))
chirp = @. sin(2π * (2 * t + (20 - 2) / (2 * 4) * t^2))
ph = DSP.unwrap(htransform(chirp).ph)
fi = diff(ph) .* fs ./ (2π)
fig, ax = lines(t[1:(end - 1)], fi; label = "instantaneous frequency",
axis = (xlabel = "Time [s]", ylabel = "Frequency [Hz]"))
lines!(ax, t, 2 .+ (20 - 2) / 4 .* t; color = :red, linestyle = :dash, label = "true frequency")
ylims!(ax, 0, 25)
axislegend(ax; position = :lt)
fig
For NEURO objects use frqinst(obj; ch) (in Hz) and plot_fi():
obj = signal2obj(chirp, fs)
fi = frqinst(obj; ch = "s")[1, :, 1]
plot_fi(fi, t; title = "Instantaneous Frequency")
The analytic signal has a single amplitude, phase and frequency at each time point. They are meaningful only for narrowband (monocomponent) signals, dominated by one oscillation. For a mixture of oscillations, the instantaneous frequency jumps between the components and can even be negative or outside the frequency range of the signal:
fs = 500
t = collect(0:(1 / fs):(2 - 1 / fs))
mix = generate_sine(5, t, 1) .+ generate_sine(20, t, 0.8)
fi_mix = diff(DSP.unwrap(htransform(mix).ph)) .* fs ./ (2π)
NeuroAnalyzer.plot_signal(t[1:(end - 1)], fi_mix;
ylabel = "Frequency [Hz]",
title = "Instantaneous frequency of 5 Hz + 20 Hz")
Therefore, broadband signals such as EEG are band-pass filtered before the Hilbert transform (the filter-Hilbert method), or decomposed into narrowband components (e.g. intrinsic mode functions, see below). After band-pass filtering, the envelope, phase and frequency describe the oscillation in the chosen band:
obj = signal2obj(mix, fs)
obj_20 = NeuroAnalyzer.filter(obj; ch = "s", fprototype = :fir, ftype = :bp, cutoff = (15, 25), bw = 4)
fi_20 = frqinst(obj_20; ch = "s")[1, :, 1]
NeuroAnalyzer.plot_signal(t, fi_20;
ylabel = "Frequency [Hz]",
title = "Instantaneous frequency after 15–25 Hz band-pass filtering")
The frequency is 20 Hz except at the edges.
The alpha (8–13 Hz) amplitude envelope of an EEG channel:
eeg = import_edf("files/eeg/eeg.edf")
remove_dc!(eeg; ch = "eeg")
eeg_alpha = NeuroAnalyzer.filter(eeg; ch = "eeg", fprototype = :fir, ftype = :bp, cutoff = (8, 13), bw = 2)
h = transform(eeg_alpha; ch = "O1", h = true)seg = (30 * sr(eeg) + 1):(35 * sr(eeg)) # 30–35 s
t_seg = eeg.time_pts[seg]
fig, ax = lines(t_seg, eeg_alpha.data[get_channel(eeg; ch = "O1")[1], seg, 1]; color = :black, label = "O1, 8–13 Hz",
axis = (xlabel = "Time [s]", ylabel = "Amplitude [μV]"))
lines!(ax, t_seg, h.a[1, seg, 1]; color = :red, label = "envelope")
axislegend(ax; position = :rt)
fig
The envelope shows how the strength of the alpha rhythm changes over time (e.g. alpha bursts), which the power spectrum of the whole recording cannot show.
The Hilbert transform gives the amplitude of one narrowband signal over time. To obtain a time-frequency representation (power as a function of time and frequency), the signal is band-pass filtered around a series of frequencies and the instantaneous power is calculated for each frequency.
The Gaussian-Hilbert spectrogram filters the signal with a series of Gaussian-shaped band-pass filters (in the frequency domain, width gw Hz) and calculates the instantaneous power of each filtered signal with the Hilbert transform.
A chirp from 2 to 40 Hz:
fs = 1000
t = collect(0:(1 / fs):(2 - 1 / fs))
chirp = @. sin(2π * (2 * t + (40 - 2) / (2 * 2) * t^2))
NeuroAnalyzer.plot_signal(t, chirp)
spec_data = ghtspectrogram(chirp; fs = fs, db = false)f_idx = spec_data.f .<= 60
plot_heatmap(spec_data.p[f_idx, :];
x = spec_data.t,
y = spec_data.f[f_idx],
xlabel = "Time [s]",
ylabel = "Frequency [Hz]",
cb_title = "Power [AU²]")
For NEURO objects use spectrogram(obj; ch, method = :gh).
The Hilbert-Huang transform (HHT) avoids predefined filters: the signal is decomposed into intrinsic mode functions (IMFs) by empirical mode decomposition (EMD), and the Hilbert transform is applied to each IMF.
The HHT is adaptive (data-driven) and suited for non-stationary and non-linear signals. See also Empirical Mode Decomposition.
Two sines (5 and 20 Hz):
fs = 200
t = collect(0:(1 / fs):5)
s = generate_sine(5, t, 0.7) .+ generate_sine(20, t, 1.0)
obj = signal2obj(s, fs)
NeuroAnalyzer.plot_signal(t, s)
IMFs (the last row is the residue):
imf = emd(obj; ch = "s", ep = 1)
plot_imf(imf; t = t)
Hilbert-Huang spectrogram:
spec_data = hhtspectrogram(s, t; fs = fs, db = false)
f_idx = spec_data.f .<= 40
plot_heatmap(spec_data.p[f_idx, :];
x = spec_data.t,
y = spec_data.f[f_idx],
xlabel = "Time [s]",
ylabel = "Frequency [Hz]",
cb_title = "Power [AU²]")
The two components appear at 5 and 20 Hz; the edge artifacts of EMD and of the Hilbert transform are visible at the beginning and the end of the signal.
The Hilbert marginal spectrum sums the Hilbert spectrum over time, giving a frequency spectrum comparable to a power spectrum:
hms = hmspectrum(obj; ch = "s")
f_idx = hms.f .<= 40
plot_psd(hms.f[f_idx], hms.p[f_idx, 1];
xlabel = "Frequency [Hz]",
ylabel = "Hilbert marginal spectrum")
For comparison, a short-time Fourier transform spectrogram (1-second windows) of the same signal:
spec_data = NeuroAnalyzer.spectrogram(obj; ch = "s", method = :stft, db = false)
f_idx = spec_data.f .<= 40
plot_heatmap(spec_data.p[f_idx, :, 1, 1];
x = spec_data.t,
y = spec_data.f[f_idx],
xlabel = "Time [s]",
ylabel = "Frequency [Hz]",
cb_title = "Power [AU²/Hz]")