NeuroAnalyzer tutorials: Analyze: Frequency domain: Power spectrum

Initialize NeuroAnalyzer

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

Spectral energy: \[ E_s = \vert X(t) \vert^2 \] \(X(t)\): Fourier transform of \(x(t)\)

Spectrum: amplitude of signal as a function of frequency.

Spectral density: describes the distribution of power into frequency components composing that signal.

Spectral density = power spectrum = power spectral density (PSD).

Units (for EEG):

amplitude \(\mathrm{μV}\)
power: \(\mathrm{μV^2}\)
PSD: \(\mathrm{μV^2/Hz}\)

Spectrogram: spectrum changing over time.

Power spectrum shows power (strength) of the signal at each of the component frequencies of the signal. Therefore, it converts the signal from time domain to frequency domain.

Methods of estimating power:

Welch periodogram
Bartlett’s method
multitaper
wavelets

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

Signal analysis uncertainty principle:

\[ \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.

Welch periodogram

Welch’s method (periodogram): split signal into smaller windows, compute PSD for each window and finally average all individual power spectra.

Periodogram has high variance = noisy spectrogram

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

Plot Welch PSD:

plot_psd(eeg,
         ch="Fp1")

Fast Fourier transform

Plot FFT PSD:

plot_psd(eeg,
         ch="Fp1",
         method=:fft)

Short time Fourier transform

For PSD, short time Fourier transform is like the Welch’s method.

Plot short time Fourier transform PSD:

plot_psd(eeg,
         ch="Fp1",
         method=:stft)

Multi-taper

Fourier transform is performed several times, each times on data tapered using a different taper.

Plot multi-taper PSD:

plot_psd(eeg,
         ch="Fp1",
         method=:mt)

(!) To calculate recommended number of tapers, use ntapers():

plot_psd(eeg,
         ch="Fp1",
         method=:mt,
         nt=ntapers(eeg, df=1))

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

Morlet wavelet

Plot Morlet wavelet PSD:

plot_psd(eeg,
         ch="Fp1",
         method=:mw)

Scaling axes

Plot PSD, scale x-axis:

plot_psd(eeg,
         ch="Fp1",
         ax=:loglin)

Multi-channel spectrum

Plot multi-channel PSD:

plot_psd(eeg,
         ch=["Fp1", "Fp2", "F3", "F4"])

Butterfly plot of multi-channel PSD:

plot_psd(eeg,
         ch=["Fp1", "Fp2", "F3", "F4"],
         type=:butterfly)

Plot multi-channel PSD, mean ± 95CI:

plot_psd(eeg,
         ch=["Fp1", "Fp2", "F3", "F4"],
         type=:mean)

Plot PSD 3d waterfall:

plot_psd(eeg,
         ch=["Fp1", "Fp2", "F3", "F4"],
         method=:mw,
         type=:w3d)

Plot PSD 3d surface:

plot_psd(eeg,
         ch=["Fp1", "Fp2", "F3", "F4"],
         type=:s3d)

Topographical map of PSD

Plot topomap of PSDs:

plot_psd(eeg,
         ch="eeg",
         type=:topo)

Phase spectral density

Plot phase spectral density:

plot_phsd(eeg,
          ch="eeg")

(!) For plot_phsd(), the same plot types are available as for plot_psd().