Tapering

Initialize NeuroAnalyzer

using NeuroAnalyzer
using CairoMakie
CairoMakie.activate!(type = "png")
using FFTW
using Statistics
NeuroAnalyzer.progress_bar = false

# 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 windowed signal (zero-padded to n_fft points),
# scaled so that a sine wave of amplitude A gives a peak of A
function amp_spectrum(s, w, fs; n_fft = 8 * length(s))
    x = vcat(s .* w, zeros(n_fft - length(s)))
    a = abs.(fft(x))[1:(n_fft ÷ 2 + 1)] .* 2 ./ sum(w)
    f = collect(0:(n_fft ÷ 2)) .* fs ./ n_fft
    return f, a
end

Tapering multiplies a signal by a window (a taper), a function that is largest in the middle and falls towards zero at both ends. It is used before spectral analysis (power spectra, spectrograms) to reduce spectral leakage.

Spectral Leakage

The Fourier transform of a finite segment treats it as one period of a periodic signal: the end of the segment is joined to its start. If an oscillation does not complete a whole number of cycles in the segment, the join is a discontinuity, and the power of the oscillation spreads (“leaks”) to other frequencies.

Cutting a segment is the same as multiplying the infinite signal by a rectangular window (1 inside the segment, 0 outside). Multiplication in time is convolution in frequency, so the spectrum of every oscillation is smeared by the spectrum of the window: a narrow main lobe around the true frequency and a series of side lobes. The side lobes of the rectangular window are high (the first is only 13 dB below the main lobe) and decay slowly; tapered windows have much lower side lobes, at the cost of a wider main lobe.

A 10 Hz and a 10.5 Hz sine wave (amplitude 2 μV), 1 s at 256 Hz: the 10 Hz wave completes 10 whole cycles, the 10.5 Hz wave does not:

fs = 256
t = collect(0:(1 / fs):(1 - 1 / fs))
rect = ones(length(t))
hann = generate_window(:hann, length(t))
fig = Figure(size = (900, 600))
for (idx, f0) in enumerate([10, 10.5])
    s = 2 .* sin.(2π * f0 .* t)
    ax1 = Axis(fig[idx, 1];
               ax_style_small...,
               title = "$f0 Hz: signal",
               xlabel = "Time [s]",
               ylabel = "Amplitude [μV]")
    lines!(ax1, t, s; color = :gray)
    lines!(ax1, t, s .* hann; color = :red)
    ax2 = Axis(fig[idx, 2];
               ax_style_small...,
               title = "$f0 Hz: amplitude spectrum",
               xlabel = "Frequency [Hz]",
               ylabel = "Amplitude [dB]")
    for (w, c, l) in [(rect, :gray, "no taper"), (hann, :red, "Hann")]
        f, a = amp_spectrum(s, w, fs)
        lines!(ax2, f, 20 .* log10.(a .+ 1e-12); color = c, label = l)
        # values at the frequencies of the FFT of the 1-s segment (1 Hz)
        f1, a1 = amp_spectrum(s, w, fs; n_fft = length(s))
        scatter!(ax2, f1, 20 .* log10.(a1 .+ 1e-12); color = c, markersize = 6)
    end
    xlims!(ax2, 0, 30)
    ylims!(ax2, -80, 15)
    axislegend(ax2; position = :rt)
end
fig

Lines: the continuous spectrum (computed with zero padding), dots: the values at the frequencies computed by the FFT of the 1-s segment (every 1 Hz). Without a taper, the 10 Hz wave appears in a single frequency bin: the other bins fall exactly on the zeros between the side lobes. The 10.5 Hz wave falls between two bins, and without a taper its power leaks to all frequencies, decaying slowly (bins far from 10.5 Hz are only 25-40 dB below the peak). With the Hann window, the leakage drops quickly away from the peak, but the main lobe is wider: the peak covers more bins.

Real EEG never contains only whole numbers of cycles, so leakage is always present: power from strong rhythms (alpha) or from large slow activity leaks into weaker neighboring frequencies, and line noise leaks over the whole spectrum.

Window Functions

generate_window() returns windows of a given length:

n = 256
windows = [:hann, :bh, :bn, :nutall, :flat, :bohman, :triangle, :exp]
fig = Figure(size = fig_size)
ax = Axis(fig[1, 1];
          ax_style...,
          title = "Windows",
          xlabel = "Sample",
          ylabel = "Weight")
# 10 distinct colors
colors = cgrad(:tab10, 10; categorical = true)
lines!(ax, 1:n, ones(n); color = :black, label = "rect")
for (idx, w) in enumerate(windows)
    lines!(ax, 1:n, generate_window(w, n);
           color = colors[idx],
           label = String(w))
end
Legend(fig[1, 2], ax)
fig

Properties of the windows (measured from their spectra, computed below):

# properties of the window w from its zero-padded spectrum
function window_props(w; pad = 64)
    n = length(w)
    W = abs.(fft(vcat(w, zeros(n * (pad - 1)))))[1:(n * pad ÷ 2)]
    db = 20 .* log10.(W ./ W[1] .+ 1e-300)
    # end of the main lobe: first local minimum
    i = findfirst(k -> db[k] < db[k - 1] && db[k] <= db[k + 1],
                  2:(length(db) - 1)) + 1
    side = maximum(db[i:end])
    # full width of the main lobe at -3 dB, in frequency bins (1/n)
    bw = 2 * findfirst(<(-3), db) / pad
    # attenuation of a frequency half-way between two bins
    scallop = db[pad ÷ 2 + 1]
    return side, bw, mean(w), scallop
end
println(rpad("window", 10), rpad("side lobe", 11), rpad("-3 dB width", 13),
        rpad("gain", 7), "scalloping")
for (name, w) in [(:rect, ones(256));
                  [(w, generate_window(w, 256)) for w in windows]]
    side, bw, cg, scallop = window_props(w)
    println(rpad(name, 10),
            rpad(string(round(side; digits = 1), " dB"), 11),
            rpad(string(round(bw; digits = 2), " bins"), 13),
            rpad(round(cg; digits = 2), 7),
            round(scallop; digits = 2), " dB")
end

window    side lobe  -3 dB width  gain   scalloping
rect      -13.3 dB   0.94 bins    1.0    -3.92 dB
hann      -31.5 dB   1.5 bins     0.5    -1.41 dB
bh        -92.0 dB   1.94 bins    0.36   -0.82 dB
bn        -97.9 dB   1.94 bins    0.36   -0.84 dB
nutall    -93.3 dB   1.97 bins    0.35   -0.81 dB
flat      -92.5 dB   3.78 bins    0.21   -0.01 dB
bohman    -46.0 dB   1.75 bins    0.4    -1.02 dB
triangle  -26.5 dB   1.31 bins    0.5    -1.84 dB
exp       -5.1 dB    2.0 bins     0.04   -0.93 dB

side lobe: the highest side lobe relative to the main lobe; lower values mean less leakage to distant frequencies,
−3 dB width: the width of the main lobe (in frequency bins, i.e. multiples of \(1/T\), where \(T\) is the segment length); wider main lobes mean lower frequency resolution: two close peaks merge,
gain (coherent gain): the mean of the window; the amplitude of a tapered sine wave is multiplied by it,
scalloping: the loss of amplitude of a frequency half-way between two bins (the worst case).

Window	Use
rectangular (no taper)	frequencies with exactly whole numbers of cycles; maximum resolution
`:hann`	general-purpose default for EEG spectra (used by `psd()`, `spectrogram()`)
`:bh`, `:bn`, `:nutall`	very low leakage (−92 to −98 dB): weak peaks next to strong ones, line noise
`:flat`	accurate amplitude of peaks between bins (scalloping ≈ 0 dB); wide main lobe
`:bohman`	low leakage (−46 dB) with a moderate main lobe
`:triangle`	simple, moderate leakage
`:exp`	poor taper (high side lobes); not recommended for spectral analysis

Other windows (Hamming, Kaiser, Tukey, Gaussian, cosine, DPSS) are available in DSP.jl and can be used in the same way.

Amplitude correction

A taper reduces the amplitude of the signal, so the spectrum must be divided by the gain of the window (the function amp_spectrum() used here divides by sum(w), i.e. by the gain times the number of samples). Without the correction, the 10 Hz sine wave of 2 μV would appear with an amplitude of about 1 μV after the Hann window (gain 0.5). For the power, the correction factor is the mean of the squared window.

Tapering NEURO Objects

taper() multiplies the selected channels by a window t, which must have the length of one epoch; each epoch is tapered separately (taper!() modifies the object in place):

eeg = import_edf("files/eeg/eeg.edf")
NeuroAnalyzer.filter!(eeg;
                      ch = "all",
                      fprototype = :butterworth,
                      ftype = :hp,
                      cutoff = 0.5,
                      order = 4)
NeuroAnalyzer.filter!(eeg;
                      ch = "all",
                      fprototype = :fir,
                      ftype = :lp,
                      cutoff = 35,
                      bw = 10)
e2 = epoch(crop(eeg; seg = (100, 160)); ep_len = 2)
e2_t = taper(e2; ch = "eeg", t = generate_window(:hann, epoch_len(e2)))

NeuroAnalyzer.plot(e2_t; ch = ["Fp1", "C3", "O1"], ep = 1, gui = false)

Each epoch now starts and ends at zero. Only the selected channels are changed (here the EEG channels; EOG and ECG are not).

Tapering and psd()

psd(), spectrogram() and the other spectral functions already apply a Hann window by default (w = true). After tapering with taper(), set w = false, otherwise the signal is tapered twice:

e2_bh = taper(e2; ch = "eeg", t = generate_window(:bh, epoch_len(e2)))
p_bh = psd(e2_bh; ch = "O1", method = :fft, w = false)
p_hann = psd(e2; ch = "O1", method = :fft)

When to taper

taper epochs (short segments) immediately before spectral analysis; do not edit or analyze the tapered signal in the time domain (the attenuated edges distort ERPs, amplitudes and connectivity),
do not taper continuous recordings before filtering or epoching; filter first, then split into epochs, reject bad epochs and taper,
Welch’s method (method = :welch, default of psd()) tapers each of its short windows; multitaper (method = :mt) uses several orthogonal (Slepian) tapers and averages their spectra, which reduces the variance of the estimate.