ELL319 Digital Signal Processing

Lab Report | Experiment 3

Vaishnavi Nandkishor Bisen

2022EE11719

March 7, 2025

DFT, IDFT, and FFT

1 Objectives
• To compute the N-point Discrete Fourier Transform (DFT) and Inverse Discrete Fourier
Transform (IDFT) of a given sequence.

• To analyze the magnitude and phase spectra of the transformed signal.

• To perform linear and circular convolution using DFT and IDFT.

• To implement and understand the Fast Fourier Transform (FFT) algorithm.

• To compare the computational efficiency of FFT and DFT.

2 Code
This is the code for all the parts, and the observation given is the output of the code.
1 #i n c l u d e <s t d i o . h>
2 #i n c l u d e <math . h>
3 #i n c l u d e < s t d l i b . h>
4 #i n c l u d e <time . h>
5 #d e f i n e PI 3 . 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6
6 typedef struct {
7 double r e a l ;
8 d o u b l e imag ;
9 } Complex ;
10 d o u b l e magnitude ( Complex c ) {
11 r e t u r n s q r t ( c . r e a l ∗ c . r e a l + c . imag ∗ c . imag ) ;
12 }
13 d o u b l e phase ( Complex c ) {
14 r e t u r n atan2 ( c . imag , c . r e a l ) ;
15 }
16 v o i d computeDFT ( d o u b l e x [ ] , Complex X [ ] , i n t N) {
17 int k , n;
18 f o r ( k = 0 ; k < N; k++) {
19 X[ k ] . r e a l = 0 ;
20 X[ k ] . imag = 0 ;
21 f o r ( n = 0 ; n < N; n++) {
22 d o u b l e a n g l e = −2.0 ∗ PI ∗ k ∗ n / N;
23 X[ k ] . r e a l += x [ n ] ∗ c o s ( a n g l e ) ;
24 X[ k ] . imag += x [ n ] ∗ s i n ( a n g l e ) ;
25 }

1
26 }
27 }
28 v o i d computeIDFT ( Complex X [ ] , d o u b l e x r e c o n s t r u c t e d [ ] , i n t N) {
29 int n , k ;
30 f o r ( n = 0 ; n < N; n++) {
31 double sum real = 0 . 0 ;
32 f o r ( k = 0 ; k < N; k++) {
33 d o u b l e a n g l e = 2 . 0 ∗ PI ∗ k ∗ n / N;
34 s u m r e a l += X[ k ] . r e a l ∗ c o s ( a n g l e ) − X[ k ] . imag ∗ s i n ( a n g l e ) ;
35 }
36 x r e c o n s t r u c t e d [ n ] = s u m r e a l / N;
37 }
38 }
39 Complex multiplyComplex ( Complex a , Complex b ) {
40 Complex r e s u l t ;
41 r e s u l t . r e a l = a . r e a l ∗ b . r e a l − a . imag ∗ b . imag ;
42 r e s u l t . imag = a . r e a l ∗ b . imag + a . imag ∗ b . r e a l ;
43 return r e s u l t ;
44 }
45 v o i d c o n v o l u t i o n ( d o u b l e x [ ] , d o u b l e h [ ] , d o u b l e y [ ] , i n t L , i n t K) {
46 int N = L + K − 1;
47 d o u b l e x padded [ 1 0 ] , h padded [ 1 0 ] ;
48 Complex X[ 1 0 ] , H[ 1 0 ] , Y [ 1 0 ] ;
49 int i ;
50 f o r ( i = 0 ; i < N; i ++) {
51 x padded [ i ] = ( i < L) ? x [ i ] : 0 . 0 ;
52 h padded [ i ] = ( i < K) ? h [ i ] : 0 . 0 ;
53 }
54 computeDFT ( x padded , X, N) ;
55 computeDFT ( h padded , H, N) ;
56 int k ;
57 f o r ( k = 0 ; k < N; k++) {
58 Y[ k ] = multiplyComplex (X[ k ] , H[ k ] ) ;
59 }
60 computeIDFT (Y, y , N) ;
61 }
62 v o i d convolutionTimeDomain ( d o u b l e x [ ] , i n t L , d o u b l e h [ ] , i n t K, d o u b l e y [ ] ) {
63 int N = L + K − 1;
64 int n;
65 f o r ( n = 0 ; n < N; n++) {
66 y [n ] = 0;
67 int k ;
68 f o r ( k = 0 ; k < L ; k++) {
69 i f ( n − k >= 0 && n − k < K) {
70 y [ n ] += x [ k ] ∗ h [ n − k ] ;
71 }
72 }
73 }
74 }
75 Complex addComplex ( Complex a , Complex b ) {
76 Complex r e s u l t ;
77 result . real = a . real + b. real ;
78 r e s u l t . imag = a . imag + b . imag ;
79 return r e s u l t ;
80 }
81 Complex subtractComplex ( Complex a , Complex b ) {
82 Complex r e s u l t ;
83 result . real = a . real − b. real ;
84 r e s u l t . imag = a . imag − b . imag ;
85 return r e s u l t ;
86 }
87 v o i d f f t ( Complex ∗X, i n t N) {
88 i f (N <= 1 ) r e t u r n ;

2
 89 Complex even [ 5 0 ] , odd [ 5 0 ] ;
90 int i ;
91 f o r ( i = 0 ; i < N / 2 ; i ++) {
92 even [ i ] = X[ i ∗ 2 ] ;
93 odd [ i ] = X[ i ∗ 2 + 1 ] ;
94 }
95 f f t ( even , N/ 2 ) ;
96 f f t ( odd , N/ 2 ) ;
97 int k ;
98 f o r ( k = 0 ; k < N / 2 ; k++) {
99 Complex t w i d d l e ;
100 d o u b l e a n g l e = −2.0 ∗ PI ∗ k / N;
101 twiddle . r e a l = cos ( angle ) ;
102 t w i d d l e . imag = s i n ( a n g l e ) ;
103 Complex temp = multiplyComplex ( t w i d d l e , odd [ k ] ) ;
104 X[ k ] = addComplex ( even [ k ] , temp ) ;
105 X[ k + N/ 2 ] = subtractComplex ( even [ k ] , temp ) ;
106 }
107 }
108 i n t main ( ) {
109 int N = 8;
110 int K = 3;
111 double x [ 8 ] = { 1 . 0 , 2 . 0 , 3 . 0 , 4 . 0 , 5 . 0 , 6 . 0 , 7 . 0 , 8 . 0 } ;
112 double h [ 3 ] = { 1 . 0 , 1 . 0 , 1 . 0 } ;
113 Complex X [ 8 ] , H [ 3 ] ;
114 clock t start dft = clock () ;
115 computeDFT ( x , X, N) ;
116 clock t end dft = clock () ;
117 d o u b l e t i m e d f t = ( ( d o u b l e ) ( e n d d f t − s t a r t d f t ) ) / CLOCKS PER SEC ;
118 computeDFT ( h , H, K) ;
119 double x r [ 8 ] ;
120 computeIDFT (X, x r , N) ;
121 d o u b l e M[ 8 ] , P [ 8 ] ;
122 int i ;
123 f o r ( i = 0 ; i < N; i ++) {
124 M[ i ] = magnitude (X[ i ] ) ;
125 P [ i ] = phase (X[ i ] ) ;
126 }
127 d o u b l e Mh[ 3 ] , Ph [ 3 ] ;
128 int k ;
129 f o r ( k = 0 ; k < K; k++) {
130 Mh[ k ] = magnitude (H[ k ] ) ;
131 Ph [ k ] = phase (H[ k ] ) ;
132 }
133 d o u b l e y c onv [ 1 0 ] ;
134 double y [ 1 0 ] ;
135 c o n v o l u t i o n ( x , h , y conv , N, K) ;
136 convolutionTimeDomain ( x , N, h , K, y ) ;
137 Complex X f f t [ 8 ] ;
138 int j ;
139 f o r ( j =0; j <N; j ++){
140 X f f t [ j ] . r e a l=x [ j ] ;
141 X f f t [ j ] . imag =0;
142 }
143 clock t s t a r t f f t = clock () ;
144 f f t ( X f f t ,N) ;
145 clock t end fft = clock () ;
146 d o u b l e t i m e f f t = ( ( d o u b l e ) ( e n d f f t − s t a r t f f t ) ) / CLOCKS PER SEC ;
147 return 0;
148 }
lab3.c

3
 Code for Part A and Part B and the graphs of this code is below.
1 #i n c l u d e <s t d i o . h>
2 #i n c l u d e <math . h>
3 #i n c l u d e < s t d l i b . h>
4
5 #d e f i n e PI 3 . 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6
6
7 typedef struct {
8 double r e a l ;
9 d o u b l e imag ;
10 } Complex ;
11
12 // t y p e d e f s t r u c t {
13 // d o u b l e magnitude ;
14 // d o u b l e phase ;
15 // } FC ;
16
17 d o u b l e magnitude ( Complex c ) {
18 r e t u r n s q r t ( c . r e a l ∗ c . r e a l + c . imag ∗ c . imag ) ;
19 }
20
21 d o u b l e phase ( Complex c ) {
22 r e t u r n atan2 ( c . imag , c . r e a l ) ;
23 }
24
25 v o i d computeDFT ( d o u b l e x [ ] , Complex X [ ] , i n t N) {
26 int k ;
27 f o r ( k = 0 ; k < N; k++) {
28 X[ k ] . r e a l = 0 ;
29 X[ k ] . imag = 0 ;
30 int n;
31 f o r ( n = 0 ; n < N; n++) {
32 d o u b l e a n g l e = −2.0 ∗ PI ∗ k ∗ n / N;
33 X[ k ] . r e a l += x [ n ] ∗ c o s ( a n g l e ) ;
34 X[ k ] . imag += x [ n ] ∗ s i n ( a n g l e ) ;
35 }
36 }
37 }
38
39 // v o i d computeFC ( Complex X [ ] , FC f [ ] , i n t N) {
40 // int j ;
41 // f o r ( j = 0 ; j < N; j ++){
42 // f [ j ] . magnitude = magnitude (X[ j ] ) ;
43 // f [ j ] . phase = phase (X[ j ] ) ;
44 // }
45 // }
46
47 v o i d computeIDFT ( Complex X [ ] , d o u b l e x r e c o n s t r u c t e d [ ] , i n t N) {
48 int n , k ;
49 f o r ( n = 0 ; n < N; n++) {
50 double sum real = 0 . 0 ;
51 f o r ( k = 0 ; k < N; k++) {
52 d o u b l e a n g l e = 2 . 0 ∗ PI ∗ k ∗ n / N;
53 s u m r e a l += X[ k ] . r e a l ∗ c o s ( a n g l e ) − X[ k ] . imag ∗ s i n ( a n g l e ) ;
54 }
55 x r e c o n s t r u c t e d [ n ] = s u m r e a l / N;
56 }
57 }
58
59 i n t main ( ) {
60 Complex X [ 5 ] ;
61 // FC f [ 5 ] ;
62 double y [ 5 ] = { 1 . 0 , 2 . 0 , 3 . 0 , 4 . 0 , 5 . 0 } ;

4
63 computeDFT ( y , X, 5 ) ;
64 // computeFC (X, f , 5 ) ;
65 d o u b l e M[ 5 ] ;
66 double P [ 5 ] ;
67 int a ;
68 f o r ( a = 0 ; a < 5 ; a++) {
69 M[ a ] = magnitude (X[ a ] ) ;
70 P [ a ] = phase (X[ a ] ) ;
71 }
72
73 double x r [ 5 ] ;
74 computeIDFT (X, x r , 5 ) ;
75 return 0;
76 }
lab3partanb.c

3 Theory
3.1 Discrete Fourier Transform (DFT)
The Discrete Fourier Transform (DFT) is a fundamental mathematical technique used to analyze
the frequency components of a discrete signal. Given a sequence x[n] of length N , the DFT is
defined as:
N
X −1
X(k) = x[n]e−j2πkn/N for k ∈ {0, 1, . . . , N − 1}
n=0
Expanding the exponential term into its real and imaginary components:
N −1   N −1  
X 2πkn X 2πkn
X(k) = x[n] cos −j x[n] sin
N N
n=0 n=0

In this equation: - The real part represents the correlation between the input signal x[n]
and a cosine function. - The imaginary part represents the correlation between x[n] and a sine
function.
In the provided C code, the function computeDFT() calculates the DFT.
The magnitude and phase of the DFT are calculated using:
p
|X(k)| = Re(X(k))2 + Im(X(k))2
 
Im(X(k))
Phase(X(k)) = arctan
Re(X(k))
In the code, the magnitude and phase are computed via:

M[a] = magnitude(X[a]);
P[a] = phase(X[a]);

Finally, to convert the frequency index k to an actual frequency in Hertz:

k
f = fs ·
N
where fs is the sampling frequency.
This implementation effectively demonstrates the DFT process and magnitude/phase ex-
traction.

5
3.2 Inverse DFT (IDFT)
The Inverse Discrete Fourier Transform (IDFT) is the inverse of the DFT operation and is
defined as follows:
N −1
1 X
xn = X(k) · ej2πkn/N ∀n ∈ {0, 1, . . . , N − 1}
N
k=0
This operation reconstructs the original time-domain signal x[n] from its frequency spectrum
X(k).
In the provided C code, the IDFT is implemented using the function computeIDFT.
This implementation accurately follows the IDFT mathematical formula, ensuring that the
reconstructed signal closely matches the original input sequence.

3.3 Convolution of Two Signals using DFT and IDFT

An interesting and useful property of the Fourier Transform is that convolution in the time
(or sample) domain is equivalent to multiplication in the frequency domain. This property
simplifies the convolution process significantly.
Given two signals x[n] and h[n], their convolution y[n] = x[n] ∗ h[n] has the corresponding
Fourier transform:

Y (k) = X(k) · H(k)

where X(k) and H(k) are the Fourier transforms of x[n] and h[n] respectively.
The steps to compute the convolution using DFT and IDFT are as follows:
• Pad both signals with zeros to make their size N = L + K − 1, where L and K are the
lengths of x[n] and h[n].
• Compute the N-point DFT of x[n] and h[n] using the function computeDFT().
• Compute Y (k) = X(k) · H(k) — an element-wise multiplication in the frequency domain.
• Compute the N-point IDFT of Y (k) using the function computeIDFT() to obtain y[n].
In addition to this method, I have also implemented convolution directly using the standard
convolution formula:
∞
X
y[n] = x[k]h[n − k]
k=−∞
This method was computed using the function convolutionTimeDomain().
By comparing these two approaches — convolution via DFT/IDFT and direct convolution
— the efficiency and accuracy of both methods can be analyzed.
The convolution via Fourier Transform is generally faster for longer sequences, while the
standard method is straightforward and intuitive.

3.4 Fast Fourier Transform (FFT)

The Discrete Fourier Transform (DFT) has a time complexity of O(N 2 ), which becomes in-
efficient for large signals such as audio recordings or datasets with thousands of samples. To
improve this, the Fast Fourier Transform (FFT) was developed as an efficient algorithm that
reduces the complexity to O(N log N ).
FFT utilizes the divide-and-conquer strategy, recursively breaking the original signal into
smaller parts until they reach trivial sizes. The most common FFT algorithm is the radix-2
Decimation In Time (DIT) FFT, which assumes the signal size N is a power of 2.

6
3.4.1 FFT Algorithm Explanation
The N-point DFT of a discrete time signal is given by:
N
X −1
X(k) = x[n]e−j2πkn/N ∀k ∈ {0, 1, . . . , N − 1}
n=0

The Radix-2 DIT algorithm splits this into two smaller DFTs:
N/2−1 N/2−1
X X
Xk = x2m e−j2π(2m)k/N + x2m+1 e−j2π(2m+1)k/N
m=0 m=0

From this step, factoring out common terms reveals that the FFT is efficiently calculated
by combining the DFT of the even-indexed and odd-indexed parts of the original signal:

Xk = Ek + e−j2πk/N Ok

Xk+N/2 = Ek − e−j2πk/N Ok
Where: - Ek is the DFT of the even-indexed part of x[n] - Ok is the DFT of the odd-indexed
part of x[n]

3.4.2 FFT Algorithm Steps (Pseudo-code)

The steps for the radix-2 FFT algorithm are as follows:

Algorithm 1 FFT Algorithm (Radix-2 DIT)

1: Function: fft(X, N)
2: if N = 1 then
3: Return single element (base case)
4: else
5: Divide the input into even-indexed and odd-indexed parts
6: Recursively compute FFT for both parts
7: for k = 0 to N/2 − 1 do
8: Compute Xk = Ek + e−j2πk/N Ok
9: Compute Xk+N/2 = Ek − e−j2πk/N Ok
10: end for
11: end if

3.4.3 Implementation in Code

In my code, the FFT is implemented using the function fft().
This implementation efficiently computes the DFT by leveraging the divide-and-conquer
strategy, significantly improving performance for larger signals. The FFT method is especially
useful in real-time applications where computational efficiency is critical.

4 Graphs
The graphs given below are for this observation-

7
 Figure 1: Obseravtion for the graphs

Figure 2: graph of y[n] i.e. the signal of which we will do DFT

8
Figure 3: graph of M[n] i.e. Magnitude of DFT of y[n]

Figure 4: graph of P[n] i.e. Phase of DFT of y[n]

9
 Figure 5: graph of xr [n] i.e. reconstructed signal from y[n]’s DFT, X[n]

5 Observations

10
11
6 Results
In this section, I present the outcomes of the implemented Fourier Transform methods, including
DFT, IDFT, convolution, and FFT. The following points summarize the results:

• First, I computed the Discrete Fourier Transform (DFT) of a given sequence using the
computeDFT() function.

• Next, I calculated the Inverse Discrete Fourier Transform (IDFT) using the computeIDFT()
function to reconstruct the original sequence. The reconstructed sequence matched the
original input, confirming the correctness of the DFT and IDFT implementation.

• I then performed convolution using two methods:

1. Convolution via DFT and IDFT using the convolution() function.

2. Direct convolution using the traditional formula y[n] = ∞
P
k=−∞ x[k]h[n − k] with the
function convolutionTimeDomain().

Both methods produced identical results, verifying the accuracy of the implemented con-
volution techniques.

12
 • To analyze efficiency, I computed the Fast Fourier Transform (FFT) using the fft()
function and compared its execution time with the DFT method.

• The results showed that the FFT method significantly reduced computation time com-
pared to the standard DFT, confirming FFT’s efficiency advantage, particularly for larger
data sizes.

The obtained results validate the correctness of the Fourier Transform implementations and
demonstrate the computational advantage of FFT over DFT.

13