NeuroAnalyzer tutorials: Analyze: Time-frequency domain: Spectrogram

Initialize NeuroAnalyzer 
using ContinuousWavelets
using NeuroAnalyzer
using Plots
eeg = load("files/eeg.hdf");

Spectrogram: spectrum changing over time.

Spectrogram: a visual representation of data in the time-frequency domain.

Brain signal is non-stationary (it changes its properties over time). Therefore spectrogram (changes of power over time: frequencies vs time, color coded for power) is a better representation than power spectrum.

Spectrograms are typically plotted with time on the X axis and frequency on the Y axis. Signal measurements (e.g. amplitude or power) at each time-frequency point are represented at the corresponding point on the plot by mapping their value to a color.

Spectrograms are usually calculated in time segments. Time segments should be tapered to reduce edge artifacts.

Time resolution of spectrogram equals segments width.

Frequency resolution of spectrogram depends on the frequency range between elements; equals half of the segments width.

\[ tr = n / r \]

\[ fr = r / n \]

\(r\): sampling rate

\(n\): segment length

Signal duration required for frequency resolution: to distinguish two signals by 0.1 Hz difference: 1 / 0.1 Hz = 10 seconds

Temporal precision = number of chunks / signal length.

Gabor limit: trade-off between time and frequency resolution, you can’t have both (longer window causes lower frequency resolution and higher time resolution).

\[ \Delta t \Delta f = \frac{1}{4 \pi} \]

  • narrow window: poor frequency resolution, good time resolution
  • wide window: good frequency resolution, poor time resolution

High frequencies required good time frequency (they happen more often), while low frequencies require frequency resolution.

Short segments gives better temporal precision, but poor frequency resolution.

Conversely, long segments gives better frequency precision (since FFT precision depends on the number of time points), but poor temporal precision.

Segment overlap: 50-90% is recommended.

Segment duration required for frequency resolution: depends on the lowest frequency to be captured, e.g. for 3 Hz:

\[ 1 / 3 = 0.333 \]

in seconds

0.1 Hz = 10 cycle per second = 1 cycle length = 10000 ms
1 Hz = 1 cycle per second = 1 cycle length = 1000 ms
3 Hz = 3 cycles per second = 1 cycle length = 333 ms

Optimally, segments length should be 2-3-times higher then recommended minimum.

Relationship between temporal and frequency precisions: \(\downarrow\) temporal precision \(\leftrightarrow\) \(\uparrow\) frequency resolution.

Solution: use different widths for different frequencies:

  • low frequencies - longer window
  • high frequencies - shorter window

ECG signal: within 1-2 Hz range.

Short-time Fourier transform

Short-time Fourier transform (STFT): compute the Fourier transform in short, successive time windows, rather than once over the entire time series. Each window is assumed to be stationary.

Invented by Dennis Gabor.

Segments are assumed to be stationary.

Hopping length usually is half of tapering window length.

Algorithm: split signal into Short windows, perform FFT of each time segment, compute PSD of this time window and rotate PSD by 90 degree. Each PSD becomes one vertical color bar of the spectrogram (along the X axis, time). Successive spectra displayed along time show the evolution of frequencies over time.

Plot STFT spectrogram:

plot_spectrogram(eeg,
                 ch="Fp1",
                 db=true)

(!) To smooth the histogram, use smooth=true option and n parameter (kernel size of the Gaussian blur; larger kernel means more smoothing).

plot_spectrogram(eeg,
                 ch="Fp1",
                 db=true,
                 smooth=true,
                 n=5)

Multi-taper spectrogram

Single taper spectrum has high variance = noisy spectrogram.

  • no tapering: highest noise (high variance)
  • single taper spectrum: less of noise, but still low resolution

Solution: use multi-taper spectral analysis

  • less noise (low variance)
  • better resolution

For many trials: average them to reduce noise.

For single trial: use multiple tapers (DPSS or Slepian sequence).

DPSS: tapers come from Slepian sequence (discrete prolate spheroidal sequence):

  • tapers that keep most of the power in the main lobe and minimize spectral leakage
  • each taper is orthogonal to each other taper = produces uncorrelated projections of data, these projections can be averaged to reduce variability
  • reduce broadband bias (tapers remove power from the side lobes)
  • very clear main lobe
  • less variance
  • some of DPSS tapers increase near the window edges (Hanning is zero at both ends): better representation of the signal edges

DPSS gives smooth spectrogram:

  • good if a lot of noise
  • poor temporal precision, especially for high frequencies (for signal > 30 Hz better signal-noise ratio)
  • difficult to interpret for frequencies < 30 Hz

Algorithms:

  • generate DPSS tapers
  • multiplying data by DPSS tapers produces tapered data (one vector per each taper)
  • estimate single-taper spectrum from each taper
  • compute the mean single-taper spectrum (average across tapers)

Individual tapers are separate representations of power at various time points.

Plot multi-taper spectrogram:

plot_spectrogram(eeg,
                 ch="Fp1",
                 db=true,
                 method=:mt)

How to determine number of tapers: \[ TW = X \]

\(T\): duration of the recording

\(2W\): desired frequency resolution (resolution bandwidth)

\[ TW = (N \times \Delta F) / 2 \]

\(N\): size of the data window (s); set by determining the length of time over which the signal is assumed to appear stationary; how many seconds of data will be assessed

\(\Delta F\): frequency resolution (bandwidth); set by assuming oscillatory structure of the data; bandwidth of the main lobe = minimum distance between frequency peaks to resolve; minimum difference between two frequencies that can be separated using FFT, smallest distance between frequency peaks that we want to observe (e.g. 1 Hz)

\(TW\): time-half-bandwidth product

\(L\): number of tapers to use

\[ L = (2 \times TW) - 1 \]

E.g. for 10 seconds signal and 2 Hz spectrum resolution:

\(N = 10\)

\(\Delta F = 2\)

\(TW = 10\)

\[ L = (2 \times 10) - 1 = 19 \]

↑ \(L\) → better variance reduction (less noise)

↑ \(\Delta F\) → wider main lobe bandwidth = more blurry spectrogram

More tapers does not improve the spectrum estimator and may lead to spurious effects.

To calculate recommended number of tapers:

ntapers(eeg, df=1)
1204

where df is frequency resolution (bandwidth) or the smallest distance between frequency peaks that we want to observe (e.g. 1 Hz).

Morlet wavelet

For more details on wavelet transformation, see the introduction to wavelet transformation.

Morlet wavelet: sine wave tapered with Gaussian.

fs = 100
t = 0:1/fs:2
f = 10
s = generate_sine(f, t, 1, 90)
plot_signal(t, s)

Generate 10 Hz Gaussian:

gt = t[end] / 2
g = generate_gaussian(fs, f, gt)
plot_signal(t, g)

(!) Gaussian length is -t:1/fs:t.

Generate Morlet wavelet:

mw = s .* g
plot_signal(t, mw)

In the frequency domain (FFT): looks like Gaussian, with peak at the wavelet frequency, width of the Gaussian used to taper the wavelet.

_, _, p, _ = ftransform(mw)
f, _ = freqs(mw, fs)
f = round.(f, digits=1)
f_idx = vsearch(20, f)
plot_psd(f[1:f_idx], p[1:f_idx])
Plots.vline!([10], ls=:dash, lc=:red)

Frequency response of a Morlet wavelet is always Gaussian.

Plot Morlet wavelet spectrogram:

plot_spectrogram(eeg,
                 ch="Fp1",
                 db=true,
                 method=:mw)

Gaussian-Hilbert

Detailed introduction to Gaussian-Hilbert spectrogram can be in the introduction to Hilbert transform.

Plot Gaussian-Hilbert spectrogram:

plot_spectrogram(eeg,
                 ch="Fp1",
                 method=:gh)

Hilbert-Huang

Detailed introduction to Hilbert-Huang spectrogram can be in the introduction to Hilbert transform.

Plot Hilbert-Huang spectrogram:

plot_spectrogram(eeg,
                 ch="Fp1",
                 method=:hht)

Continuous wavelet transformation

For more details on wavelet transformation, see the introduction to wavelet transformation.

Plot spectrogram:

plot_spectrogram(eeg,
                 ch="Fp1",
                 method=:cwt)

(!) Default wavelet is wavelet(Morlet(2π), β=2).

Multi-channel spectrogram

In multi-channel spectrograms, X axis is the frequency range and Y axis is the list of channels. For each channel its spectrogram is calculated and averaged to a single row in the powers matrix.

Plot multi-channel spectrogram:

plot_spectrogram(eeg,
                 ch="eeg",
                 db=true,
                 frq_lim=(0, 50))

Thresholding

To mark a thresholded region:

plot_spectrogram(eeg,
                 ch="Fp1",
                 db=true,
                 frq_lim=(0, 50),
                 smooth=true,
                 n=5,
                 threshold=0,
                 threshold_type=:geq)

(!) threshold parameters sets the region (0 dB in that case), the following types of thresholding are available:

  • :eq: draw region if values are equal to threshold
  • :neq: draw region if values are not equal to threshold
  • :geq: draw region if values are ≥ to threshold
  • :leq: draw region if values are ≤ to threshold
  • :g: draw region if values are > to threshold
  • :l: draw region if values are < to threshold