Recursive and Iterative FFT
Algorithms for N = 2m
Objective:
At the end of this experiment, we will know about:
1. How to implement the recursive and iterative FFT algorithms?
2. Which of the two is more accurate?
3. How to calculate the runtime of a code block in Matlab?
4. What is the time complexity and run-time of the two algorithms,
as well as that of the DFT calculated using brute force (standard
formula)?
5. Why do we use FFT algorithms to calculate DFT?
6. What is circular convolution?
7. Why do we use circular convolution over linear convolution?
What is their time complexity?
Implementation:
a) Recursive FFT algorithm:
1. Define a recursive condition with return input(signal) if size of
input(signal) is one.
2. Then split the signal into two parts depending on even odd
nature of n.
3. Recursively perform FFT on both the even and odd parts.
4. Now, combine the two using the condition,
If k<n/2
X[k] = E[k] + (WN) k.O[k]
X[k+n/2] = E[k] - (WN) k.O[k]
b) Iterative FFT Algorithm:
1. Rearrange the input in bit-reversed order.
2 πk
2. Define the twiddle factors, W k M =e− j M .
3. Group the entries together by using butterflies of size 2,4,8,
……,N one after the while multiplying odd parts with the proper
twiddle factor.
c) Naive DFT Function:
1. Calculate DFT using the formula,
N −1 2π
−j kn
X [ k ]= ∑ x [n]e N
n=1
d) Inverse FFT algorithm:
1. In the previous iterative fft algorithm, change the twiddle factor
2 πk
to, W k M =e j M .
2. Normalize the output by dividing it by N.
e) Circular Convolution:
1. Pad zeros to the two signals( x 1 and x 2) to ensure that both
signals have length N, where N is the next power of two after (N 1
+N 2 −1), such that N 1 and N 2 are the lengths of the two signals.
2. Perform iterative fft on the new zero padded signals.
3. Multiply the two ffts.
4. Perform inverse fft on the signal obtained from multiplication.
f) Linear Convolution
1. Pad zeros to the two signals( x 1 and x 2) to ensure that both
signals have length N, where N is the next power of two after (N 1
+N 2 −1), such that N 1 and N 2 are the lengths of the two signals.
2. Perform linear convolution using the formula,
n
X [ n ]=∑ x1 (k )× x 2 (n−k )
k=0
g) Main script:
1. Calculate the DFT of a bandlimited signal using the 3 algorithms.
2. Calculate the error of the recursive and iterative algorithms
when N=64 using the formula,
√
N−1
∑ |X [ k ] −X mat [k ]|
2
k=1
Error= N−1
∑ |X mat|
2
k=1
3. Calculate the runtime of recursive FFT, iterative FFT, and naive
DFT for N=2 m , where m=3:1:10.
4. Plot the runtime against m for all 3 algorithms.
5. Calculate the runtime of linear and circular convolution for n,
where n=2 7 to 2 16 .
6. Calculate the error of linear and circular convolution.
Code:
1. Main Script:
clc;
clear all;
close all;
Fs=5000;% Sampling Frequency
runtime_navie=zeros(1,8);
runtime_recursive=zeros(1,8);
runtime_iterative=zeros(1,8);
neo=3:10;
for i=1:8 %
N=2^(i+2);
t=0:1/Fs:(N-1)/Fs;
x=3*sin(2*pi*500*t)+sin(2*pi*1500*t);
tic;
y_navie=(dft_navie(x,N));
runtime_navie(i)=toc;
% Naive DFT
tic;
y_recursive=(FFT_recursive(x,N));
runtime_recursive(i)=toc;
% Recursive FFT
tic;
y_iterative=(FFT_iterative(x,N));
runtime_iterative(i)=toc;
% Iterative FFT
y=fft(x,N);% FFT function in MATLAB
% Errors
error_navie=sqrt(sum((abs(y_navie-y)).^2)/sum((abs(y)).^2));
error_recursive=sqrt(sum((abs(y_recursive-y)).^2)/sum((abs(y)).^2));
error_iterative=sqrt(sum((abs(y_iterative-y)).^2)/sum((abs(y)).^2));
end
n=(-N/2:N/2-1)*Fs/N;
y_navie=fftshift(y_navie);
y_recursive=fftshift(y_recursive);
y_iterative=fftshift(y_iterative);
y=fftshift(y);
plot(t,x,LineWidth=2)
title('Signal')
xlabel('Time')
ylabel('Amplitude')
figure
plot(n,abs(y),LineWidth=2)
title('Matlab fft function')
xlabel('frequency')
ylabel('|X(f)|')
figure
plot(n,abs(y_navie),LineWidth=2)
title('Navie dft function')
xlabel('frequency')
ylabel('|X(f)|')
figure
plot(n,abs(y_recursive),LineWidth=2)
title('Recursive fft function')
xlabel('frequency')
ylabel('|X(f)|')
figure
plot(n,abs(y_iterative),LineWidth=2)
title('Iterative fft function')
xlabel('frequency')
ylabel('|X(f)|')
figure
plot(neo,runtime_navie,LineWidth=2)
title('Time Complexity of Naive DFT')
xlabel('m')
ylabel('Time')
grid on
figure
plot(neo,runtime_recursive,LineWidth=2)
title('Time Complexity of Recursive FFT')
xlabel('m')
ylabel('Time')
grid on
figure
plot(neo,runtime_iterative,LineWidth=2)
title('Time Complexity of Iterative FFT')
xlabel('m')
ylabel('Time')
grid on
figure
plot(neo,runtime_navie,LineWidth=2)
hold on
plot(neo,runtime_recursive,LineWidth=2)
hold on
plot(neo,runtime_iterative,LineWidth=2)
xlabel('m')
ylabel('Time')
title('Time Complexity')
legend('Time Complexity of Naive DFT','Time Complexity of Recursive
FFT','Time Complexity of Iterative FFT')
grid on
N=[128,256,512,1024,2048,4096,8192,16384,32768,65536];
lin_conv_runtime=zeros(1,10);
fft_conv_runtime=zeros(1,10);
for i=1:10
t=0:1/N(i):1-1/N(i);
x1=sin(2*pi*100*t);
x2=sin(2*pi*200*t);
tic;
X_f=conv_fft(x1,x2);
fft_conv_runtime(i)=toc;
tic;
X_l=lin_conv(x1,x2);
lin_conv_runtime(i)=toc;
end
plot(N,lin_conv_runtime,LineWidth=2);
hold on
plot(N,fft_conv_runtime,LineWidth=2);
xlabel('N=N1=N2')
ylabel('Run-time(in sec)')
legend('Direct Linear convolution','Circular Convolution')
title('Run-time vs N')
grid on
N1=25;
N2=30;
t1=0:1/N1:1-1/N1;
t2=0:1/N2:1-1/N2;
x1=sin(2*pi*10*t1);
x2=sin(2*pi*20*t2);
X_l=lin_conv(x1,x2);
X_f=conv_fft(x1,x2);
N=size(X_f,2);
n=0:N-1;
N_neo=size(X_l,2);
n_neo=0:N_neo-1;
figure
plot(n_neo,X_l,LineWidth=2);
xlabel('n');
ylabel('Amplitude')
hold on
plot(n,X_f,LineWidth=2);
xlabel('n');
ylabel('Amplitude')
title('Time Complexity')
legend('dft_naive','Recursive FFT','Iterative FFT')
grid on
error=(sum((abs(X_l-X_f)).^2)/sum((abs(X_l)).^2));
disp(error)
2. Recursive FFT function:
function d=FFT_recursive(x,n)
d=zeros(1,n);
if n==1
d=x(n);
else
ex=zeros(1,n/2);
ox=zeros(1,n/2);
for i=1:n
if rem(i-1,2)==0
ex((i+1)/2)=x(i);
else
ox((i)/2)=x(i);
end
end
e_j=FFT_recursive(ex,n/2);
o_j=FFT_recursive(ox,n/2);
for j=1:n/2
temp=o_j(j)*exp((-1i)*2*pi*(j-1)/n);
d(j)=e_j(j)+temp;
d(j+n/2)=e_j(j)-temp;
end
end
end
3. Iterative FFT function:
function x1=FFT_iterative(x,N)
x1=bitrevorder(x);
m=log2(N);
for s=1:m
M=2^s;
W=exp((-1i)*2*pi/M);
for k=0:M:N-1
w=1;
for p=0:M/2-1
t=w*x1(k+p+M/2+1);
u=x1(k+p+1);
x1(k+p+1)=t+u;
x1(k+p+M/2+1)=u-t;
w=w*W;
end
end
end
4. Navie DFT function:
function X=dft_navie(x,n)
X=zeros(1,n);
for i=1:n
for j=1:n
X(i)=X(i)+x(j)*exp((-1i)*2*pi*(i-1)*(j-1)/n);
end
end
5. Inverse Iterative FFT function:
function x1=iFFT_iterative(x,N)
x1=bitrevorder(x);
m=log2(N);
for s=1:m
M=2^s;
W=exp((1i)*2*pi/M);
for k=0:M:N-1
w=1;
for p=0:M/2-1
t=w*x1(k+p+M/2+1);
u=x1(k+p+1);
x1(k+p+1)=t+u;
x1(k+p+M/2+1)=u-t;
w=w*W;
end
end
end
x1=x1/N;
6. Circular Convolution function:
function d=conv_fft(x1,x2)
N1=size(x1,2);
N2=size(x2,2);
N=N1+N2-1;
m=ceil(log2(N));
N=2^m;
x1_=zeros(1,N-N1);
x2_=zeros(1,N-N2);
x1=[x1,x1_];
x2=[x2,x2_];
X1=FFT_iterative(x1,N);
X2=FFT_iterative(x2,N);
X=X1.*X2;
d=iFFT_iterative(X,N);
7. Linear Convolution function:
function d=lin_conv(x1,x2)
N1=size(x1,2);
N2=size(x2,2);
N=N1+N2-1;
d=zeros(1,N);
for i=0:N-1
for j=0:N1-1
if (i-j)>=0 && (i-j)<= (N2-1)
d(i+1)=d(i+1)+x1(j+1)*x2(i-j+1);
end
end
end
Results:
Graphs:
Original Signal:
Matlab FFT function:
Naive DFT:
Spectrum:
Run-time:
Recursive FFT:
Spectrum:
Runtime:
Iterative FFT:
Spectrum:
Runtime:
Run-Time Comparison:
Linear Convolution:
Runtime:
Circular Convolution:
Runtime:
Convolution Run-time:
Errors:
a) Mean Squared Error of dft_naive = 1.07226498677876 ×10 -28
b) Mean Squared Error of FFT_iterative = 2.87327893672804 ×10 -31
c) Mean Squared Error of FFT_recursive = 7.05819275328769 ×10 -32
d) Mean Squared Error between zero-padded circular convolution
and linear convolution = 8.80138890912222×10 -31
Interpretation:
In this experiment, we performed DFT using 3 different algorithms,
namely the DFT naive(Standard formula), Recursive FFT, and Iterative
FFT. We then compared the runtime and error of the 3 algorithms.
The DFT naive algorithm is just the simple implementation of the
standard formula for DFT. The recursive FFT works on the principle of
recursion using repeated function calls to calculate the DFT. The
iterative FFT uses loops in order to implement the butterfly operators.
Out of the 3, DFT naive has the largest runtime as it has time
complexity O(n 2 ). The runtime of both the recursive and iterative
algorithms is similar at smaller values of n, with the runtime of the
recursive algorithm being a bit smaller. At larger values of n, it is clear
that the iterative algorithm is significantly faster than the recursive
one. The theoretically calculated time complexity of both algorithms
is O(nlogn), though the hidden constant of the iterative algorithm is
much smaller than that of the recursive one, making it faster. From
the plot of runtimes of the two algorithms, we can see that the
theoretical time complexity holds for large values of n but fails for the
smaller values. This is because time complexity is calculated with the
assumption that n is large, which lets us ignore the lower powers of n
in our calculations. In our case, the runtime of the iterative FFT is
c1×nlogn + c 2 n (the bitrevorder function has time complexity O(n) ).
At small values of n(c 2 >c 1 logn), the time complexity is O(n), and when
it changes to O(nlogn), the runtime falls off, creating a maxima as
seen at n=2 7 . The same this happens with recursive FFT. Another
reason for the deviation in runtime is that for small vectors, practical
details like In-place computation, Nested Loop Structure, and Memory
Access patterns overshadow the time complexity formula.
Observing the errors of the 3 algorithms, the DFT naive has the largest
error. This is primarily because it requires more number of arithmetic
operations, which causes more rounding off error and leads to lower
numerical stability than the other algorithms. The errors of the
recursive and iterative algorithms are almost equal, with the recursive
algorithm having slightly less error than the iterative algorithm. Due
to this, in practical circumstances, both are considered to have almost
the same accuracy.
In practice, we usually use iterative FFT as recursive FFT needs more
memory. Also, recursive FFT can cause stack overhead at large values
of N.
In the second part of this experiment, we performed linear
convolution by taking the circular convolution of zero-padded signals.
This is done because if we directly perform linear convolution, it has
time complexity O(n 2 ), but the time complexity of circular convolution
is O(nlogn), making it much faster. This can be easily verified by
looking at the plot of runtime against n at large values of n. Circular
convolution uses iterative FFT, which is why it also has deviation at
smaller values of n.
If we calculate the error between direct linear convolution and
circular convolution with zero padding, it is of the order 10 -31 , which is
negligible. This is why we say that the circular convolution of zero-
padded signals is the same as the linear convolution of the two, with
the exception of some zeros being added at the end.
