Introduction to Signal Theory

Initialize NeuroAnalyzer

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

A practical guide to signals, systems, and the mathematical foundations of neural signal processing.

Signals

A signal is a quantity that varies over time (or space) and carries information, e.g. the voltage recorded by an EEG electrode.

A continuous-time signal $x(t)$ is defined at every instant $t$ (e.g. the electrical potential on the scalp).
A discrete-time signal $x[n]$ is defined only at integer sample indices $n$; it is obtained by sampling a continuous signal every $T_s$ seconds: $x[n] = x(n T_s)$.
A digital signal is a discrete-time signal whose values are also quantized (stored with a finite number of levels by the analog-to-digital converter).

Below, a continuous sine (line) and its samples (stems):

t = linspace(0, 1, 1000)       # "continuous" time
ts = 0:(1 / 20):1              # sampled every 50 ms (20 Hz)
fig, ax = lines(t, sin.(2π * 2 .* t); axis = (xlabel = "Time [s]", ylabel = "Amplitude"))
stem!(ax, ts, sin.(2π * 2 .* ts))
fig

Time series vs. signal

\[ \text{time series} = \text{signal} + \text{noise} \]

A time series is the recorded sequence of samples; the signal is the component of interest within it, and noise is everything else (other brain activity, artifacts, amplifier noise). In EEG, the time series is the voltage (μV) recorded over time (s or ms).

Key Terminology

Period ($T$): duration of one cycle of an oscillation (s).

Frequency ($f$): number of cycles per second (Hz):

\[ f = 1 / T \]

For a sine wave, the frequency is also half the number of zero-crossings per second (a sine crosses zero twice in each cycle).

Sampling rate ($f_s$): number of samples per second (Hz).

Sampling interval ($T_s$): time between consecutive samples: $T_s = 1 / f_s$.

Nyquist frequency: the highest frequency that can be represented at a given sampling rate: $f_s / 2$ (see Sampling and aliasing).

Fundamental frequency ($f_0$): the frequency of a periodic signal, $1 / T$; the lowest frequency of its harmonic components.

Harmonics: integer multiples of the fundamental frequency ($2 f_0, 3 f_0, \dots$).

Octave: a doubling of frequency (e.g. 4 → 8 Hz).

Dominant (peak) frequency: the frequency component with the highest amplitude (power), e.g. the alpha peak frequency.

Sampling and Aliasing

The sampling theorem (Nyquist–Shannon) states that a signal containing no frequencies above $f_{max}$ can be reconstructed exactly from its samples if the sampling rate is higher than twice this frequency:

\[ f_s > 2 f_{max} \]

A component above the Nyquist frequency ($f_s / 2$) is not lost, but appears at a wrong, lower frequency: it is aliased. A component of frequency $f$ appears at $|f - k f_s|$ for the integer $k$ that brings it into $[0, f_s / 2]$.

A 40 Hz sine sampled at 50 Hz (Nyquist frequency: 25 Hz) is indistinguishable from a 10 Hz sine:

t = linspace(0, 0.5, 2000)
ts = 0:(1 / 50):0.5
fig, ax = lines(t, sin.(2π * 40 .* t); color = :grey, label = "40 Hz",
                axis = (xlabel = "Time [s]", ylabel = "Amplitude"))
lines!(ax, t, -sin.(2π * 10 .* t); color = :red, label = "aliased: 10 Hz")
scatter!(ax, ts, sin.(2π * 40 .* ts); color = :black, label = "samples (50 Hz)")
axislegend(ax; position = :rb)
fig

Therefore, the analog signal is low-pass filtered by the amplifier (anti-aliasing filter) before it is digitized, and a signal must be low-pass filtered before it is downsampled (resample() does this).

Systems

A system transforms an input signal into an output signal:

  x  ──►  H  ──►  y

For a continuous-time system: $y(t) = H\{x(t)\}$; for a discrete-time system: $y[n] = H\{x[n]\}$. Filters, amplifiers and re-referencing are examples of systems.

Linear Systems

A system is linear if it obeys the superposition principle, i.e. both:

Additivity: $H\{x_1[n] + x_2[n]\} = H\{x_1[n]\} + H\{x_2[n]\}$
Homogeneity (scaling): $H\{a \, x[n]\} = a \, H\{x[n]\}$

For example, $y[n] = 2x[n] + x[n - 1]$ is linear; $y[n] = x[n]^2$ and $y[n] = x[n] + 1$ are not.

Time-Invariant Systems

A time-invariant system behaves identically regardless of when the input arrives: a shift of the input by $k$ samples produces the same shift of the output.

\[ H\{x[n]\} = y[n] \Rightarrow H\{x[n - k]\} = y[n - k] \]

For example, $y[n] = x[n] - x[n - 1]$ is time-invariant; $y[n] = n \, x[n]$ is not.

LTI Systems

A system that is both linear and time-invariant is an LTI system. LTI systems are completely described by their impulse response (see below); all filters used in EEG analysis are LTI systems.

Memoryless Systems

A system is memoryless when its output at time $n$ depends only on its input at that same instant.

Memoryless system:

\[ y[n] = x[n]^2 \]

System with memory:

\[ y[n] = x[n - 1] \]

Causal Systems

A causal system’s output depends only on present and past inputs, never on future ones.

Causal system:

\[ y[n] = x[n] - 2 x[n - 1] \]

Non-causal system:

\[ y[n] = x[n + 3] \]

A causal system cannot anticipate future inputs, so only causal systems can work in real time. For example, a 3-sample moving average of samples $n-1$, $n$ and $n+1$ is non-causal: it requires a future sample.

Recorded data can be processed with non-causal systems: zero-phase filtering (filtering forward and backward, the default in NeuroAnalyzer.filter()) is non-causal, which is why it does not shift the signal in time.

Stable Systems

A system is (BIBO) stable if every bounded input produces a bounded output. Recursive (IIR) filters can be unstable; FIR filters are always stable.

Feedback

A feedback loop routes part of a system’s output back to its input (e.g. recursive IIR filters). Feedback can stabilize a system, but it can equally introduce instability if the loop gain is too large.

Special Signals

Unit Impulse

The discrete-time unit impulse (Kronecker delta) $\delta[n]$ is a single sample of value 1 at $n = 0$:

\[ \delta[n] = \begin{cases} 1 & \text{if $n = 0$}\\ 0 & \text{otherwise} \end{cases} \]

Its continuous-time counterpart is the Dirac delta $\delta(t)$: zero everywhere except $t = 0$, with a unit area.

Any discrete signal is a sum of scaled and shifted impulses:

\[ x[n] = \sum_{k} x[k] \, \delta[n - k] \]

Unit Step

The unit step $u[n]$ switches from 0 to 1 at $n = 0$:

\[ u[n] = \begin{cases} 1 & \text{if $n \geq 0$}\\ 0 & \text{otherwise} \end{cases} \]

It is a running sum of impulses (and the impulse is the difference of two steps, $\delta[n] = u[n] - u[n - 1]$):

\[ u[n] = \sum_{k=0}^{\infty} \delta[n - k] = \sum_{k=-\infty}^{n} \delta[k] \]

Impulse Response and Convolution

The impulse response $h[n]$ is the output of a system when the input is a unit impulse. Because any input is a sum of shifted impulses, the output of an LTI system is the sum of the correspondingly scaled and shifted impulse responses: the convolution of the input with the impulse response:

\[ y[n] = (x * h)[n] = \sum_{k} x[k] \, h[n - k] \]

The impulse response therefore completely characterizes an LTI system. For example, a 5-sample moving average has the impulse response $h = [0.2, 0.2, 0.2, 0.2, 0.2]$:

h = fill(0.2, 5)
x = zeros(30)
x[5] = 1                                  # impulse
x[15:30] .= 1                             # step
y = NeuroAnalyzer.DSP.conv(x, h)[1:30]
fig, ax = stem(0:29, x; label = "input",
               axis = (xlabel = "Sample", ylabel = "Amplitude"))
stem!(ax, 0:29, y; color = :red, label = "output")
axislegend(ax; position = :lt)
fig

The impulse at sample 4 is spread over 5 samples (the impulse response) and the step rises gradually (the step response, the running sum of the impulse response).

Elementary Signals

t = -1:0.01:1
us = Float64.(t .>= 0)                    # unit step
imp = Float64.(t .== 0)                   # impulse
ramp = t .* us                            # ramp
quadratic = t .^ 2 .* us                  # quadratic
expd = exp.(-5 .* t) .* us                # decaying exponential (t ≥ 0)
hwr = sin.(2π * 2 .* t) .* (sin.(2π * 2 .* t) .> 0)  # half-wave rectified sine
fwr = abs.(sin.(2π * 2 .* t))             # full-wave rectified sine

fig = Figure(size = (900, 500))
for (i, (y, name)) in enumerate(zip([us, imp, ramp, quadratic, expd, hwr, fwr],
                                    ["unit step", "impulse", "ramp", "quadratic",
                                     "decaying exponential", "half-wave rectified", "full-wave rectified"]))
    lines!(Axis(fig[(i - 1) ÷ 4 + 1, (i - 1) % 4 + 1]; title = name), t, y)
end
fig

A step can be used to switch (gate) another signal on and off:

gated = sin.(2π * 5 .* t) .* us
NeuroAnalyzer.plot_signal(collect(t), gated; title = "Gated 5 Hz sine")

Signal Classification

Signals can be classified as:

Continuous-time or discrete-time (and analog or digital)
Periodic or aperiodic
Even, odd or neither
Energy or power signals
Real or complex
Deterministic or random (stochastic)
Stationary or non-stationary
Causal (zero for $t < 0$) or non-causal

Periodic signal: repeats with a period $T$: $x(t + T) = x(t)$ for all $t$ (the smallest such $T$ is the fundamental period). A discrete-time sinusoid $\sin(2\pi f n / f_s)$ is periodic only if $f / f_s$ is a rational number.

Even signal: symmetric about $t = 0$: $x(t) = x(-t)$ (e.g. cosine).

Odd signal: antisymmetric about $t = 0$: $x(t) = -x(-t)$ (e.g. sine). Every signal is the sum of an even and an odd part: $x_e(t) = \frac{x(t) + x(-t)}{2}$, $x_o(t) = \frac{x(t) - x(-t)}{2}$.

Deterministic signal: its values are completely determined by a formula (e.g. a sine). Random signal: its values can be described only statistically (e.g. noise, the background EEG).

Stationary signal: its statistical properties (mean, variance, spectrum) do not change over time. EEG is non-stationary in the long term (sleep stages, events), but is often treated as stationary within short segments (a few seconds), which is the assumption behind spectral analysis of epochs.

Operations on Signals

Amplification: $k \, x(t)$ for $|k| > 1$ - increases the amplitude.

Attenuation: $k \, x(t)$ for $|k| < 1$ - decreases the amplitude.

Time shift: $x(t - t_0)$ - delays the signal (shift to the right) for $t_0 > 0$; $x(t + t_0)$ - advances it (shift to the left).

Time scaling: $x(k t)$ compresses the signal in time for $k > 1$ (it becomes faster) and stretches it for $0 < k < 1$ (it becomes slower).

Time reversal: $x(-t)$ mirrors the signal about $t = 0$.

Addition and multiplication: two signals are added, $a(t) + b(t)$, or multiplied, $a(t) \times b(t)$, sample by sample (multiplication by a sine is used in amplitude modulation and in the Fourier transform).

t = 0:0.001:2
wave(t) = sin(2π * 1 * t)
fig, ax = lines(t, wave.(t); color = :black, label = "x(t)",
                axis = (xlabel = "Time [s]", ylabel = "Amplitude"))
lines!(ax, t, wave.(2 .* t); color = :red, label = "x(2t): compressed")
lines!(ax, t, wave.(0.5 .* t); color = :blue, label = "x(0.5t): stretched")
lines!(ax, t, wave.(t .- 0.25); color = :green, linestyle = :dash, label = "x(t - 0.25): delayed")
axislegend(ax; position = :rb)
fig

Sine Wave

The general formula of the sine wave is:

\[ y(t) = a \sin(2\pi f t + \theta) \]

where:
$a$ is the amplitude,
$f$ is the frequency (Hz),
$t$ is the time (s),
$\theta$ is the phase (rad); $\theta = 0$: sine, $\theta = \pi/2$: cosine.

Any oscillation is characterized by three parameters:

Frequency: how fast does it oscillate?
Amplitude (or power, the squared amplitude): how strong is it?
Phase: where in its cycle is the oscillation at $t = 0$?

Sine and cosine of the same frequency are orthogonal: their product integrates to zero over a full period. Both are related to complex exponentials by Euler’s formula:

\[ \sin(t) = \frac{e^{i t} - e^{-i t}}{2i} = \Im(e^{i t}), \quad \cos(t) = \frac{e^{i t} + e^{-i t}}{2} = \Re(e^{i t}) \]

Generating a Sine Wave

generate_sine(f, t, a, p) generates a sine of frequency f (Hz), amplitude a and phase p (degrees) for the time vector t:

fs = 256
t = collect(0:(1 / fs):2)            # 2 seconds
s = generate_sine(3, t, 10, 90)      # 3 Hz, amplitude 10, phase 90° (cosine)
# or: s = @. 10 * sin(2π * 3 * t + deg2rad(90))
NeuroAnalyzer.plot_signal(t, s)

Superposition of Sine Waves

Signals composed of several frequencies are sums of sine waves (this is the basis of the Fourier transform):

fs = 256
t = collect(0:(1 / fs):2)
a = [10, 5, 2]          # amplitudes
f = [2, 5, 10]          # frequencies (Hz)
p = [0, 180, 90]        # phases (degrees)
s = zeros(length(t))
for idx in eachindex(a)
    s .+= generate_sine(f[idx], t, a[idx], p[idx])
end
NeuroAnalyzer.plot_signal(t, s)

Complex Sine Wave

Euler’s formula links trigonometric functions and complex exponentials:

\[ e^{i k} = \cos(k) + i \sin(k) \]

Substituting the argument of the sine wave gives the complex sine wave:

\[ y(t) = a e^{i 2 \pi f t} = a \underbrace{\cos(2 \pi f t)}_{\text{real part}} + i \, a \underbrace{\sin(2 \pi f t)}_{\text{imaginary part}} \]

The real and imaginary parts are a cosine and a sine (phase difference: $\pi / 2$); the complex sine wave is a point rotating on a circle of radius $a$ in the complex plane. Complex sine waves are used in the Fourier transform and in complex wavelets, because the magnitude and the angle of the complex result give the amplitude and the phase directly.

t = collect(0:0.001:2)
s = generate_csine(2, t, 1)
# or: s = @. exp(1im * 2π * 2 * t)

fig, ax = lines(t, real.(s); label = "real part (cosine)",
                axis = (xlabel = "Time [s]", ylabel = "Amplitude"))
lines!(ax, t, imag.(s); color = :red, label = "imaginary part (sine)")
axislegend(ax; position = :rb)
fig

lines(t, real.(s), imag.(s);
      axis = (type = Axis3, xlabel = "Time [s]", ylabel = "Real", zlabel = "Imaginary"))

Amplitude

Amplitude measures the magnitude of the oscillation:

Code

t = range(0, 1, 501)
s = sin.(2π .* t)
rms = 1 / sqrt(2)
fig, ax = lines(t, s; color = :red, linewidth = 2,
                axis = (xlabel = "Time", ylabel = "Amplitude",
                        yticks = ([-1, 0, rms, 1], ["-û", "0", "û/√2", "û"])))
hlines!(ax, 0; color = :black)
lines!(ax, [0, 0.75], [1, 1]; color = :gray, linestyle = :dash)
lines!(ax, [0, 0.25], [rms, rms]; color = :gray, linestyle = :dash)
lines!(ax, [1, 1], [0, -1.25]; color = :gray, linestyle = :dash)
# double-headed arrows: ③ RMS, ① peak, ② peak-to-peak, ④ period
arrows2d!(ax, [Point2f(0.15, 0), Point2f(0.55, 0), Point2f(0.75, -1), Point2f(0, -1.25)],
          [Point2f(0.15, rms), Point2f(0.55, 1), Point2f(0.75, 1), Point2f(1, -1.25)];
          argmode = :endpoint, color = :gray, tipwidth = 10, tiplength = 10,
          tail = Point2f[(0, 0), (1, -0.5), (1, 0.5)], tailwidth = 10, taillength = 10)
text!(ax, [0.17, 0.57, 0.77, 0.52], [rms / 2, 0.5, 0.5, -1.18];
      text = ["③", "①", "②", "④"], fontsize = 24, color = :blue,
      align = (:left, :center))
fig

Peak amplitude ($\hat{u}$): the maximum absolute deviation of the signal from its baseline (or zero).
Peak-to-peak amplitude ($2\hat{u}$ for a sine): the difference between the maximum and the minimum of the signal.
RMS amplitude (root mean square): $\hat{u} / \sqrt{2}$ for a sine.
Wave period ($T$): a measure of time, not amplitude.

Peak amplitude is used for transient events (e.g. ERP components, epileptic spikes). Peak-to-peak amplitude is used for oscillations (e.g. alpha waves) and for artifact rejection (e.g. epochs exceeding ±100 μV). RMS amplitude summarizes the magnitude of a varying signal over a time window, accounting for all samples; for a discrete signal $s$ of $N$ samples:

\[ \text{RMS} = \sqrt{\frac{1}{N} \sum_{i=1}^{N} s_i^2} \]

For a sine wave of amplitude 10:

s = generate_sine(10, collect(0:(1 / 256):(1 - 1 / 256)), 10)
(peak = maximum(abs.(s)),
 peak_to_peak = maximum(s) - minimum(s),
 rms = sqrt(mean(s .^ 2)),
 expected_rms = 10 / sqrt(2))

(peak = 10.0, peak_to_peak = 20.0, rms = 7.071067811865475, expected_rms = 7.071067811865475)

Phase and Time Shift

Shifting a sine wave in time is equivalent to changing its phase: replacing $t$ with $t - t_0$ (a delay by $t_0$) gives

\[ a \sin(2\pi f (t - t_0)) = a \sin(2\pi f t - 2\pi f t_0) \]

so a delay by $t_0$ decreases the phase by $2\pi f t_0$, and an advance increases it. The same time shift corresponds to a larger phase shift at higher frequencies.

In Fourier analysis the phase is usually measured relative to a cosine: phase 0 means that a peak of the wave is at $t = 0$.

t = 0:0.01:4
s1 = sin.(2π * 0.5 .* t)
s2 = sin.(2π * 0.5 .* (t .+ 0.5))      # advanced: phase + π/2
s3 = sin.(2π * 0.5 .* (t .- 0.5))      # delayed: phase - π/2
fig, ax = lines(t, s1; color = :black, label = "original",
                axis = (xlabel = "Time [s]", ylabel = "Amplitude"))
lines!(ax, t, s2; color = :red, label = "t + 0.5 s")
lines!(ax, t, s3; color = :blue, label = "t - 0.5 s")
axislegend(ax; position = :rt)
fig

Noise

Noise is characterized by its amplitude distribution and by its spectrum (how its power is distributed over frequencies).

White Noise

White noise has equal power at all frequencies (a flat spectrum). Its samples are independent and can have any distribution, e.g. normal (Gaussian) or uniform.

Gaussian white noise:

fs = 100
t = collect(0:(1 / fs):10)
gwn = generate_noise(length(t); type = :whiten)
NeuroAnalyzer.plot_signal(t, gwn; title = "Gaussian White Noise")

plot_histogram(gwn;
               title = "Amplitude Distribution",
               xlabel = "Amplitude",
               ylabel = "Count",
               draw_median = false)

X = ftransform(gwn; db = true)
f, _ = freqs(length(gwn), fs)
plot_psd(f, X.p; flim = (1, 45), xlabel = "Frequency [Hz]", ylabel = "Power [dB]")

Uniform white noise:

uwn = generate_noise(length(t); type = :whiteu)
plot_signal(t, uwn; title = "Uniform White Noise")

plot_histogram(uwn;
               title = "Amplitude Distribution",
               xlabel = "Amplitude",
               ylabel = "Count",
               draw_median = false)

X = ftransform(uwn; db = true)
plot_psd(f, X.p; flim = (1, 45), xlabel = "Frequency [Hz]", ylabel = "Power [dB]")

Both spectra are flat (apart from random fluctuations); only the amplitude distributions differ.

Pink Noise

Pink noise ($1/f$ noise) has a power that decreases with frequency: $P(f) \propto 1 / f$, i.e. the power falls by 3 dB per octave. Signals with $1/f$-like spectra are common in biology; the background (aperiodic) EEG spectrum follows $1 / f^{\beta}$ with $\beta$ typically between 1 and 2.

pn = generate_noise(length(t); type = :pink)
plot_signal(t, pn; title = "Pink Noise")

plot_histogram(pn;
               title = "Amplitude Distribution",
               xlabel = "Amplitude",
               ylabel = "Count",
               draw_median = false)

On logarithmic axes the $1/f$ spectrum is a straight line:

X = ftransform(pn; db = true)
plot_psd(f, X.p; flim = (1, 45), frq = :log, xlabel = "Frequency [Hz]", ylabel = "Power [dB]")

Brown (red) noise ($1/f^2$) is the running sum (integral) of white noise, a random walk.

Resolution

Amplitude resolution of an analog-to-digital converter (ADC) with $b$ bits: $2^b$ distinct levels (a 16-bit ADC: 65 536 levels, a 24-bit ADC: 16 777 216 levels).

Quantization: each sample is rounded to the nearest level, which adds quantization noise (at most half a level).

Dynamic range: about $6.02 \times b$ dB (16 bits: ~96 dB, 24 bits: ~144 dB).

Smallest step (least significant bit): input range ÷ number of levels. Example: an input range of ±100 mV (200 mV) digitized with 16 bits gives $200 / 65536 \approx 0.003$ mV = 3 μV; with 24 bits, 0.012 μV. Modern EEG amplifiers use 24-bit ADCs.

Computations use floating-point numbers, whose relative precision is $\approx 1.2 \times 10^{-7}$ for 32-bit (eps(Float32)) and $\approx 2.2 \times 10^{-16}$ for 64-bit numbers (eps(Float64), used by NeuroAnalyzer).

Temporal resolution: the sampling interval, $1 / f_s$. Sampling rates of 250–2000 Hz cover most EEG analyses; 1000 Hz is convenient (1 sample = 1 ms). Higher rates (e.g. 5000 Hz) are used for brainstem evoked potentials and for simultaneous EEG–fMRI (to capture the MRI gradient artifact).

Frequency resolution: the spacing of the frequency bins of the Fourier transform, determined by the duration $D$ of the analyzed segment:

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

where $N$ is the number of samples. A 2-second segment gives 0.5 Hz bins, regardless of the sampling rate.

The Sinc Function

The sinc function is the Fourier transform of a rectangular function (and vice versa); it is the impulse response of an ideal low-pass filter and the basis of ideal reconstruction of a signal from its samples.

\[ \mathrm{sinc}(x) = \frac{\sin(\pi x)}{\pi x} \]

(the normalized sinc, as sinc() in Julia; the unnormalized form is $\sin(x) / x$).

As a function of time, with frequency $f$ and the peak at time $m$ (generate_sinc()):

\[ s(t) = \frac{\sin(2 \pi f (t - m))}{\pi (t - m)} \]

Real Sinc

t = collect(-1:0.001:1)
s = generate_sinc(t; f = 20)
plot_signal(t, s)

Complex Sinc

$\sin(z) / z$ for complex $z$ (magnitude):

xs = LinRange(-10, 10, 100)
ys = LinRange(-4, 4, 100)
zs = [abs(sin(complex(x, y)) / complex(x, y)) for x in xs, y in ys]
surface(xs, ys, log10.(zs); axis = (type = Axis3, xlabel = "Re", ylabel = "Im", zlabel = "log₁₀|sinc|"))

The Gaussian Function

The Gaussian function is a symmetric, bell-shaped curve:

\[ g(x) = a \, e^{-\frac{(x - b)^2}{2c^2}} \]

where:
$a$ is the peak height,
$b$ is the position of the peak,
$c$ is the standard deviation (the width).

Gaussians are used in EEG analysis as:

windows in time-frequency analysis: the Gabor transform (a short-time Fourier transform with a Gaussian window) and the Morlet wavelet (a complex sine wave multiplied by a Gaussian)
smoothing kernels (Gaussian filters)
models of the shape of ERP components and spectral peaks (e.g. in FOOOF)
the probability density of the normal distribution

The Fourier transform of a Gaussian is a Gaussian: a narrow Gaussian in time is wide in frequency and vice versa (the time-frequency trade-off).

Time-Domain Gaussian

Centered at 0 ($b = 0$), with the width parameter $w = 2c^2$:

\[ g(t) = a \, e^{-\frac{t^2}{w}} \]

Gaussian with FWHM

The full width at half maximum (FWHM, $h$) is often more intuitive than the standard deviation:

\[ g(t) = a \, e^{-4 \ln(2) \frac{(t - b)^2}{h^2}}, \qquad h = 2 \sqrt{2 \ln 2} \, c \approx 2.355 \, c \]

Frequency-Domain Gaussian

Centered at the frequency $f_c$, with the standard deviation $c$ derived from the FWHM $h$ (Hz):

\[ g(f) = e^{-\frac{(f - f_c)^2}{2c^2}}, \qquad c = \frac{h}{2 \sqrt{2 \ln 2}} \]

fs = 256
g = generate_gaussian(fs, 10, 0.5; ncyc = 5)       # σ = 5 / (2π × 10 Hz)
tg = collect(-0.5:(1 / fs):0.5)
plot_signal(tg, g; title = "Gaussian (σ = $(round(5 / (2π * 10); digits = 3)) s)")

Common Waveforms

fs = 100
t = collect(0:(1 / fs):10)

Square Wave

A pulse wave alternates between two levels; its duty cycle is the fraction of the period spent at the high level (pulse width ÷ period). A square wave has a duty cycle of 50%:

s = Float64.(mod.(t, 2) .< 1)      # period 2 s, duty cycle 50%
plot_signal(t, s)

Triangle Wave

s = @. abs(mod(t, 2) - 1)          # period 2 s
plot_signal(t, s)

Chirp (Frequency Sweep)

A chirp is a non-stationary signal whose instantaneous frequency changes over time. In a linear chirp the frequency increases from $f_1$ to $f_2$ over the duration $D$:

\[ s(t) = \sin\left(2\pi \left(f_1 t + \frac{f_2 - f_1}{2D} t^2\right)\right) \]

f1, f2 = 2, 10
D = t[end]
s = @. sin(2π * (f1 * t + (f2 - f1) / (2D) * t^2))
plot_signal(t, s)

Linearly Increasing Amplitude

a = linspace(1, 10, length(t))
s = @. a * sin(2π * 3 * t)
plot_signal(t, s)

Abrupt Frequency Changes

A signal composed of non-overlapping segments of different frequencies and amplitudes:

a = [2, 5, 10]
f = [5, 10, 20]
edges = round.(Int64, range(0, length(t); length = length(a) + 1))
s = zeros(length(t))
for idx in eachindex(a)
    seg = (edges[idx] + 1):edges[idx + 1]
    s[seg] = a[idx] .* sin.(2π * f[idx] .* t[seg])
end
plot_signal(t, s)

Energy and Power

Instantaneous power is the squared signal value:

\[ p(t) = |x(t)|^2 \]

Energy is the instantaneous power integrated over time (for a discrete signal: summed over samples):

\[ E = \int_{-\infty}^{\infty} |x(t)|^2 \, \mathrm{d}t, \qquad E = \sum_{n} |x[n]|^2 \]

Average power is the energy per unit time:

\[ P = \lim_{T \to \infty} \frac{1}{T} \int_{-T/2}^{T/2} |x(t)|^2 \, \mathrm{d}t, \qquad P = \frac{1}{N} \sum_{n=1}^{N} |x[n]|^2 \]

An energy signal has finite energy (and zero average power), e.g. a transient pulse; a power signal has finite, non-zero average power (and infinite energy), e.g. a sine wave or a stationary EEG rhythm. The average power is the squared RMS amplitude; for a sine of amplitude $a$ it is $a^2 / 2$.

Parseval’s theorem: the energy of a signal is the same in the time and in the frequency domain, which is why the power spectrum distributes the signal power over frequencies:

x = generate_sine(10, collect(0:(1 / 256):(2 - 1 / 256)), 3)
X = FFTW.fft(x)
(energy_time = sum(abs2, x),
 energy_frequency = sum(abs2, X) / length(x),
 average_power = mean(abs2, x),
 expected_power = 3^2 / 2)

(energy_time = 2304.000000000001, energy_frequency = 2304.0, average_power = 4.500000000000002, expected_power = 4.5)

Decibels (dB)

Decibels express ratios on a logarithmic scale:

\[ \text{Power (dB)} = 10 \log_{10}(P_1 / P_2) \]

\[ \text{Amplitude (dB)} = 20 \log_{10}(A_1 / A_2) \]

The factor of 20 for amplitudes arises because power is proportional to the squared amplitude: $10 \log_{10}(A^2) = 20 \log_{10}(A)$. Doubling the power is +3 dB, doubling the amplitude is +6 dB; +10 dB is a 10× increase of power, +20 dB a 10× increase of amplitude.

Conversion between power and dB:

pow2db(x) = 10 .* log10.(x)
db2pow(x) = 10 .^ (x ./ 10)
(pow2db([1, 2, 10, 100]), db2pow([0, 3, 10, 20]))

([0.0, 3.010299956639812, 10.0, 20.0], [1.0, 1.9952623149688795, 10.0, 100.0])

(DSP.pow2db(), DSP.db2pow(), DSP.amp2db() and DSP.db2amp() are also available.)