Sampling

Initialize NeuroAnalyzer 
using NeuroAnalyzer
using CairoMakie
CairoMakie.activate!(type = "png")
using FFTW
# Makie figure and axis settings matching NeuroAnalyzer plots
fig_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 figures
ax_style_small = merge(ax_style, (titlesize = 14,
                                  xticks = LinearTicks(5),
                                  xminorticks = IntervalsBetween(5)))
# amplitude spectrum of a real signal
function amp_spectrum(x::AbstractVector, fs::Real)
    n = length(x)
    return rfftfreq(n, fs), 2 .* abs.(rfft(x)) ./ n
end

Sampling converts a continuous signal (e.g. the voltage on the scalp) into a discrete signal: the signal is measured at regular time intervals by an analog-to-digital converter (ADC). Digitization has two aspects:

  1. Sampling (discretization in time): the signal is measured every \(\Delta t\) seconds,
  2. Quantization (discretization in amplitude): each measurement is rounded to one of a finite set of values.

Sampling rate and sampling interval

The sampling rate \(f_s\) (Hz) is the number of samples per second; the sampling interval is the time between consecutive samples:

\[ \Delta t = \frac{1}{f_s} \]

e.g. 256 Hz → 3.9 ms, 1000 Hz → 1 ms. The sampling interval is also the precision with which the timing of events (e.g. ERP latencies) can be measured.

Nyquist Frequency and Aliasing

Nyquist–Shannon sampling theorem: a signal containing no frequencies above \(f_{max}\) is completely determined by its samples if

\[ f_s > 2 f_{max} \]

The highest frequency that can be represented at a given sampling rate is the Nyquist frequency:

\[ f_N = \frac{f_s}{2} \]

e.g. 128 Hz for \(f_s = 256\) Hz.

Aliasing: a frequency \(f\) above the Nyquist frequency is not lost, but appears in the sampled signal as a lower frequency (an alias): the frequency is “folded” back into the range 0 to \(f_N\):

\[ f_{alias} = \vert f - k f_s \vert \]

where \(k\) is the integer for which \(f_{alias} \le f_N\). For example, a 100 Hz sine sampled at 80 Hz appears as a 20 Hz sine, and a 90 Hz sine sampled at 64 Hz as a 26 Hz sine.

A 9 Hz sine sampled at 50 Hz (above the Nyquist rate) and at 10 Hz (below it):

f = 9
# "continuous" signal: very fine time grid
t_c = 0:0.0005:1
fig = Figure(; size = (900, 600))
for (idx, fs) in enumerate((50, 10))
    t_s = 0:(1 / fs):1
    f_alias = round(Int, abs(f - round(f / fs) * fs))
    ttl = "f_s = $fs Hz" *
          (fs > 2 * f ? "" : ": the samples show a $f_alias Hz sine")
    ax = Axis(fig[idx, 1];
              ax_style_small...,
              title = ttl,
              xlabel = "Time [s]",
              ylabel = "Amplitude")
    lines!(ax, t_c, sin.(2π * f .* t_c); color = :gray,
           label = "$f Hz signal")
    fs < 2 * f && lines!(ax, t_c, sin.(2π * (f - fs) .* t_c);
                         color = :red, linestyle = :dash,
                         label = "alias ($f_alias Hz)")
    scatter!(ax, t_s, sin.(2π * f .* t_s); color = :black,
             markersize = 8, label = "samples")
    axislegend(ax; position = :rt)
end
fig

At 10 Hz, the samples of the 9 Hz sine are indistinguishable from the samples of a 1 Hz sine. Once a signal has been sampled, an alias cannot be distinguished from a real low-frequency component: aliasing is irreversible.

Anti-aliasing filter: to prevent aliasing, frequencies above the Nyquist frequency must be removed before sampling, with an analog low-pass filter. EEG amplifiers have such a filter built in; its cutoff is usually set well below the Nyquist frequency (e.g. 1/3–1/4 of the sampling rate), because analog filters have a gradual roll-off.

Choosing the sampling rate

  • The sampling rate must exceed twice the highest frequency of interest; in practice it is set to at least 2.5–4× this frequency, to leave room for the roll-off of the anti-aliasing filter.
  • For EEG analyses up to ≈ 40–50 Hz, 200–256 Hz is sufficient; high-frequency activity (high gamma, ripples) requires 1000 Hz or more, and spike recordings tens of kHz.
  • Higher sampling rates also give a finer timing of events (e.g. 1 ms at 1000 Hz), at the cost of larger files and longer computations.

Quantization

The ADC represents each sample with a finite number of bits: an \(N\)-bit converter distinguishes \(2^N\) levels. The step between levels (the least significant bit, LSB) is:

\[ \text{LSB} = \frac{\text{input range}}{2^N} \]

Rounding to the nearest level introduces a quantization error of at most ±½ LSB. For a full-scale sine, the resulting signal-to-quantization-noise ratio is ≈ \(6.02 N + 1.76\) dB.

A sine quantized with 3 bits (8 levels):

t = 0:0.001:1
x = sin.(2π * 2 .* t)
n_bits = 3
lsb = 2 / 2^n_bits
# 2^n_bits levels, at the centers of the quantization steps
x_q = clamp.((floor.(x ./ lsb) .+ 0.5) .* lsb, -1 + lsb / 2, 1 - lsb / 2)
fig, ax = lines(t, x; color = :gray, label = "signal",
                figure = (size = fig_size,),
                axis = (ax_style...,
                        xlabel = "Time [s]",
                        ylabel = "Amplitude"))
stairs!(ax, t, x_q; color = :red, step = :center,
        label = "quantized ($n_bits bits)")
lines!(ax, t, x .- x_q; color = :blue, label = "quantization error")
axislegend(ax; position = :rt)
fig

Modern EEG amplifiers use 16–24-bit converters: with a 24-bit ADC and an input range of ±400 mV, 1 LSB is ≈ 0.05 μV, far below the amplitude of the EEG and the noise of the electrodes.

Frequency Resolution

The spectrum of a signal of \(N\) samples (duration \(T = N / f_s\)) contains frequencies from 0 to the Nyquist frequency in steps of:

\[ \Delta f = \frac{f_s}{N} = \frac{1}{T} \]

so it has \(N/2 + 1\) frequency bins (for even \(N\)). The frequency resolution depends only on the duration of the signal, not on the sampling rate: 2 s of data give 0.5 Hz, 10 s give 0.1 Hz. The lowest non-zero frequency, \(1/T\), is also called the Rayleigh frequency: a signal must contain at least one full cycle of a frequency for that frequency to be resolved.

freqs() returns the frequencies of the spectrum and the Nyquist frequency, for a number of samples and a sampling rate or for a NEURO object:

fs = 256
n = 2 * fs
hz, nqf = freqs(n, fs)
println("number of frequencies: ", length(hz))
println("frequency step:        ", hz[2] - hz[1], " Hz")
println("Nyquist frequency:     ", nqf, " Hz")
number of frequencies: 257
frequency step:        0.5 Hz
Nyquist frequency:     128.0 Hz

The index of a frequency is found with vsearch() (index of the nearest value):

idx = vsearch(10.3, hz)
println("frequency nearest to 10.3 Hz: ", hz[idx], " Hz (index ", idx, ")")
frequency nearest to 10.3 Hz: 10.5 Hz (index 22)

Resampling

Resampling changes the sampling rate of a recorded signal, e.g. to:

  • reduce the size of the data and the computation time (downsampling),
  • give recordings from different devices a common sampling rate,
  • match the sampling rate expected by a method or a model.

Downsampling

Keeping every \(k\)-th sample reduces the sampling rate \(k\) times, and with it the Nyquist frequency. Frequencies between the new and the old Nyquist frequency are then aliased, exactly as in sampling. Downsampling must therefore be preceded by a digital anti-aliasing low-pass filter with the cutoff at (or below) the new Nyquist frequency.

A signal with a 10 Hz and a 90 Hz component, sampled at 256 Hz and downsampled to 64 Hz (new Nyquist frequency: 32 Hz):

fs = 256
t = collect(0:(1 / fs):(10 - 1 / fs))
x = sin.(2π * 10 .* t) .+ 0.5 .* sin.(2π * 90 .* t)
# naive downsampling: every 4th sample, no filter
x_naive = x[1:4:end]
# resample(): anti-aliasing filter + downsampling
obj = create_object(; data_type = "eeg")
create_data!(obj; data = reshape(x, 1, :, 1), fs = fs, type = "eeg")
x_rs = resample(obj; new_sr = 64).data[1, :, 1]
fig = Figure(; size = (900, 600))
spectra = ((x, fs, "original, 256 Hz"),
           (x_naive, 64, "every 4th sample, 64 Hz"),
           (x_rs, 64, "resample(), 64 Hz"))
for (idx, (s, f_s, ttl)) in enumerate(spectra)
    f, a = amp_spectrum(s, f_s)
    ax = Axis(fig[idx, 1];
              ax_style_small...,
              title = ttl,
              xlabel = "Frequency [Hz]",
              ylabel = "Amplitude")
    lines!(ax, f, a)
    xlims!(ax, 0, 100)
end
fig

Without the filter, the 90 Hz component appears at \(\vert 90 - 64 \vert = 26\) Hz, a frequency that is not present in the signal; resample() removes it.

Upsampling

Upsampling inserts new samples between the existing ones (in practice: zeros are inserted and the result is low-pass filtered, which interpolates the signal). Upsampling does not add information: the spectrum above the original Nyquist frequency stays empty.

Resampling in NeuroAnalyzer

resample() changes the sampling rate of all channels to new_sr, calling upsample() or downsample() as needed. Resampling is done with a polyphase filter for the exact ratio of the sampling rates (e.g. 256 → 250 Hz is the ratio 125/128), with a built-in anti-aliasing low-pass filter. The time axes are recalculated; markers are stored in seconds and do not change.

eeg = import_edf("files/eeg/eeg.edf")
eeg_128 = resample(eeg; new_sr = 128)
eeg_512 = resample(eeg; new_sr = 512)
eeg_250 = resample(eeg; new_sr = 250)
for obj in (eeg, eeg_128, eeg_512, eeg_250)
    println(lpad(sr(obj), 4), " Hz: ",
            lpad(signal_len(obj), 6), " samples, ",
            signal_duration(obj), " s, Nyquist frequency ",
            freqs(obj)[2], " Hz")
end
 256 Hz: 308480 samples, 1205.0 s, Nyquist frequency 128.0 Hz
 128 Hz: 154240 samples, 1205.0 s, Nyquist frequency 64.0 Hz
 512 Hz: 616960 samples, 1205.0 s, Nyquist frequency 256.0 Hz
 250 Hz: 301250 samples, 1205.0 s, Nyquist frequency 125.0 Hz

Up- and downsampling back to the original rate restores the signal (except for the frequencies above the lower Nyquist frequency):

x = eeg.data[1, :, 1]
x_back = resample(eeg_512; new_sr = 256).data[1, :, 1]
# middle of the signal (edges excluded)
m = 1000:(length(x) - 1000)
println("relative error 256 → 512 → 256 Hz: ",
        round(sqrt(sum(abs2, x_back[m] .- x[m]) / sum(abs2, x[m]));
              sigdigits = 2))
relative error 256 → 512 → 256 Hz: 0.00029

resample!(), upsample!() and downsample!() modify the object in place.

Practical notes

  • Resample the continuous signal, before epoching: each epoch is resampled separately, and the filter edges distort the beginning and the end of short epochs.
  • When downsampling, check that the frequencies of interest stay well below the new Nyquist frequency (the anti-aliasing filter attenuates frequencies near it).
  • Downsampling early in a pipeline (e.g. from 1000 Hz to 250 Hz for analyses below 50 Hz) speeds up all subsequent steps, e.g. ICA.
  • When combining recordings with different sampling rates, downsample to the lowest rate (upsampling the others adds no information) or resample all to a common rate.
  • Event timing is limited by the sampling interval of the resampled signal: marker times are kept, but samples exist only every \(1 / f_s\).