Introduction to Fourier Transform

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

The Fourier transform decomposes a signal into its frequency components. It transforms a time-domain signal into its frequency-domain representation: for each frequency, how strongly it is present in the signal (amplitude or power) and how it is aligned in time (phase).

Fourier Theorem

According to the Fourier theorem, a signal can be represented as a sum of sine waves, each with its own frequency, amplitude and phase. A periodic signal is a sum of sine waves at multiples of its fundamental frequency (a Fourier series); a finite, sampled signal of \(N\) samples is exactly a sum of \(N\) complex sine waves (the discrete Fourier transform).

The Fourier transform of a continuous-time signal \(s(t)\):

\[ X(f) = \int_{-\infty}^{\infty} s(t) e^{-i 2\pi ft} \, \mathrm{d}t \]

The discrete Fourier transform (DFT) of a signal \(s[n]\) of \(N\) samples:

\[ X[k] = \sum_{n=0}^{N - 1} s[n] e^{-i 2\pi kn/N}, \qquad k = 0, 1, \dots, N - 1 \]

The coefficient \(X[k]\) corresponds to the frequency \(f_k = k f_s / N\), where \(f_s\) is the sampling rate.

The fast Fourier transform (FFT) is an efficient algorithm to compute the DFT and its inverse (in \(N \log N\) instead of \(N^2\) operations), exploiting the symmetry and periodicity of the complex exponentials.

How it works

Each DFT coefficient is the dot product of the signal with a complex sine wave \(e^{-i 2\pi kn/N}\) of the frequency \(f_k\): the signal is multiplied by the sine wave sample by sample and the products are summed. If the signal contains an oscillation at this frequency, the products add up to a large value; otherwise they cancel out.

The DFT computed “by hand” with dot products equals the result of fft():

fs = 100
n = 200
t = (0:(n - 1)) ./ fs
s = 2 .* cos.(2π * 5 .* t) .+ sin.(2π * 12 .* t)

X_dft = zeros(ComplexF64, n)
for k in 0:(n - 1)
    csw = exp.(-1im * 2π * k .* (0:(n - 1)) ./ n)   # complex sine wave of frequency k × fs / n
    X_dft[k + 1] = sum(s .* csw)                   # dot product
end
maximum(abs.(X_dft .- fft(s)))
5.3799825102094e-12

Lossless Property

The Fourier transform is lossless: the signal can be reconstructed exactly from its Fourier coefficients (by the inverse Fourier transform):

  • A signal of \(N\) samples has \(N\) Fourier coefficients. For a real signal, the coefficients of negative and positive frequencies are complex conjugates, so \(N\) real numbers in the time domain correspond to \(N\) independent real numbers (real and imaginary parts) in the frequency domain.
  • The transformation is one-to-one: every signal has exactly one set of Fourier coefficients and vice versa.
  • The complex sine waves of the DFT are orthogonal, so the signal energy is preserved (Parseval’s theorem) and the coefficients explain 100% of the signal variance.

Complex Fourier Coefficients

The Fourier transform returns a complex number \(X[k] = a + bi\) (a Fourier coefficient) for each frequency, where \(a\) is the real part (the similarity with a cosine) and \(b\) is the imaginary part (the similarity with a sine).

x = rand(5)
X = fft(x)
5-element Vector{ComplexF64}:
   2.129112604501158 + 0.0im
 -0.2505075822321687 - 0.5962743907819374im
 -0.1552506856378541 - 0.9041568752648635im
 -0.1552506856378541 + 0.9041568752648635im
 -0.2505075822321687 + 0.5962743907819374im

Magnitude (absolute value): the strength of the frequency component:

\[ \vert X \vert = \sqrt{a^2 + b^2} \]

mag = abs.(X)
5-element Vector{Float64}:
 2.129112604501158
 0.6467589951892261
 0.9173889199678323
 0.9173889199678323
 0.6467589951892261

Phase (angle, in radians): the timing of the frequency component, i.e. where its peaks are relative to the start of the signal:

\[ \theta = \operatorname{atan2}(b, a) \]

ph = angle.(X)
5-element Vector{Float64}:
  0.0
 -1.9685274393216505
 -1.740845735212179
  1.740845735212179
  1.9685274393216505

angle() (or atan(b, a)) must be used instead of \(\arctan(b / a)\): the two-argument arctangent returns angles in the correct quadrant (\(-\pi\) to \(\pi\)), and does not fail for \(a = 0\).

Phase is independent of amplitude, but when the amplitude is 0 the phase is undefined, and for small amplitudes it is dominated by noise. Coefficients of (near) zero magnitude should therefore be ignored (or their phase set to 0) when phases are interpreted.

Phases are used in neuroscience:

  • in time-frequency analysis: the consistency of the phase over repeated trials (inter-trial phase clustering, ITPC)
  • in connectivity analysis: the consistency of phase differences between electrodes (ISPC, PLV, PLI) as a measure of functional coupling
  • in spike-field coherence: whether action potentials are locked to a phase of the local field potential

Phase angles are circular and must not be averaged directly (the mean of 350° and 10° is 0°, not 180°). Instead, the angles are represented as unit vectors \(e^{i\theta}\) and the vectors are averaged: the angle of the mean vector is the mean phase and its length (0–1) shows how consistent the phases are.

θ = deg2rad.([350, 10, 5, 355])
m = mean(exp.(1im .* θ))
(arithmetic_mean = rad2deg(mean(θ)), circular_mean = rad2deg(angle(m)), vector_length = abs(m))
(arithmetic_mean = 180.0, circular_mean = -3.0103647900187057e-15, vector_length = 0.9905012255519767)

Frequencies and Frequency Resolution

The DFT of \(N\) samples recorded at the sampling rate \(f_s\) gives coefficients at the frequencies

\[ f_k = k \frac{f_s}{N}, \qquad k = 0, 1, \dots \]

starting at 0 Hz (the mean of the signal), regardless of the sampling rate and the number of samples. The distance between the frequencies, the frequency resolution (bin width), is:

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

where \(T\) is the duration of the signal. It depends only on the duration: 1 second of data gives 1 Hz bins, 5 seconds give 0.2 Hz bins. To distinguish two oscillations whose frequencies differ by \(\Delta f\), the signal must be at least \(1 / \Delta f\) seconds long.

freqs() returns the frequencies of the one-sided spectrum (0 Hz to the Nyquist frequency) and the Nyquist frequency:

fs = 10
n = 10
f, nyquist = freqs(n, fs)
([0.0, 1.0, 2.0, 3.0, 4.0, 5.0], 5.0)

Positive and Negative Frequencies

The DFT returns \(N\) coefficients: for \(k\) up to \(N / 2\) they correspond to frequencies from 0 to the Nyquist frequency (\(f_s / 2\)); the remaining coefficients (above the Nyquist frequency) correspond to negative frequencies, in reverse order.

Negative frequencies are a consequence of representing signals with complex exponentials: a real cosine is the sum of two complex sine waves rotating in opposite directions, \(\cos(2\pi f t) = \frac{1}{2}(e^{i 2\pi f t} + e^{-i 2\pi f t})\). For a real signal (e.g. EEG), the coefficients of negative frequencies are the complex conjugates of the positive ones: the spectrum is symmetric and the amplitude of each component is split equally between the positive and the negative frequency.

Two-sided amplitude spectrum of a signal with 50 Hz (amplitude 0.7) and 120 Hz (amplitude 1) components and noise, sampled at 1000 Hz:

fs = 1000
t = (0:(1500 - 1)) ./ fs
s = 0.7 .* sin.(2π * 50 .* t) .+ sin.(2π * 120 .* t) .+ 2 .* randn(length(t))
X = fft(s)
f2 = fftshift(fftfreq(length(s), fs))       # -500 … 500 Hz
lines(f2, fftshift(abs.(X)) ./ length(s);
      axis = (xlabel = "Frequency [Hz]", ylabel = "|X| / N"))
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/dVJIx/src/scenes.jl:264

Therefore, for real signals only the positive frequencies are used (rfft() computes only these), and the amplitudes are doubled to account for the negative frequencies, except for 0 Hz and the Nyquist frequency, which have no negative counterpart.

The negative frequencies are needed for the inverse Fourier transform with ifft() (irfft() reconstructs a real signal from the positive frequencies only) and they differ from the positive ones for complex signals.

Amplitude Spectrum

The amplitude spectrum is the magnitude of the Fourier coefficients, normalized to the units of the signal: divided by \(N\) and doubled for the positive frequencies (one-sided spectrum). A sine of amplitude \(A\) then has the amplitude \(A\) at its frequency.

fs = 100
t = (0:(200 - 1)) ./ fs
s = 2 .* cos.(2π * 5 .* t) .+ sin.(2π * 12 .* t)
n = length(s)
X = rfft(s)
A = abs.(X) ./ n
A[2:(end - 1)] .*= 2              # all except 0 Hz and Nyquist (even n)
f, _ = freqs(n, fs)
plot_psd(f, A; xlabel = "Frequency [Hz]", ylabel = "Amplitude [AU]")
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/dVJIx/src/scenes.jl:264

The amplitudes at 5 and 12 Hz are 2 and 1, as in the signal. The units of the amplitude spectrum are the units of the signal (e.g. μV).

Power Spectrum

The power spectrum is the squared amplitude spectrum: how the power (mean square) of a signal is distributed over the frequencies:

\[ P[k] = \frac{\vert X[k] \vert^2}{N^2} \]

(doubled for the positive frequencies of the one-sided spectrum, except 0 Hz and the Nyquist frequency). With this normalization a sine of amplitude \(A\) has the power \(A^2 / 2\) (its mean square), and the powers of all frequencies sum to the mean square of the signal (Parseval’s theorem). The units are the squared units of the signal (e.g. μV²).

P = abs2.(X) ./ n^2
P[2:(end - 1)] .*= 2
(power_5Hz = P[f .== 5][1], power_12Hz = P[f .== 12][1],
 sum_of_powers = sum(P), mean_square = mean(s .^ 2))
(power_5Hz = 2.0, power_12Hz = 0.5000000000000002, sum_of_powers = 2.5, mean_square = 2.4999999999999996)

ftransform() returns the Fourier coefficients, amplitudes, powers and phases with this normalization (one-sided by default).

Amplitude and power spectra are appropriate for signals composed of distinct oscillations (a discrete spectrum); their values do not depend on the signal duration. For signals with a continuous spectrum (noise, EEG), the power is spread over the frequencies and the power per bin decreases with the bin width: the power spectral density should be used.

Power Spectral Density

The power spectral density (PSD) is the power per unit of frequency (e.g. μV²/Hz): the power spectrum divided by the bin width \(f_s / N\):

\[ PSD[k] = \frac{\vert X[k] \vert^2}{f_s N} \]

(doubled for the one-sided PSD, except 0 Hz and the Nyquist frequency). The integral of the PSD over frequency (the sum of the PSD values × the bin width) equals the mean square of the signal. The PSD of a noise does not depend on the signal duration, which makes PSDs of different recordings comparable.

The PSD computed in this way (from the DFT of the whole signal) is the periodogram:

fs = 250
N = 1000                                   # 4 seconds
t = (0:(N - 1)) ./ fs
x = generate_sine(10, t, 1) .+ generate_sine(20, t, 0.5) .+ 0.5 .* randn(N)

X = rfft(x)
f, _ = freqs(N, fs)
PSD = abs2.(X) ./ (fs * N)
PSD[2:(end - 1)] .*= 2

plot_psd(f, PSD;
         xlabel = "Frequency [Hz]",
         ylabel = "PSD [AU²/Hz]",
         title = "Periodogram")
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/dVJIx/src/scenes.jl:264

(integral_of_psd = sum(PSD) * fs / N, mean_square = mean(x .^ 2))
(integral_of_psd = 0.8510371551626658, mean_square = 0.8510371551626658)

The periodogram of a noisy signal is noisy itself: its variance does not decrease with a longer signal. Better estimates:

  • Welch’s method: the signal is divided into overlapping segments, each segment is windowed and its periodogram is calculated; the periodograms are averaged (less variance, lower frequency resolution).
  • Multitaper method: the periodograms of the signal multiplied by several orthogonal tapers (Slepian sequences) are averaged.
p_welch = welch_pgram(x, 250; fs = fs)          # 1-second segments, 50% overlap
plot_psd(Float64.(freq(p_welch)), power(p_welch);
         xlabel = "Frequency [Hz]",
         ylabel = "PSD [AU²/Hz]",
         title = "Welch's method",
         flim = (0, 40))
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/dVJIx/src/scenes.jl:264

In NeuroAnalyzer, psd() calculates the PSD of a NEURO object using Welch’s method (default), the multitaper method and others.

Spectral Leakage and Windowing

The DFT treats the signal as one period of an infinitely repeated signal. If an oscillation does not complete an integer number of cycles within the signal, the repeated signal has a discontinuity at its edges, and the power of the oscillation leaks into neighbouring frequencies: the peak is lower and wider and the spectrum shows side lobes.

fs = 100
t = (0:(100 - 1)) ./ fs                    # 1 second: 1 Hz bins
s_int = sin.(2π * 10 .* t)                 # 10 cycles
s_frac = sin.(2π * 10.5 .* t)              # 10.5 cycles
f, _ = freqs(length(t), fs)
fig, ax = stem(f, ftransform(s_int).a; label = "10 Hz (integer number of cycles)",
               axis = (xlabel = "Frequency [Hz]", ylabel = "Amplitude"))
stem!(ax, f, ftransform(s_frac).a; color = :red, label = "10.5 Hz (leakage)")
xlims!(ax, 0, 25)
axislegend(ax; position = :rt)
fig
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/dVJIx/src/scenes.jl:264

Leakage is reduced by windowing (tapering): the signal is multiplied by a window that falls smoothly to zero at the edges (e.g. Hann), which removes the discontinuity. The window reduces the side lobes, but widens the main peak and reduces the amplitude (for a Hann window by half, the coherent gain of 0.5, which can be corrected by dividing by the mean of the window):

w = generate_window(:hann, length(t))
a_rect = ftransform(s_frac).a
a_hann = ftransform(s_frac .* w).a ./ mean(w)
fig, ax = lines(f, 20 .* log10.(a_rect); label = "rectangular (no window)",
                axis = (xlabel = "Frequency [Hz]", ylabel = "Amplitude [dB]"))
lines!(ax, f, 20 .* log10.(a_hann); color = :red, label = "Hann window")
axislegend(ax; position = :rt)
fig
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/dVJIx/src/scenes.jl:264

Zero Padding

Zero padding appends zeros to the signal before the FFT. A signal of \(L\) samples padded with \(n\) zeros has \(L + n\) Fourier coefficients, so the frequency bins are denser: \(\Delta f = f_s / (L + n)\).

Zero padding does not add information and does not improve the frequency resolution (the ability to separate close frequencies, which depends only on the signal duration): it interpolates the spectrum, which makes it smoother and helps to locate peaks between the original bins.

fs = 100
s = sin.(2π * 10.5 .* (0:99) ./ fs)        # 1 second, 10.5 Hz
f1, _ = freqs(100, fs)
f4, _ = freqs(400, fs)
fig, ax = lines(f4, ftransform(s; pad = 300).a; color = :red, label = "padded to 400 samples",
                axis = (xlabel = "Frequency [Hz]", ylabel = "Amplitude"))
stem!(ax, f1, ftransform(s).a; label = "100 samples")
xlims!(ax, 0, 25)
axislegend(ax; position = :rt)
fig
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.

└ @ Makie ~/.julia/packages/Makie/dVJIx/src/scenes.jl:264

To keep the amplitudes in the units of the signal, the coefficients are normalized by the original signal length \(L\), not by \(L + n\) (ftransform(s; pad = n) does this).

Zero padding is also used:

  • to make the FFT faster: FFT algorithms are fastest for lengths that are powers of two (nextpow2() returns the next power of two, fft2() pads and transforms)
  • for convolution in the frequency domain: to compute the linear convolution of signals of lengths \(N_1\) and \(N_2\), both are padded to \(N_1 + N_2 - 1\) samples

NeuroAnalyzer functions: pad0(s, n) appends n zeros, fft0(s, n) and rfft0(s, n) pad and transform (the coefficients are not normalized), ifft0(X, n) inverts and removes the n padded samples.

s = rand(10)
l = nextpow2(length(s))
X = fft0(s, l - length(s))
length(X)
16

Tip: Always report padding (e.g. “zero-padded from 1000 to 2048 samples”).

Euler’s Formula

Euler’s formula connects exponential, trigonometric and complex number representations:

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

For a real angle \(\theta\), \(e^{i \theta}\) is a unit vector (magnitude 1) in the complex plane at the angle \(\theta\) from the positive real axis:

  • real part: \(\Re = \cos(\theta)\) (x-coordinate)
  • imaginary part: \(\Im = \sin(\theta)\) (y-coordinate)

A vector of magnitude \(m\) and angle \(\theta\):

\[ m e^{i \theta} = m \cos(\theta) + i \, m \sin(\theta) \]

Conversely, any complex number \(z = x + iy\) can be written in polar form:

\[ z = \vert z \vert e^{i\theta}, \qquad \vert z \vert = \sqrt{x^2 + y^2}, \qquad \theta = \operatorname{atan2}(y, x) \]

For \(\theta = \pi\) Euler’s formula gives Euler’s identity \(e^{i \pi} + 1 = 0\).

Euler’s formula is the basis of the Fourier transform: a complex sine wave \(e^{i 2\pi f t}\) is a vector rotating \(f\) times per second, and a Fourier coefficient is a vector whose magnitude is the amplitude and whose angle is the phase of the frequency component. The same representation is used by the Hilbert transform and complex wavelets to extract the instantaneous amplitude (or power, \(\vert z \vert^2\)) and phase.

θ = π / 4                                  # 45°
m = 1.2
z = m * exp(1im * θ)                       # or: m * (cos(θ) + 1im * sin(θ))
plot_polar([angle(z) abs(z)]; ticks = true)

Example

A 2-second cosine of 10 Hz frequency, amplitude 2 (arbitrary units) and phase +30° (0.52 rad), sampled at 1000 Hz:

fs = 1000
t = collect(0:(1 / fs):(2 - 1 / fs))      # 2000 samples: exactly 20 cycles
s = generate_cosine(10, t, 2, 30)
# X contains: Fourier coefficients (X.c), amplitudes (X.a), powers (X.p), phases (X.ph)
X = ftransform(s)
f, _ = freqs(length(s), fs)
plot_signal(t, s; ylabel = "Amplitude [AU]")

Amplitude spectrum (amplitude 2 at 10 Hz):

plot_psd(f, X.a; flim = (0, 40), xlabel = "Frequency [Hz]", ylabel = "Amplitude [AU]")

Power spectrum (power \(2^2 / 2 = 2\) at 10 Hz):

plot_psd(f, X.p; flim = (0, 40), xlabel = "Frequency [Hz]", ylabel = "Power [AU²]")

The phases of frequencies without signal are given by numerical noise; they are set to zero:

X.ph[X.a .< 0.01 * maximum(X.a)] .= 0
plot_phase(rad2deg.(X.ph[1:50]), f[1:50]; unit = :deg, type = :stem)

The phase at 10 Hz is 30°, as expected: the DFT measures phase relative to a cosine, so a sine has the phase −90°. This holds when the signal contains an integer number of cycles; otherwise leakage shifts the phase as well.

Inverse Fourier Transform

Because the Fourier transform is lossless, the inverse Fourier transform reconstructs the original signal from its Fourier coefficients:

\[ s[n] = \frac{1}{N} \sum_{k = 0}^{N - 1} X[k] e^{i 2\pi kn/N} \]

The inverse transform is a sum of complex sine waves (one per frequency) weighted by the Fourier coefficients.

X = ftransform(s; nf = true)              # all (positive and negative) frequencies
s_ifft = real(ifft(X.c))
maximum(abs.(s_ifft .- s))                 # reconstruction error
1.3322676295501878e-15

Tip: ftransform() returns coefficients for positive frequencies only by default (it uses rfft()). For the inverse transform with ifft(), the negative frequencies are needed: use ftransform(s; nf = true). Alternatively, irfft(X.c, length(s)) reconstructs a real signal from the positive frequencies.

X = ftransform(s)
maximum(abs.(irfft(X.c, length(s)) .- s))
1.3322676295501878e-15

Modifying the Fourier coefficients before the inverse transform is a simple way of filtering: e.g. setting the coefficients of unwanted frequencies to zero removes them from the signal (in practice, filters with smooth frequency responses are preferred, see Filtering).

Non-Stationary Signals

The Fourier transform describes which frequencies are present in the whole signal, but not when: the coefficients are computed from all samples. For a stationary signal (whose frequency content does not change over time) this is sufficient; for non-stationary signals the spectrum can be misleading:

  • smeared peaks: a frequency that changes over time (e.g. a chirp) is spread over a range of frequencies
  • lower amplitudes: the power of an oscillation present only for a short time is averaged over the whole signal
  • no timing: an oscillation at the beginning and one at the end of the signal give the same spectrum

A chirp (2 → 20 Hz) and two consecutive oscillations (5 Hz, then 15 Hz) have broad spectra that do not show how their frequency changes:

fs = 200
t = collect(0:(1 / fs):(4 - 1 / fs))
chirp = @. sin(2π * (2 * t + (20 - 2) / (2 * 4) * t^2))
steps = [ti < 2 ? sin(2π * 5 * ti) : sin(2π * 15 * ti) for ti in t]
f, _ = freqs(length(t), fs)
fig = Figure(size = (900, 450))
lines(fig[1, 1], t, chirp; axis = (title = "chirp 2–20 Hz", xlabel = "Time [s]"))
lines(fig[1, 2], f, ftransform(chirp).a; axis = (title = "amplitude spectrum", xlabel = "Frequency [Hz]"))
lines(fig[2, 1], t, steps; axis = (title = "5 Hz, then 15 Hz", xlabel = "Time [s]"))
lines(fig[2, 2], f, ftransform(steps).a; axis = (title = "amplitude spectrum", xlabel = "Frequency [Hz]"))
fig

EEG is non-stationary: its frequency content changes with the state of the subject and in response to events. It is therefore analyzed in short segments (epochs, which are assumed to be approximately stationary) or with time-frequency methods (short-time Fourier transform, wavelets, Hilbert transform of band-pass filtered signals), which estimate the spectrum as a function of time; see Hilbert Transform, Wavelet Transform and Spectrogram.

Summary

Quantity Formula (one-sided, \(0 < f_k < f_s / 2\)) Units ftransform()
Fourier coefficient \(X[k]\) c
Amplitude \(2 \vert X[k] \vert / N\) signal units (e.g. μV) a
Power \(2 \vert X[k] \vert^2 / N^2\) units² (e.g. μV²) p
Phase \(\operatorname{atan2}(\Im X[k], \Re X[k])\) rad ph
PSD \(2 \vert X[k] \vert^2 / (f_s N)\) units²/Hz (e.g. μV²/Hz) (psd())
Frequency \(k f_s / N\) Hz (freqs())
Frequency resolution \(f_s / N = 1 / T\) Hz