Discrete Fourier Transform (DFT) Along With IDFT With MatLAB Code

BASIC

The discrete Fourier transform (DFT) converts a finite sequence of equally-spaced samples of a function into a same-length sequence of equally-spaced samples of the discrete-time Fourier transform (DTFT), which is a complex-valued function of frequency.

DFT of x(n) is defined by,

MATLAB CODE

To evaluate a DFT code sometimes values of x(n) may be given as sample value (i.e. x=[3, 4, 5, 6]). But most of the case it will be given as a sine function or cosine function or a combination of both. If the problem is given with sample values then it is ok. But if the problem is given with a function then we need a sampling frequency. That can be fixed or may be needed to take input. Here we are declaring it.

fs=8000;

In the function, we need the sampling time. We can find that by inversing the sampling frequency.

ts=1/fs;

We need the number of points for the function x(n). It is also the number of point of DFT. Here that is N.

N=8;

Now we need to find the index of x(n). That is the n.

n=0:N-1;

Consider there is a function is given. You should to input that here. I will compute the function below.

x=sin(2*pi*1000*n*ts) + 0.5*sin(2*pi*2000*n*ts + 3*pi/4)

Now, we need a variable X. That will hold the value of X(m). As it will sum continuously multiplied values so it cannot be empty. Declare it with 0.

X = zeros(N,1);

We need two loop for calculating DFT. One loop is nested in another loop. The outer loop will run from 1 to N. That is the index of X(m). The inner loop will run from 1 to N. That is the index of x(n). In inner loop the process of finding dft will run as the formula above.

for m = 1:N
    for n = 1:N
        X(m) = X(m) + x(n)*exp(-2j*pi*(n-1)*(m-1)/N);
    end
end

As the value of n has been converted into the loop from array to variable. It should be declare again as array with the same values above.

n=0:N-1;

Our DFT is done. Now we will plot that. We need to plot 6 items. SO divide it to 6 plots (3,2).

Plot x(n) on first subplot.

subplot(3,2,1);
stem(n,x)
title('x(n)')

Plot X(m) on second subplot.

subplot(3,2,2);
stem(n,real(X))
title('X(m)')

Amplitude Spectrum is the absolute value or modulus of the complex number (DFT result is complex). Show the amplitude spectrum on third subplot.

subplot(3,2,3);
stem(n,abs(X))
title('Amplitude Spectrum')

Power Spectrum is the square of amplitude spectrum and divided by the n. To show power spectrum we need to find it.

y = abs(X).^2/N;

Now show the power spectrum on fourth subplot.

subplot(3,2,4);
stem(n,y)
title('Power Spectrum')

Phase spectrum is the angle or argument of complex number (DFT result). The phase is by default in radian. Convert it to degree. Show the phase spectrum on fifth subplot.

subplot(3,2,5);
stem(n,angle(X)*180/pi)
title('Phase Spectrum')

To find inverse DFT use the build in function idft. Show the idft result on sixth subplot.

subplot(3,2,6);
stem(n,real(ifft(X)))
title('Inverse DFT')

FULL CODE

fs=8000;
ts=1/fs;
N=8;
n=0:N-1;
x=sin(2*pi*1000*n*ts) + 0.5*sin(2*pi*2000*n*ts + 3*pi/4);
X = zeros(N,1);
for m = 1:N
    for n = 1:N
        X(m) = X(m) + x(n)*exp(-2j*pi*(n-1)*(m-1)/N);
    end
end
n=0:N-1;
subplot(3,2,1);
stem(n,x)
title('x(n)')
subplot(3,2,2);
stem(n,real(X))
title('X(m)')
subplot(3,2,3);
stem(n,abs(X))
title('Amplitude Spectrum')
y = abs(X).^2/N;
subplot(3,2,4);
stem(n,y)
title('Power Spectrum')
subplot(3,2,5);
stem(n,angle(X)*180/pi)
title('Phase Spectrum')
subplot(3,2,6);
stem(n,real(ifft(X)))
title('Inverse DFT')

OUTPUT

I think your experience with this tutorial is quite good. If you have to seek anything you can comment below.

Thanks.