using DSP
using Statistics
using CairoMakie
using NeuroAnalyzer
eeg = load("files/eeg.hdf")
NeuroAnalyzer.filter!(eeg,
fprototype = :fir,
ftype = :bp,
cutoff = (5, 25),
order = 91,
ch = "all")
e2 = epoch(eeg, ep_len = 2)
keep_epoch!(e2, ep = 11:20)
detrend!(e2,
type = :mean,
ch = "all")Objective
The Hilbert transform is used to derive the analytic signal from a real-valued signal. This allows you to calculate instantaneous properties of the signal, such as amplitude, phase, and frequency, at any point in time.
Analytic Signal
An analytic signal is a complex-valued signal that contains no negative-frequency components. It is constructed from a real signal using the Hilbert transform.
Key Property: The analytic signal’s Z-transform is zero on the lower half of the unit circle, meaning it contains only positive frequencies.
For a real signal \(f(t)\), the analytic signal \(z(t)\) is created using the Hilbert transform \(H\{f(t)\}\):
\[ z(t) = f(t) + iH\{f(t)\} \]
where:
\(z(t)\) is the analytic signal (complex-valued),
\(f(t)\) is the original real signal,
\(iH\{f(t)\}\) is the imaginary part of the analytic signal, derived from the Hilbert transform.
# example real signal (e.g., EEG data)
x = randn(1000)
# compute the analytic signal
z = hilbert(x)
# extract instantaneous amplitude and phase
instantaneous_amplitude = abs.(z)
instantaneous_phase = DSP.angle.(z)fig = Figure()
ax = Axis(fig[1, 1],
title = "Instantaneous Amplitude")
lines!(ax,
instantaneous_amplitude,
label = "Instantaneous Amplitude")
fig
fig = Figure()
ax = Axis(fig[1, 1],
title="Instantaneous Phase")
lines!(ax,
instantaneous_phase,
label="Instantaneous Amplitude")
fig
Hilbert transform works best for narrowband signals (signals with all energy centered about a single frequency).
Power and phase at each time point are from the frequency with the most power at each time point. Therefore, for changes of phase angles to be interpretable, the signal must be narrowband.
The procedure requires that the signal be mono-component: signal has to be band-pass filtered for better estimation of instantaneous frequency.
NeuroAnalyzer transform() and htransform() provide Hilbert transformation and return instantaneous amplitudes, powers and phases.
The amplitude of the analytic signal is the envelope of the original signal.
fs = 1000
t = collect(0:1/1000:0.1)
s = @. (1 + cos(2pi * 20 * t)) * sin(2pi * 200 * t)
# Hilbert transform
a = htransform(s).a
# plot the instantaneous amplitudes over time
fig = Figure()
ax = Axis(fig[1, 1],
title = "Instantaneous Amplitude")
lines!(ax,
t,
s,
color = :black)
lines!(ax,
t,
a,
color = :red)
lines!(ax,
t,
-a,
color = :red)
fig
For other signals, envelopes will be inaccurate. Also, since the lower envelope is symmetric to the upper one with respect to the X axis, signal should not be shifted along the Y axis.
fs = 1000
t = collect(0:(1 / fs):1)
s = generate_cosine(10, t, 1)
w = generate_window(:bh, length(s))
sw = s .* w
# plot the signal
NeuroAnalyzer.plot(t, sw)
# Hilbert transform
ph = htransform(s).ph
# plot the instantaneous phases over time
NeuroAnalyzer.plot_phase(ph,
t,
unit = :rad,
type = :line,
title = "Instantaneous Phase",
xlabel = "Time [s]")
Phase unwrapping ensures that all appropriate multiples of \(2 \pi\) have been included in phase.
NeuroAnalyzer.plot_phase(DSP.unwrap(ph),
t,
unit = :rad,
type = :line,
title = "Unwrapped Instantaneous Phase",
xlabel = "Time [s]")
Definition
Instantaneous frequency represents the time-varying frequency of a signal at each point in time. It is derived from the temporal derivative of the phase angle of the analytic signal:
\[ \omega(t) = \frac{\mathrm{d}\varphi(t)}{\mathrm{d}t} \]
where:
\(\omega(t)\) is the instantaneous frequency (angular velocity, in radians per second),
\(\varphi(t)\) is the instantaneous phase angle (in radians).
Phase Unwrapping
Before computing the derivative, the phase angles must be unwrapped to avoid discontinuities:
Interpretation and Caveats
fi = frqinst(sw) .* fs
NeuroAnalyzer.plot_fi(fi,
t,
title = "Instantaneous Frequency")
Tip: frqinst() returns frequencies in in cycles per sample. To obtain frequencies in Hz, multiply by the sampling rate.
The Hilbert transform itself does not directly provide a frequency spectrum like the Fourier transform. Instead, it provides time-varying amplitudes and phases, which can be used to reconstruct or analyze the spectrum.
Estimating the Spectrum from Instantaneous Amplitudes
The instantaneous amplitude \(A(t)\) represents the envelope of the signal at each time point.
To estimate the spectrum from \(A(t)\), you can:
This gives you an estimate of how power is distributed across frequencies.
Plotting Hilbert instantaneous powers:
t = e2.epoch_time
X = NeuroAnalyzer.transform(e2, ch = "F3", h = true, db = true)
NeuroAnalyzer.plot(t,
X.p[1, :, 1], # there is only one channel in the data
xlabel = "t",
ylabel = "P(t)")
Computing the power spectral density (PSD) using Fourier transform:
X = NeuroAnalyzer.transform(e2, ch = "F3", h = true, db = true)
psd_data = psd(X.p[1, :, 1],
fs=sr(e2))
NeuroAnalyzer.plot_psd(psd_data.f, # frequencies
psd_data.p, # powers
xlabel = "Frequency (Hz)",
ylabel = "Power")
Limitations and Considerations
For a more accurate time-frequency representation, consider using the Hilbert-Huang Transform.
For a more robust spectral estimate, the Hilbert-Huang Transform (HHT) combines the Empirical Mode Decomposition (EMD) with the Hilbert transform:
The HHT is particularly useful for non-stationary and non-linear signals, such as EEG or biological time series.
Plotting signal (channel Fp3, 1st epoch):
ep = 1
ch = "F3"
NeuroAnalyzer.plot(e2,
ch = ch,
ep = ep,
gui = false)
Computing and plotting Hilbert-Huang spectrogram:
spec_data = NeuroAnalyzer.spectrogram(e2,
ch = ch,
method = :hht,
db = true)
NeuroAnalyzer.plot_heatmap(spec_data.p[:, :, 1, ep],
x = spec_data.t,
y = spec_data.f,
xlabel = "Time [s]",
ylabel = "Frequency [Hz]",
cb_title = "Power [dB AU^2/Hz]")
Tip: As other techniques based on the Hilbert transform, Hilbert-Huang spectrogram works best for a narrow-band signal.
Hilbert marginal spectrum is obtained by integrating the Hilbert spectrum over time, providing a way to analyze non-linear and non-stationary signals by reducing the Hilbert spectrum (a time-frequency map) to a single frequency-amplitude plot.
hms, f = hmspectrum(e2, ch = ch)
NeuroAnalyzer.plot_psd(f,
hms[:, ep],
xlabel = "Frequency [Hz]",
ylabel = "Hilbert marginal spectrum")
Example signal, two sines (5 and 20 Hz) mixed:
fs = 200
t = 0:(1 / fs):5
s1 = generate_sine(5, t, 0.7)
s2 = generate_sine(20, t, 1.0)
s = s1 .+ s2
NeuroAnalyzer.plot(t, s)
Computing and plotting IMFs:
imf = emd(s, t)
NeuroAnalyzer.plot_imf(imf,
t = t)Plotting Hilbert-Huang spectrogram:
spec_data = hhtspectrogram(s, t, fs = fs, db = false)
# plot up to 40 hz
f_idx = vsearch(40, spec_data.f)
NeuroAnalyzer.plot_heatmap(spec_data.p[1:f_idx, :];
x = spec_data.t,
y = spec_data.f[1:f_idx],
xlabel = "Time [s]",
ylabel = "Frequency [Hz]",
cb_title = "Power [AU^2/Hz]")
Please note the edge artifacts in the reconstructed instantaneous frequency time series.
Computing and plotting a short-time Fourier transform spectrogram:
spec_data = NeuroAnalyzer.spectrogram(s,
fs = fs,
db = false,
method = :stft)
# plot up to 40 hz
f_idx = vsearch(40, spec_data.f)
NeuroAnalyzer.plot_heatmap(spec_data.p[1:f_idx, :],
x = spec_data.t,
y = spec_data.f[1:f_idx],
xlabel = "Time [s]",
ylabel = "Frequency [Hz]",
cb_title = "Power [AU^2/Hz]")
The Gaussian-Hilbert Spectrogram is a powerful tool for time-frequency analysis, combining the Gaussian window with the Hilbert transform to provide a high-resolution, localized representation of a signal’s frequency content over time.
Gaussian-Hilbert Spectrogram uses a series of Gaussian tapers to extract individual frequencies and calculate power per each frequency using the Hilbert transform.
It is particularly useful for non-stationary signals (such as EEG), where frequency content changes over time.
Chirp from 10 to 20 Hz:
fs = 1000
t = collect(0:(1 / fs):2)
n = length(t)
# frequencies from 2 to 40 Hz
f = [2, 40]
ff = linspace(f[1], f[2] * Statistics.mean(f) / f[2], n)
# generate chirp
s = @. sin(2 * pi * ff * t)
NeuroAnalyzer.plot(t, s)
Computing Gaussian-Hilbert spectrogram:
spec_data = ghtspectrogram(s, fs = fs, db = false)Plotting Gaussian-Hilbert spectrogram for the 0 to 100 Hz frequency range:
f1 = vsearch(0, spec_data.f)
f2 = vsearch(100, spec_data.f)
plot_heatmap(spec_data.p[f1:f2, :],
x = spec_data.t,
y = spec_data.f[f1:f2],
xlabel = "Time [s]",
ylabel = "Frequency [Hz]",
cb_title = "Power [AU^2]")