Sampling, Nyquist and aliasing

Why this matters

Downsampling is the first processing step in many pipelines, because it makes everything after it faster, and it is the step most often done wrong. Take every fifth sample of a 500 Hz recording and you have not thrown away the content above 50 Hz; you have folded it into the band you kept, where it is indistinguishable from brain activity. The same arithmetic explains why a 60 Hz mains harmonic can turn up at 40 Hz, why a 160 Hz recording cannot show you its own line-noise harmonics, and why a headset’s sampling rate says nothing about the frequencies it can record. Twenty minutes with the formula saves years of chasing artifacts that were manufactured by a resampling call.

Concepts

The sampling theorem

A continuous signal that contains no frequencies above B hertz is completely determined by its samples taken at any rate fs greater than 2B; from those samples the continuous signal can be reconstructed exactly (by sinc interpolation). Stated the other way round: a recording sampled at fs can faithfully represent frequencies only up to fs/2, the Nyquist frequency. At 160 Hz (ds-eegbci) the Nyquist frequency is 80 Hz; at 500 Hz (ds-iowapd) it is 250 Hz; at 2500 Hz (ds-lemon raw) it is 1250 Hz.

The condition is strict. A component at the Nyquist frequency is not represented properly (its amplitude depends on its phase relative to the sampling instants), and a real anti-alias filter needs a transition band, so usable content ends somewhat below fs/2 in every practical system.

Aliasing: where a frequency above Nyquist goes

A component above the Nyquist frequency does not vanish when sampled; the samples are exactly consistent with a lower frequency, and that lower frequency is what every later analysis will see. The folded frequency is

f_alias = | f − round(f / fs) · fs |

that is, the distance from f to the nearest integer multiple of the sampling rate. Two worked cases:

• A 140 Hz tone sampled at 200 Hz: the nearest multiple of 200 is 200, so the alias is |140 − 200| = 60 Hz. The samples of a 140 Hz sinusoid and of a 60 Hz sinusoid are identical.
• A 70 Hz component in a 500 Hz recording decimated to 100 Hz without filtering: the nearest multiple of 100 is 100, so it appears at |70 − 100| = 30 Hz, squarely inside the beta band. This is the notebook’s demonstration on ds-iowapd.

Line noise gives a third case that matters in ds-eegbci: at 160 Hz the mains harmonics at 120 Hz and 180 Hz lie above the 80 Hz Nyquist frequency. If they reached the converter unfiltered they would alias to |120 − 160| = 40 Hz and |180 − 160| = 20 Hz. Whether any such aliased residue is present in the data is TODO(confirm) from the recordings themselves (§6 L1.6); the arithmetic tells you where to look.

Anti-alias filtering, in hardware and in software

Because aliasing happens at the moment of sampling, the only defence is to remove the offending content before sampling. In the amplifier this is a hardware low-pass ahead of the converter (L0.3). In software the same rule applies to decimation: to reduce a 500 Hz recording to 100 Hz you must first low-pass filter it below the new Nyquist frequency of 50 Hz, and only then keep every fifth sample. raw.resample() in MNE does both steps (an anti-alias low-pass followed by resampling), and so do scipy.signal.decimate and scipy.signal.resample_poly; taking data[::5] does neither. The site’s DSP library offers both resample (filtered) and decimateNaive (unfiltered) so that the widget can show you the difference on real data.

Resampling versus decimation. Decimation keeps every n-th sample and so only allows integer rate ratios; resampling interpolates to an arbitrary rate. Both need the anti-alias filter when the rate goes down; going up needs no filter (there is nothing above the old Nyquist frequency to fold) but creates no new information either. Resampling also shifts event samples: MNE resamples the stimulus channel and annotations along with the data, but any event table you keep separately must be converted with the same ratio, and rounding can move an event by a sample.

Choosing a sampling rate

The rate is set by the highest frequency you will analyse plus room for the anti-alias filter’s transition band, not by habit and not by the maximum your amplifier offers:

• For ERPs and oscillations up to the beta band, 250 Hz is generous; for gamma and for high-frequency oscillations, higher rates and, above all, an amplifier whose hardware bandwidth actually extends that far.
• Record at a higher rate than you need and downsample properly afterwards; you can always throw information away, never recover it.
• Keep enough rate that filter transition bands (L1.5) have room: a 45 Hz low-pass with a wide transition band needs the Nyquist frequency comfortably above it.

And read the hardware bandwidth before believing the rate. ds-arithmetic is nominally 500 Hz but was hardware low-pass filtered near 30 Hz (the descriptor’s wording is ambiguous, and the catalog reads it as a low-pass); the consumer datasets in the directory output 128 or 256 Hz from hardware with about 43 Hz of bandwidth and built-in mains notches (§10.9). In all of them the band between the hardware ceiling and the Nyquist frequency contains attenuated signal and noise, not brain activity.

The data behind this lesson

• ds-iowapd sub-001, channel Oz, 10 s at 500 Hz (raw, 0.1 Hz online high-pass, 60 Hz mains with harmonics): the widget’s real snippet and the notebook’s decimation subject. CC0.
• ds-eegbci S001 R02, channel O1, 10 s at 160 Hz: the counterpart snippet, to show a recording whose usable band is already narrow. ODC-By 1.0.
• ds-arithmetic (500 Hz, ~30 Hz hardware cut-off, 50 Hz hardware notch): the counter-example that a sampling rate says nothing about bandwidth. ODC-By 1.0.

Explore

Type 140 Hz at a 200 Hz rate and read the alias; then sweep the tone across the Nyquist frequency and watch the alias reflect back. Switch to the real 500 Hz snippet, decimate to 100 Hz with the anti-alias filter off, and find the mains line that has moved: 60 Hz folds to |60 − 100| = 40 Hz. Turn the filter on and it is gone rather than moved.

Practice

The notebook compares raw.resample() with naive decimation on ds-iowapd, shows a 70 Hz component folding to 30 Hz when 500 Hz is decimated to 100 Hz without filtering, and then plots the PSD of ds-arithmetic up to its Nyquist frequency to find the hardware knee.

Exercises