NeuroAnalyzer tutorials: Introduction to Fourier transform

Initialize NeuroAnalyzer 
using Plots
using FFTW
using DSP
using NeuroAnalyzer;

Fourier transform

Fourier transform is a mathematical model that decomposes a signal into its constituent frequencies. It helps to transform the signal between the frequency domain to the time domain.

Fourier theorem: any signal can be represented as a combination of sine waves, each with its own frequency, amplitude and phase.

Fourier transform works by comparing the signal \(s(t)\) to the sinusoids oscillating at frequency \(f_j\); when the data and sinusoid at frequency \(f_j\) align (their dot product is not 0), the summation in the Fourier transform is large and the resulting \(X_j\) component is a large number.

Fourier transform and many other frequency-analysis methods assume that signal is stationary. For example, for a chirp (which is a non-stationary signal): FFT(chirp) = FFT(reverse(chirp)). This shows how FFT is unable to locate changes of frequency over time.

               ^      ___
    ___^___ __/ \__ _/   \_
       1       2       3
  1. sharp peaks: frequency stationary, no changes of the frequency over time
  2. broad peaks: non-stationary amplitude
  3. flat plateau: non-stationary frequency

Edge artifacts: on sudden changes of amplitude.

FFT of non-stationary signal:

  • smeared out peaks
  • side-band ripples around peaks
  • amplitudes smaller than in signal

For repeated measurements we may use:

mean(fft(s))

Discrete Fourier transform (DFT): procedure of comparing similarities (dot products) of the signal with a series of different complex sine waves of the same length as the signal - therefore temporal data is lost. We do not use DFT, we use FFT (fast Fourier transform), FFT is much faster (FFT: \(O(N \times \log_2(N))\), DFT: \(O(N^2)\)).

\[ X[k] = \sum_n={0}^{N-1} x_n \times e^{-\frac{i 2 \pi} {N}k_n} = \sum_n={0}^{N-1} x_n \times e^{-\frac{i 2 \pi} {N}k_n} \]

We are multiplying a discrete signal (\(x_n\)) by a series of sines and cosines at different frequencies (\(k_n\)).

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

Fourier transformation is lossless:

  • n time points = n sine waves
  • degrees of freedom = 0
  • 100% variance explained
  • each Fourier coefficient = one frequency

Hence, inverse Fourier transform (IFT) is possible.

Fourier transform returns a vector of complex numbers (Fourier coefficients) that represent the phase and the magnitude of frequencies of the signal.

x = rand(10)
mag = abs.(fft(x))
10-element Vector{Float64}:
 4.164129696124421
 1.651060155948692
 1.2332503183107248
 0.5668076594884424
 0.5147418573777112
 0.09511951328427826
 0.5147418573777112
 0.5668076594884424
 1.2332503183107248
 1.651060155948692

But this does not relate to physical properties of the signal, just shows the composition of frequencies in the signal.

Spectrum: some parameter (e.g. amplitude or power) as a function of frequency.

x = rand(10)
fs = 5
n = length(x)
10

Amplitude spectrum:

amp = abs.(fft(x)) ./ n
10-element Vector{Float64}:
 0.5361584241922743
 0.16980307747942175
 0.12003395109219954
 0.09209127510373195
 0.02236074568472433
 0.04216559468411152
 0.02236074568472433
 0.09209127510373195
 0.12003395109219954
 0.16980307747942175

To get the Fourier coefficients in the same scale as the original data, we must divide the coefficients by the number of samples.

To reconstruct amplitudes from negative frequencies we multiple amplitudes of positive frequencies by 2. However: amplitude at 0 Hz and Nyquist frequency also gets doubled this way; correctly all frequencies above 0 and below Nyquist should be doubled.

amp[2:end - 1] .*= 2
8-element view(::Vector{Float64}, 2:9) with eltype Float64:
 0.3396061549588435
 0.24006790218439908
 0.1841825502074639
 0.04472149136944866
 0.08433118936822304
 0.04472149136944866
 0.1841825502074639
 0.24006790218439908

Power spectrum:

pow = amp.^2
10-element Vector{Float64}:
 0.2874658558323428
 0.11533234048593002
 0.0576325976592182
 0.03392321180092496
 0.002000011790307671
 0.007111749500259094
 0.002000011790307671
 0.03392321180092496
 0.0576325976592182
 0.028833085121482506

Power spectrum shows power (strength) of the signal at each of the component frequencies of the signal. Therefore, it converts the signal from time domain to frequency domain.

\[ PS = \frac{2\vert X \vert^2}{N^2} \] \(X\): Fourier-transformed (two-sided) signal

These methods are valid if we are using a rectangular window, the signal is stationary and we calculate discrete power spectrum (which means that individual frequencies are distanced; otherwise spectra of these frequencies will merge and nothing can be said about individual frequencies).

Power spectral density (PSD): power spectrum normalized for the sampling frequency.

\[ PSD = \frac{2\vert X \vert^2}{f_s \times N^2} \] \(X\): Fourier-transformed (two-sided) signal

psd = 1/(fs * n) * abs.(fft(x)).^2
10-element Vector{Float64}:
 0.5749317116646855
 0.05766617024296502
 0.028816298829609104
 0.01696160590046248
 0.0010000058951538356
 0.0035558747501295467
 0.0010000058951538356
 0.01696160590046248
 0.028816298829609104
 0.05766617024296502

Amplitude/power spectrum are useful for signal of discrete spectrum, for a signal with continuous spectrum, PSD should be used.

Bin width: distance between to samples in the spectrum.

\[ bw = \frac{f_s}{n} \] \(f_s\): sampling rate

\(n\): number of samples

Therefore, to have a more narrow distance between separate frequencies in the spectrum, we need longer data (more samples) = more frequency resolution.

Same may be obtain by padding the signal with 0s. However, zero-padding only increases visualization, but not frequency resolution.

Frequencies:

k = 0:n/2     # 0.0:1.0:5.0
bw = fs / n   # 0.5
freq = collect(k .* bw)
6-element Vector{Float64}:
 0.0
 0.5
 1.0
 1.5
 2.0
 2.5

Units (for EEG):

  • amplitude \(\mathrm{μV}\)
  • power: \(\mathrm{μV^2}\)
  • PSD: \(\mathrm{μV^2/Hz}\)

Positive and negative frequencies

Fourier transform returns coefficients above Nyquist frequency; these are negative frequencies - below frequencies going backward.

Negative frequencies: sine waves with opposite directions (positive frequencies: counterclockwise angular motion, negative frequencies: clockwise angular motion), they go faster and faster, but due to limits set by Nyquist frequency we cannot measure them precisely (aliasing).

For real-valued signal (e.g. EEG) the negative frequencies mirror the positive frequencies and the amplitude gets split between the positive and negative frequencies. The negative frequencies are important only for signal of complex numbers (they are different than positive signal then), in neuroscience we usually ignore negative frequencies and double positive frequencies instead (we use rfft() instead of fft(), which requires FFTW.jl). However, negative frequencies are necessary for inverse Fourier transform.

Euler’s formula

Euler’s formula for a unit vector: \[ e^{1i \times \theta} = \cos(\theta) + 1i \times \sin(\theta) \] \(\theta\): phase

These can be used interchangeably: \[ \Re = \cos(\theta) \]

\[ \Im = \sin(\theta) \]

Raising a complex number to the natural log base gives a vector with absolute value of 1 and an angle of the real component of the number.

\(e^{1i \times \theta}\) (exp(1im * theta)) is an efficient way to represent oscillatory information using polar notation.

Euler’s formula for \(\theta = \pi\): \[ e^{\pi i} + 1 = 0 \]

Euler’s formula for vectors larger than \([1, 1]\): \[ m \times e^{1i \theta} = m \times (\cos(\theta) + 1i \times \sin(\theta)) \] \(m\): magnitude of the vector

Used for extracting power and phase information; phase is a temporal information (when a frequency occurred).

\[ \Re = z \times \cos(\theta) \] \[ \Im = z \times \sin(\theta) \] \(z\): line length (amplitude)

\(\theta\): angle (phase)

theta = pi / 4
z = 1.5
plot_polar([theta z], ticks=true)
euler = z * exp(1im * theta)
# or
euler = z * (cos(theta) + 1im * sin(theta))
plot_polar([DSP.angle(euler) abs(euler)], ticks=true)

Fourier coefficients

Fourier coefficients: a set of n complex sine waves.

For each coefficient we check how much it is in the signal plane (dot product): \[ a = \vert s \vert \vert C \vert \cos(\theta) \] \(s\): signal

\(C\): Fourier coefficient

\(\theta\): angle between these

For each Fourier coefficient:

\[ C = m \times e^{ik} \] \(m\): amplitude

\(k\): phase

Therefore: \[ \vert s \vert \vert C \vert \cos(\theta) = m \times e^{ik} \]

Phase is independent of amplitude, BUT: when amplitude is 0, phase is undefined. Therefore, for small amplitude, phase estimation may be less accurate.

To calculate amplitude (magnitude): abs.(C)

\[ \vert 3+2i \vert = \sqrt{3^2 + 2^2} \]

To calculate phase:

\[ \theta = \arctan(\Im(X) / \Re(X)) \] \(\theta\): phase angle in radians

\(X\): Fourier coefficient

\(\Im(X)\): imaginary part

\(\Re(X)\): real part

Due to the division, it is practically useful to remove very small values of \(\vert X \vert\) since they will affect calculations. We also use the angle() function (theta = DSP.angle.(X)) instead of theta = atan.(imag.(X), real.(X)).

Phase angles have several uses in neuroscience:

  • in time-frequency analyses: to quantify the consistency of the timing of rhythmic activity over repeated trials
  • in synchronization analyses: are compared across different electrodes to determine whether two brain regions are coupled (functionally interacting)
  • in spike-field coherence: used to determine whether the timing of action potentials is locked to the population-level field potentials

Phase angles MUST NOT be averaged; instead, calculate an average vector, its length represents how close the angles are (e.g. /_ or _._).

Example

Generate 2-second cosine wave of 10 Hz frequency, amplitude of 2 (arbitrary units) and phase shift of \(+30 \degree\) (\(\approx 0.52\) radians) and perform the Fourier transform to calculate amplitude and phase:

fs = 1000
t = collect(0:1/fs:2)
s = generate_cosine(10, t, 2, 30)
# X will contain Fourier coefficients (X.c), amplitudes (X.a), powers (X.p) and phases (X.ph)
X = ftransform(s)
# f will contain vector of frequencies
f, _ = freqs(s, fs)
([0.0, 0.5, 1.0, 1.499, 1.999, 2.499, 2.999, 3.498, 3.998, 4.498  …  495.252, 495.752, 496.252, 496.752, 497.251, 497.751, 498.251, 498.751, 499.25, 499.75], 500.0)

Plot the signal:

plot_signal(t,
            s,
            ylabel="Amplitude [AU]")

Plot amplitude spectrum:

plot_psd(f,
         X.a,
         frq_lim=(0, 50),
         ylabel="Amplitude [AU]")

Plot power spectrum:

plot_psd(f,
         X.p,
         frq_lim=(0, 50),
         ylabel="Power [AU^2]")

Same results can be obtained using a periodogram:

S = DSP.periodogram(s, fs=1000, nfft=length(s), window=rect(length(s)))
plot_psd(Float64.(S.freq),
         S.power,
         frq_lim=(0, 50),
         ylabel="Power [AU^2]")

The frequencies where there is no signal (which is all except for the 10 Hz frequency), the phase will be given by noise. You can set the phase here to zero:

X.ph[X.a .< maximum(X.a)] .= 0;

Plot phases:

plot_phase(rad2deg.(X.ph[1:100]),
           f[1:100],
           unit=:deg,
           type=:stem)

There is a \(30 \degree\) phase shift at 10 Hz, as expected.

Remember that sine has initial phase of \(- \pi /2\) (\(-90 \degree\)). Also, if the input signal does not have an integer of periods, it will cause some phase shift.

To avoid edge artifacts, windowing (e.g. Hann) could be applied:

w = generate_window(:hann, length(s))
sw = s .* w
plot_signal(t, sw)

Re-run the Fourier transform:

Xw = ftransform(sw)
(c = ComplexF64[-6.378231276471524e-14 + 0.0im, -1.632669818663241e-5 - 6.359892484473871e-5im, -6.696689848452953e-5 - 0.0001331899928802468im, -0.0001571723923801749 - 0.00021555111107163786im, -0.0002966704345595964 - 0.00031919647452030374im, -0.0005013939636830131 - 0.0004557090906599369im, -0.000796549617555449 - 0.000641801517967411im, -0.00122198248958868 - 0.0009027203281941399im, -0.001841722868043523 - 0.0012781668535491503im, -0.0027615544080235987 - 0.00183309825819405im  …  -3.83574101663198e-8 + 3.2578698232879354e-8im, -3.84477185504071e-8 + 2.914899142756848e-8im, -3.852766157638983e-8 + 2.5719239482807973e-8im, -3.859763391076215e-8 + 2.2289776118313937e-8im, -3.865781277830343e-8 + 1.8860272669211753e-8im, -3.870746794087089e-8 + 1.543082415804777e-8im, -3.874763613613868e-8 + 1.2002194984234675e-8im, -3.877796240263695e-8 + 8.57254886917625e-9im, -3.87976162388169e-8 + 5.143543647008492e-9im, -3.8807702215364205e-8 + 1.7145071732005799e-9im], a = [3.1875218772971134e-17, 6.562831829441401e-8, 1.490031281999846e-7, 2.666352008022118e-7, 4.3555709003733853e-7, 6.772059948910398e-7, 1.0224259706642823e-6, 1.5184988576077471e-6, 2.2406766257733337e-6, 3.3129222410177625e-6  …  5.0300414799029475e-11, 4.822410825476877e-11, 4.6300277909052173e-11, 4.454914522655595e-11, 4.299171576319334e-11, 4.164904958816837e-11, 4.0543652916817954e-11, 3.9694368937906417e-11, 3.911752170081906e-11, 3.882614376649069e-11], p = [1.0160295718247715e-33, 4.307076162152917e-15, 2.2201932213381045e-14, 7.10943303068358e-14, 1.8970997868179423e-13, 4.58607959516363e-13, 1.0453548654887998e-12, 2.305838780556033e-12, 5.020631741286972e-12, 1.0975453775030154e-11  …  2.5301317289544236e-21, 2.325564616967657e-21, 2.143715734455465e-21, 1.9846263404167727e-21, 1.8482876242632064e-21, 1.734643331597708e-21, 1.643787791839401e-21, 1.5756429253786297e-21, 1.5301805040140501e-21, 1.5074694397762037e-21], ph = [3.141592653589793, -1.8220834640678587, -2.0366752536497286, -2.200829349559275, -2.319634638604541, -2.4038907379823655, -2.4633693927330635, -2.505340693315323, -2.5349007975847275, -2.5555827108601727  …  2.437478614660359, 2.4928983348405285, 2.552976983430602, 2.617888532283946, 2.687690082457542, 2.7622485587869816, 2.8412123802331344, 2.924024317811847, 3.0097875572525767, 3.0974418045977408])

Plot amplitudes:

plot_psd(f,
         Xw.a,
         frq_lim=(0, 50),
         ylabel="Amplitude [AU]")

Plot phases:

Xw.ph[Xw.a .< maximum(Xw.a)] .= 0
plot_phase(rad2deg.(Xw.ph[1:100]),
           f[1:100],
           unit=:deg,
           type=:stem)

Inverse Fourier transform:

X = ftransform(s, nf=true)
s_ifft = real(ifft0(X.c)) .* length(s)
plot_signal(t,
            s)
Plots.plot!(t,
            s_ifft,
            ls=:dot,
            lc=:red,
            label=false)

(!) Multiply by the length of the signal to re-normalize the reconstructed signal.

(!) By default NeuroAnalyzer transform() and ftransform() functions use Fourier transform and returns coefficients for positive frequencies only (internally it uses rfft()). For inverse Fourier transform we need positive and negative frequencies in the Fourier components. To do so, set nf=true for transform() or ftransform().

Frequency resolution

Frequency step size: determines minimal difference between two frequencies that can be separated using Fourier transform.

The frequencies in the Fourier transform are in normalized units, where 1 is the data sampling rate.

The frequencies if the Fourier transform will always start at 0 Hz, regardless of the data sampling rate or the number of time points.

Frequency resolution:

\[ df = 1/T \] \(T\): length of the signal (in seconds)

If we want 0.2 Hz frequency resolution, the signal must be

\[ T = 1/0.2 = 5.0 \]

seconds long.

With longer signal there are more frequencies of which we estimate powers, but individual estimates are not becoming more accurate.

Zero padding

Zero padding: adding zeros at the end of the signal before FFT - adding zeros after a time series produces a smoother (interpolated) Fourier representation:

FFT(signal, LENGTH(signal) + number of zeros)

When normalizing, we use the length of the original signal, not the padded one:

FFT(signal, LENGTH(signal) + number of zeros) / LENGTH(signal)

In NeuroAnalyzer there is fft0(s, n) function available, which will pad the signal s with n zeros.

Signal length:

length(s)
2001

Fourier coefficients of the padded signal:

X = fft0(s, 100) ./ length(s)
2101-element Vector{ComplexF64}:
  0.0008655926074806658 + 1.1096681905298916e-19im
  0.0006807916986280217 + 0.0025291539891134306im
  0.0001208586152052872 + 0.005037839531709647im
 -0.0008295773090957877 + 0.007499324595565706im
 -0.0021921936960369623 + 0.009875156250158463im
  -0.003989500574161703 + 0.01211123874388396im
  -0.006238005805728693 + 0.014136028706051686im
  -0.008940756783071168 + 0.015861224904430135im
  -0.012080080113014664 + 0.01718508738535727im
  -0.015611346750998853 + 0.017998262458189727im
  -0.019458511740593994 + 0.018191740799362317im
   -0.02351201545263038 + 0.017666362102158648im
  -0.027629387829917924 + 0.01634312862579165im
                        ⋮
  -0.027629387829917904 - 0.01634312862579163im
  -0.023512015452630375 - 0.017666362102158665im
  -0.019458511740593935 - 0.018191740799362317im
  -0.015611346750998864 - 0.017998262458189696im
  -0.012080080113014677 - 0.01718508738535725im
  -0.008940756783071173 - 0.01586122490443014im
   -0.00623800580572868 - 0.014136028706051726im
  -0.003989500574161696 - 0.01211123874388396im
  -0.002192193696036957 - 0.009875156250158452im
 -0.0008295773090957807 - 0.007499324595565722im
 0.00012085861520524464 - 0.0050378395317096335im
  0.0006807916986280158 - 0.002529153989113422im

Padding gives an estimation how the frequency would be if there were more time-points in the signal.

Zero padding allows one to use a longer FFT, which will produce a longer FFT result vector; a longer FFT result has more frequency bins that are more closely spaced in frequency.

Remember: always report padding if used.

When to use zero padding:

  1. to match FFT length for convolution (both signals must have same length)
  2. to obtaining specific frequencies for analysis: increases frequency resolution not frequency precision (double the signal → double the frequency resolution)
  3. to interpolate or smoother plots
  4. FFT works fastest for signal length corresponding to a power of two (use nextpow2() to get the next power of 2 for a given number)

Fourier series

Fourier series: representation of a function as a series of cosines and sines:

\[ \begin{align*} f(t) = a_0 &+ a_1\cos(t) + a_2\cos(2t) + \ldots \\ &+ b_1\sin(t) + b_2\sin(2t) + \ldots \end{align*} \]

\(a\), \(b\): weights, how much of these frequencies is in the signal

For any \(m \neq 0\):

\[ \int_0^{2\pi} \sin(m \times t)\ dt = 0 \]

\[ \int_0^{2\pi} \cos(m \times t)\ dt = 0 \]

Parseval’s theorem

Parseval’s theorem: average power of the signal equals the power of Fourier transformation (sum of squared magnitudes of each of the complex Fourier series components).

\[ P = \sum \vert C \vert^2 \] \(P\): average power

\(C\): Fourier coefficients

Total energy for all frequencies:

\[ TE = \int \vert s(f) \vert^2\ \mathrm{d}f \] \(s\): signal

\(f\): frequency

Energy spectral density: energy per frequency:

\[ E(a, b) = \int_a^b \vert s(f) \vert^2\ \mathrm{d}f \] \(s\): signal

\(f\): frequency

fs = 100
t = 0:1/fs:10
n = length(t)
y = generate_sine(2, t, 2)
yX = fft(y) ./ length(y);
yX_power = sum(@. abs(yX).^2)
1.9980019980019956
y_power = sum(@. abs(y).^2) / n
1.9980019980019956