Discrete Fourier Transform - a Linear Algebra Perspective

DISCRETE FOURIER TRANSFORM - A LINEAR ALGEBRA PERSPECTIVE

                        Harikrishnan NB

OBJECTIVE

  • In this tutorial we will see the interpretation of Fourier representations using the help of Linear algebra.
  • We will learn how to code Discrete Fourier Transform from scratch using MATLAB
  • We will see some of its application in solving differential equations.

INTRODUCTION

Any point (x0,y0) in a 2d plane needs two axis X-axis and Y-axis to represent it.

Fig 1. Cartesian Coordinate representation

When we move to linear algebra the point [x0, y0] is represented as a vector because it has both direction and magnitude. So to represent any point in a 2-d plane we need two axis. The point x R2  . Similarly  a point in 3d space x R3 needs 3-axis to represent it. This we can abstract to n- dimensional space where a point x Rn requires n - axis to represent the point. Let’s move back to our 2-d example provided in Fig 1. The Linear Algebra point of view of the 2d example is as follows

Fig 2. A vector represented as the linear combination of basis vectors

Here the black color x0  and y0 represents the point in the 2-d space and red color represents x0  and y0  represents the projection of the vector x = [x0, y0] to the x and y axis respectively. Here we represent x- axis as b1 and y-axis as b2. Where b1 and b2 are defined as follows:-

b1 and b2  are called basis. These vectors b1  and b2 are called basis because these are linearly independent and span the entire 2d space. This means our point [x0, y0] can be represented as some linear combination of basis vectors. The multiplying factor x0 and y0 are called coefficients.

Now the question is how do I get the right coefficient if I only have the vector [x0, y0] ? For this we will use the help of linear algebra.  Let's rewrite the equation in Fig 2 as follows

Fig 3. Matrix representation

Now we need to find the right coefficients x0 and y0 given the vector [x0, y0]

So we multiply the left side of B matrix (Basis vectors stacked as columns) by its inverse B-1 so as to get identity matrix. This in turn leads to multiplying the left side of x with  B-1. So the final matrix multiplication will be as follows:

        B-1Bc = Ic = c = B-1x

c = B-1x

But luckily in our case we choose the basis as orthonormal (i.e b1Tb2 = 0 =>dot product between b1 and b2 is zero and their norm is 1 (norm =1 =>)),  the inverse of an orthonormal matrix is its transpose. In reality our basis vectors  need not be orthonormal. It should only satisfy linear independence and should span the vector space.

What is the advantage of orthonormal basis?

 For an orthonormal matrix, its inverse is its transpose. The calculation of the inverse of matrix has a complexity of order O(N3). If we you use orthonormal basis the computational complexity is reduced to a large extent and results in faster and efficient computation. So let's write down the equations:

BTBc = Ic = c = BTx

c = BTx

Let's see how does the matrix equation for finding coefficients c = BTx  and matrix equation for representing the vector x as a linear combination of basis vectors x = Bc pictorially. For abstraction we denote the coefficients ‘c’ with red dots and vector x with black dots.

style="width: 452.72px; height: 291.50px; margin-left: 0.00px; margin-top: 0.00px; transform: rotate(0.00rad) translateZ(0px); -webkit-transform: rotate(0.00rad) translateZ(0px);" title="">

Fig 4. This figure shows the computation of coefficients.

This can also be interpreted as the projection of vector x  to different basis b1, b2, b3….bn  yielding n coefficients. This also tells how much of vector x  is present in each of the basis. Based on the coefficients value we can tell which basis is the most contributing factor. (Note : I have given red dots for coefficients. This does not mean all the coefficient values are same. They all are different) 

In signal processing how the signal and vector x are related?

In signal processing, our x vector is nothing but the sampled version of original signal x(t) with a finite duration  of T seconds. The moment the signal x(t) is sampled uniformly into N samples, we are moving into discrete domain.  Now the vector x has N samples so this vector x is a point in a N- dimensional space. So we need N basis vectors to find the coefficients. So the B matrix size will be N x N. This abstraction can be understood by the 2d example we took above. If vector x R2. Then B matrix has a size 2 x 2. This means we need 2 basis vectors.

How does the DFT and the example we took are related?

 Computing Fourier coefficients means we are computing the c.  The only difference in DFT is that the basis vectors are complex exponential. The reason why we took complex exponential as basis and how it is made will be addressed later. As mentioned above the complex exponential basis vectors in DFT are orthonormal( For computational easiness). So we can follow the same step as we mentioned earlier to find c. We project signal/ vector to the complex exponential basis to find the coefficients ‘c’.

Now lets see, given a signal / function

  1. Sample the signal
  2. What is the fundamental frequency of the signal?
  3. Define the DFT Basis
  4. Compute the Fourier Coefficients
  5. How to interpret the DFT plot

Fig 5. Sampling a signal

Since we sampled x(t) into N samples we will use linear algebra to represent the N samples as a vector x = [ x0 , x1, x2,......,xN].  Now lets visualize this.

Fig 6. Matrix vector representation of sampled signal

Now lets see how should we sample our x(t) and what is the criteria it has to follow. If we take only few samples then it gives a poor reconstruction of original signal x(t). So we need to take more number of samples. This is provided in the Shannon- Nyquist sampling criteria i.e, sampling frequency fs 2 * fm (fm is the maximum frequency of the signal x(t)). This can be even simplified as follows. If we take 1024 samples then the maximum frequency that we use to represent our sampled signal x is 1024/2 * fundamental frequency(f0) = 512*f0. If we choose the number of samples as 256 then the maximum frequency that we use to represent our sampled signal x is 256/2*f0 = 128*f0. So it is up to us to use our discretion to define the number of samples to be taken (always keep in mind it should satisfy shannon Nyquist sampling theorem criteria). The next important thing is the fundamental frequency. The fundamental frequency is given by 1/T where T is the duration  (see figure 5). The fundamental frequency has nothing to do with the number of samples. It is 1/(duration of signal).

The next important question is how do we define our basis?

Fig 7. Basis creation- Programming Intuition. Each row of the matrix B  is our basis 

Here the outer product gives a matrix and exponential is applied to each element of the matrix. The fig 7 is provides an easier way of implementing the B matrix in MATLAB.(Note: Fig 7 is an easier way to create Bmatrix programmatically. But mathematically we should not write like this).

Now where does the picture of fundamental frequency comes?

Till now we have not used the information of fundamental frequency. Refer figure 5 for fundamental frequency. Fundamental frequency comes while interpreting the results. If we simply plot the DFT coefficients we will not get these coefficients corresponds which frequency. To know this we need the information of fundamental frequency. So after computing c = Bx. (Note: Here I used c = Bx instead of BT x because already my basis functions defined in fig 7 are row wise so I don’t need to transpose it again. The key idea is we are projecting our sampled signal into different basis functions. In fig 7 the basis are already row wise). Similarly the Inverse Discrete Fourier Transform is computed as follows

x = BTc. Here the meaning is the vector x is represented as the linear combination of the basis vectors. In the case of BT  the columns represent the basis vectors(Note: Don’t be confused by the earlier notation and the current notation(B and BT). The reason is provided above with blue color)

Now for proper interpretation we have to plot the fourier coefficients with respect to w defined in fig 8. Now we will see how will plot  Fourier coefficients with respect to w.

style="width: 446.50px; height: 368.32px; margin-left: 0.00px; margin-top: 0.00px; transform: rotate(0.00rad) translateZ(0px); -webkit-transform: rotate(0.00rad) translateZ(0px);" title="">

Fig 8. w for interpretation

This is the real angular frequency corresponding to the DFT coefficients

 We need to  plot the coefficients with respect to w (defined in  fig 8) as per figure 9.

Fig 9. Plotting DFT coefficients with respect to their angular frequency

Here if we divide w (angular frequency) by 2*we get frequency.

LAB EXERCISE


MATLAB CODE for DFT from scratch

%% Author: Harikrishnan NB

% This code does DFT (Discrete Fourier Transform) from scratch

 

 

close all;

clear all;

clc;

t = 0+0.001:0.001:2;

x = sin(2*pi*13*t) + cos(2*pi *25*t) ; % Input Signal

%x = [zeros(1,975) ones(1,50) zeros(1,975)];% Give a square wave input  

 f0 = 1/t(end);% Fundamental Frequency

 N = length(t);

k1 = [0: N/2];

k2= [-(N-(N/2 + 1)):-1];

k = [k1, k2];

theta = [0:N-1]*(2*pi/N);

 

B = 1/sqrt(N)*exp(-i*k'*theta); % Basis

X = B*x';% X is the DFT Coeffcients

 subplot(311);

h = plot(t,x,'Linewidth',1.2);

grid;

xlabel('Time');

ylabel('Amplitude');

title('Original Signal');

 

subplot(312)

h = plot([k2,k1]*f0, abs([X(N/2 + 1 : N);X(1 : N/2)]),'Linewidth',1.1);

grid;

xlabel('Frequency f ');

ylabel('Magnitude Spectrum');

title('DFT');

 

% Inverse Discrete Fourier Transform

yhat = B'*X;

subplot(313);

h = plot(t, yhat,'Linewidth',1.2);

grid;

xlabel('Time');

ylabel('Amplitude');

title('IDFT');


How Fourier transform helps to compute the derivative of a function easily?

Find the derivative of sech(x) analytically and plot using MATLAB. Also find the derivative using the help of DFT.

If F(w) represents the Fourier Transform of f(t). Then the Fourier transform of f ‘(t)  (derivative of f(t)) is given by i*w*F(w).  Similarly the fourier transform of  nth order derivative of f(t) is given by (i*w)n * F(w). Now if we take the Inverse DFT of (i*w)n * F(w) we get . This helps to reduce the computational complexity.  Here w represents the angular frequency corresponding to a wave number ( w is defined in fig 8). So while doing computational experiments w should be treated as a vector. More clarity of this will be provided in fig 8. Let’s see the code. In the code I used k instead of w. Note that k is a vector.

MATLAB CODE 

Find the derivative of f(t) = sech(x) where t ranging from 0 to 2 analytically and plot using MATLAB. Also find the derivative using the help of DFT?

close all;

clear all;

clc;

 

L = 20;

n = 128;

x2 = linspace(-L/2,L/2,n+1);

x = x2(1:n);

u = sech(x);

ut = fft(u);

% DFT of f'(t) = i*k*DFT of f(t)  

% k represents the angular frequency

k = (2*pi/L) *[0 :(n/2 - 1) (-n/2):-1];

ut1 = i * k.*ut;

ut2 = - k.*k.*ut;

ut3 = -i * k.*k.*k.*ut;

u1 = ifft(ut1);

u2 = ifft(ut2);

u3 = ifft(ut3);

 

% u1exact is the analytical solution

u1exact = -sech(x) .* tanh(x);

figure(1);

plot(x,u1,'r',x,u1exact,'go');

title('Derivative computed using analytical method and DFT');

xlabel('time');

ylabel('amplitude');

grid;


Find the derivative of f(t) = sin(2*pi*13*t) + cos(2*pi*25*t) where t ranging from 0 to 2 analytically and plot using MATLAB. Also find the derivative using the help of DFT?

                                

                                

MATLAB CODE

close all;

clear all;

clc;

 

 

t = 0+0.001:0.001:2;

u = sin(2*pi*13*t) + cos(2*pi *25*t) ; % Input Signal

n = length(u);

%x = [zeros(1,975) ones(1,50) zeros(1,975)];% Give a square wave input  

f0 = 1/t(end);% Fundamental Frequency

ut = fft(u);

% DFT of f'(t) = i*k*DFT of f(t)  

% k represents the angular frequency

k = (2*pi*f0) *[0 :(n/2 - 1) (-n/2):-1];% Refer fig 8

ut1 = i * k.*ut;

ut2 = - k.*k.*ut;

ut3 = -i * k.*k.*k.*ut;

 

u1 = ifft(ut1);

u2 = ifft(ut2);

u3 = ifft(ut3);

 

% u1exact is the analytical solution

u1exact = 2*pi*13*cos(2*pi*13*t) - 2*pi*25*sin(2*pi*25*t);

figure(1);

plot(t,u1,'r',t,u1exact,'go');

title('Derivative computed using analytical method and DFT');

xlabel('time');

ylabel('amplitude');

grid;


Problem Solving

Find the first derivative of the following functions f(t) where t ranging from 0 to 2 analytically and plot using MATLAB. Also find the derivative using the help of DFT?

  1. f(t) = sin(2*pi*10*t) where 0 t 2
  2. f(t) = cos(2*pi*50*t) where 0 t 2
  3. f(t) = sin(2*pi*50*t2) where 0 t 2
  4. f(t) = sin(2*pi*20* t3 ) where 0 t 2
  5. f(t) = e2tsin(2*pi*5*t) where 0 t 2

Find the second derivative of the following functions f(t) where t ranging from 0 to 2 analytically and plot using MATLAB. Also find the derivative using the help of DFT?

  1. f(t) = sin(2*pi*10*t) where 0 t 2
  2. f(t) = cos(2*pi*50*t) where 0 t 2
  3. f(t) = sin(2*pi*50*t2) where 0 t 2
  4. f(t) = sin(2*pi*20* t3 ) where 0 t 2
  5. f(t) = e2tsin(2*pi*5*t) where 0 t 2


Solving differential equations using Discrete Fourier Transform

From previous section we learned to find derivative of function by using Fourier Transform. Now we will learn to solve a differential equation using Fourier Transform.

y’’ + a2 y = - f(t) --------------------  (1)

Equation (1) is a differential equation. So what is a differential equation? Differential equation is an equation involving the derivative of an unknown function. The solution of differential equation  in (1) is that unknown function y(t) that satisfies the equation (1).

Now using DFT we can solve this:

(iw)2 * Y(w) + a2 Y(w) =  - F(w)

(-w2 + a2) Y(w) = -F(w)

Y(w) = 

Now to find y(t), take the inverse fourier transform of Y(w). Here F(w) is known because F(w) is the Fourier Transform of known function/signal f(t). Also ‘a’ is a scalar and w is a vector.  We have seen w in detail in fig 8.

y(t) = Inverse Discrete Fourier Transform {}

Comments

  1. Always hated DSP stuff. I think this made me understand what numerous related courses of Engineering couldn't!

    ReplyDelete

Post a Comment

Popular posts from this blog