Euler in the brain

Many analyses in brain imaging (and physics, engineering, etc) rely on a complex representation of a time signal in the form e^ix and we (neuroscientists) are told this characterizes the data fully by combining amplitude and phase. How? there is a weird constant e to the power of x, worst times i! what on Earth is the power of an imaginary value (ix) anyway?

Frequency decomposition as a sum 

Frequency analysis of a signal allows examining the relative contributions of various frequencies.

The classic signal decomposition is the Fourier series, and I think it is best understood using the Sine-Cosine form:

Let's examine this equation by making up a signal, say s(t) is a 6Hz cosine plus a 40Hz sine that is 1/2 in amplitude and shifted in time (that is a phase difference .. but I'll get back to that later).

Figure 1. Made up signal s (in blue) as a sum of the green and red signals

The equation above for the signal s(t) is here simply 0 + acos(2pi*6Hz) + bsin(2pi*40Hz). We have 0 because there is no signal offset, no sums over many n because we know it is 6Hz and 40Hz, overwise we would iterate over many n (see below), we have to find a and b (we already now a=b/2) and yes there is pi, and we'll get back to that.

First Interlude: vectors, triangles and circles

Cartesian space

Using a Euclidian space (using x and y as the name of the 2D axis), we can draw a line characterized by its beginning (x1,y1) and ending (x2,y2) called a vector. A convenient way to represent a vector is to make it starts at (0,0).

Back to our signal: let's break the amplitude axis into a new orthonormal basis, so a given value of my signal is now assigned equally to x and y. Instead of a 2D time by amplitude graph, we have a 3D time by x by y. OK why not? this also allows to make a dynamic 2D representation of the amplitude as a vector that grows and shrinks over time, reverses to negative, etc .. that's fun.

Figure 2: signal s over time in cartesian space

Triangles

If we now decide to look at a given time point, we have a vector with its x and y components. Since the origin is (0,0), it makes a nice triangle. Remember triangles? this naturally leads to the Pythagorean theorem and trigonometry

Pythagorean theorem: hypothenus^2 = adjacent^2+opposite^2

cos(theta) = adjacent/hypotenuse 

sin(theta) = opposite/hypotenuse 

tang(theta) = opposite/adjacent = sin(theta)/cos(theta) 

Circles

When you hear circle, you have to think pi. 

𝝅 = circumference / diameter 

The area of a circle is also  𝝅 * radius^2. This means that if a circle has a radius of 1, the area = 𝝅. Back to our signal: we can plot s on the unit circle, for each time point, we can describe the signal amplitude at a time t as a*hypotenuse or a*(cos(theta)^2 + sin(theta)^2) with the value a being a scalar. Oh wait a minute, this starts looking a lot like our Fourier decomposition!

Representation of S in the unit circle (S on the left, s1 and s2 on the right)

About Degrees and radian: 1 radian, is the angle theta between 2 vectors such as the length on the unit circle is 1 (or on any circle it is equal to the length of the radius) which implies there are 2𝝅 radians in a circle, 360 degrees is thus 2𝝅, 180 degrees = 𝝅, 90 degrees = 𝝅/2, etc .. Back to our signal: This notation is useful because it forces us to think in terms of a circle. Our signal s can now be written at time t as a*(cos(𝝅/L)^2 + sin(𝝅/L)^2) with L a constant value related to the signal length.

Digression

About pi: Sir Isaac Newton came up with a way to compute 𝝅, based on infinite series. I am just dropping that here, but we’ll come back to that soon.

𝝅 = 4*(1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 + 1/13 - 1/15 + 1/17 … ) 

check that cool video from Veritasium discussing this.

Discrete Fourier Transform

Let say we want to know the strengh (a) of the cosine at the frequency m. We can take advantage of the fact (given here without proof) that the integral of cos^2 over the unit circle -pi to +pi is pi.

 

This means that at a given time point x, we evaluate the above equation over the unit circle - but since now the signal s  is described for each time point 'drawn' on the circle, we actually evaluate over the interval 2L: 

The same applies to sine, allowing to get the coefficients a and b from our Fourier equation. In practice, we do not know which frequencies are of interest. Given the signal length (1 sec) and sampling frequency (1000 Hz), the Nyquist and Shannon sampling theorem tells us the lowest frequency we can use here is 2Hz and the largest frequency is 500Hz  (you need at least 2 cycles to quantify, being other the whole time course,  or between 2 samples). The discrete Fourier transform is thus an iteration for each frequency, computing sums over half the signal length of those cosines and sines. 

% Matlab code for an iterative search - discrete Fourier analysis

% ----------------------------------------------------------------------------

L             = range(time_points)/2; % how many time points we can use

dt            = 1/frequency_sampling; % step size

Nyquist_limit = (1/dt)/2;             % max frequency we can compute

a             = zeros(Nyquist_limit-1,1);

b             = zeros(Nyquist_limit-1,1);

for m = 2:Nyquist_limit

    c      = cos(pi/L * m * time_points) * dt;

    a(m-1) = round(sum(S.* c) / L);

    s      = sin(pi/L * m * time_points) * dt;

    b(m-1) = round(sum(S.* s) / L);

end

The result is what we expect:

  frequency    cosine strength    sine strength
    _________    _______________    _____________

        2               0                 0      
        3               0                 0      
        4               0                 0      
        5               0                 0      
        6               2                 0      
        7               0                 0      
        8               0                 0      
        9               0                 0      
       10               0                 0      
       11               0                 0      
       12               0                 0      
       13               0                 0      
       14               0                 0      
       15               0                 0      
       16               0                 0      
       17               0                 0      
       18               0                 0      
       19               0                 0      
       20               0                 0      
       21               0                 0      
       22               0                 0      
       23               0                 0      
       24               0                 0      
       25               0                 0      
       26               0                 0      
       27               0                 0      
       28               0                 0      
       29               0                 0      
       30               0                 0      
       31               0                 0      
       32               0                 0      
       33               0                 0      
       34               0                 0      
       35               0                 0      
       36               0                 0      
       37               0                 0      
       38               0                 0      
       39               0                 0      
       40               0                 1      
       41               0                 0      
       42               0                 0      
       43               0                 0      
       44               0                 0      
       45               0                 0      
       46               0                 0      
       47               0                 0      
       48               0                 0      
       49               0                 0      
       50               0                 0      

The amplitude spectrum, is now easy to understand, it is simply the strengh of a and b at each frequency (figure 3).

 

Figure 3. Amplitude Spectrum of s

Second interlude: Imaginary numbers and rotations in the plane

If we take a real negative number and square it, say -2*-2 we obtain a positive value, here 4. This means that sqrt(4) has two solutions, 2 and -2. What if we look for the sqrt(-1)? Well instead of having no solution since a square number is supposed to be positive, we can ‘make one up’ using an imaginary number called i such as sqrt(-1)= i.

Because those new numbers are not on the number line of real numbers, this is typically represented as living orthogonally in the Euclidian space

Complex numbers

A complex number is a number that lives into the two-dimensional planes of real and imaginary numbers, for instance, 4+i and -2+2i. To add them, we proceed by adding values from each axis as we do with vectors of real values, i.e. 4-2=2 for the real part and 1i+2i=3i for the imaginary part (so the answer is 2+3i). 

What if we want to multiply them? Here something cool happens again.

Let’s step back for a second, and look at a coordinate (4,1). If I want to rotate this by -90 (𝝅/2) degrees on a plane, what will it be? (-1,4) and again (-4, -1) yes that’s a 180 degree, and again (1,-4) and once more (4,1) back to the initial position.

Let’s do a multiplication of 4+1i by i, ie i*(4+1i) = 4i+i^2 = -1+4i. That looks a lot like a -90 degree rotation in the plane. This means that multiplying by i is the same as rotating by -pi/4 and thus it is also totally logical to define imaginary numbers as living orthogonality to the real number line.  

Multiplication by i

Euler and infinite series

The Euler formula is: 

Here, e should not be understood as log ie 2.71828... but as an exponential. 

Oh this is so cool, this is an infinite series, which is how we can compute pi! ok, let's plug pi as our angle theta

and what now if we compute exp(1)

Last but not least, what if theta = i?

remember that multiplying by i is rotating by -pi/2, which means i^2 is pi or 180 degree, i^3 is 270 degrees, i^4 is 360 degrees (back to initial position), i^5 is the same as i, i^6 is the same as i^2, etc ... we simply rotate but because of the denominator the rotation or 'circle' is getting ever smaller at each quarter turn. 

Let's finish this feast, taking theta = pi as for our initial Euler formula

 

now instead of the circle getting smaller, the amount of rotation is getting smaller, and one keeps rotating on the unit circle! 

Another take on Euler formula (still rotating though)? watch this mind-blowing video from 3blue1brown

The complex Fourier transform

Thanks to the unit circle and trigonometry, for a given vector in the complex plane, we know the length in x (real numbers) is a*cos(theta) and the length in y (imaginary numbers) is b*sin(theta). 

If a wave is characterized by its cosine and sine, one can represent it as a location in the complex plane with cosine and sine basis (cos(x)+i*sin(x)). The length of the vector from (0,0) is the amplitude of the signal (cos+sin) and the temporal position of this wave in the original signal is given by the angle theta. 

Since multiplying by i, is equivalent to rotating in the plane, i*sin(x) is equivalent to moving the signal in time and the amount of shit in time is the amount of rotation (theta). How neat is that! We can finally understand the complex Fourier series notation. Since one computes this for each frequency, the time shift between waves at different frequencies is the difference in theta, also called phase.