Decimation in frequency 16 Point FFT Python code | scipy FFT

# 16 POINT FFT with Bit reversed OUTPUT Decimation in Frequency
import matplotlib.pyplot as plt
import numpy as np
from scipy.fftpack import fft

# Input Coefficients
x0 = 1 + 1j
x1 = 2 + 1j
x2 = 1 - 2j
x3 = 2 - 1j
x4 = 2 + 3j
x5 = 3 + 2j
x6 = 1 + 3j
x7 = 3 + 1j
x8 = -3 + 3j
x9 = 3 - 3j
x10 = -1 - 1j
x11 = -3 - 3j
x12 = 3 - 3j
x13 = -1 - 1j
x14 = -3 - 3j
x15 = 1 + 1j
x = [x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15]

# Twiddle Factors
# 3rd stage
tc0 = 1.0000
tc1 = 0.9239 - 0.3827j
tc2 = 0.7071 - 0.7071j
tc3 = 0.3827 - 0.9239j
tc4 = 0.0000 - 1.0000j
tc5 = -0.3827 - 0.9239j
tc6 = -0.7071 - 0.7071j
tc7 = -0.9239 - 0.3827j

# Twiddle Factors
# 2nd stage
tb0 = 1.0000
tb1 = 0.7071 - 0.7071j
tb2 = 0.0000 - 1.0000j
tb3 = -0.7071 - 0.7071j

# 1st stage
ta0 = 1.0000
ta1 = 0.0000 - 1.0000j

# Stage 3 outputs
s3_0 = x0 + x8
s3_1 = x1 + x9
s3_2 = x2 + x10
s3_3 = x3 + x11
s3_4 = x4 + x12
s3_5 = x5 + x13
s3_6 = x6 + x14
s3_7 = x7 + x15
s3_8 = (x0 - x8) * tc0
s3_9 = (x1 - x9) * tc1
s3_10 = (x2 - x10) * tc2
s3_11 = (x3 - x11) * tc3
s3_12 = (x4 - x12) * tc4
s3_13 = (x5 - x13) * tc5
s3_14 = (x6 - x14) * tc6
s3_15 = (x7 - x15) * tc7

Stage3 = [s3_0, s3_1, s3_2, s3_3, s3_4, s3_5, s3_6, s3_7, s3_8, s3_9, s3_10, s3_11, s3_12, s3_13, s3_14, s3_15]

# Stage 2 outputs
s2_0 = s3_0 + s3_4
s2_1 = s3_1 + s3_5
s2_2 = s3_2 + s3_6
s2_3 = s3_3 + s3_7
s2_4 = (s3_0 - s3_4) * tb0
s2_5 = (s3_1 - s3_5) * tb1
s2_6 = (s3_2 - s3_6) * tb2
s2_7 = (s3_3 - s3_7) * tb3
s2_8 = s3_8 + s3_12
s2_9 = s3_9 + s3_13
s2_10 = s3_10 + s3_14
s2_11 = s3_11 + s3_15
s2_12 = (s3_8 - s3_12) * tb0
s2_13 = (s3_9 - s3_13) * tb1
s2_14 = (s3_10 - s3_14) * tb2
s2_15 = (s3_11 - s3_15) * tb3

Stage2 = [s2_0, s2_1, s2_2, s2_3, s2_4, s2_5, s2_6, s2_7, s2_8, s2_9, s2_10, s2_11, s2_12, s2_13, s2_14, s2_15]
# stage 1
s1_0 = s2_0 + s2_2
s1_1 = s2_1 + s2_3
s1_2 = (s2_0 - s2_2) * ta0
s1_3 = (s2_1 - s2_3) * ta1
s1_4 = s2_4 + s2_6
s1_5 = s2_5 + s2_7
s1_6 = (s2_4 - s2_6) * ta0
s1_7 = (s2_5 - s2_7) * ta1
s1_8 = s2_8 + s2_10
s1_9 = s2_9 + s2_11
s1_10 = (s2_8 - s2_10) * ta0
s1_11 = (s2_9 - s2_11) * ta1
s1_12 = s2_12 + s2_14
s1_13 = s2_13 + s2_15
s1_14 = (s2_12 - s2_14) * ta0
s1_15 = (s2_13 - s2_15) * ta1

# stage 0
s0_0 = s1_0 + s1_1
s0_1 = s1_0 - s1_1
s0_2 = s1_2 + s1_3
s0_3 = s1_2 - s1_3
s0_4 = s1_4 + s1_5
s0_5 = s1_4 - s1_5
s0_6 = s1_6 + s1_7
s0_7 = s1_6 - s1_7
s0_8 = s1_8 + s1_9
s0_9 = s1_8 - s1_9
s0_10 = s1_10 + s1_11
s0_11 = s1_10 - s1_11
s0_12 = s1_12 + s1_13
s0_13 = s1_12 - s1_13
s0_14 = s1_14 + s1_15
s0_15 = s1_14 - s1_15

Y = [s0_0, s0_1, s0_2, s0_3, s0_4, s0_5, s0_6, s0_7, s0_8, s0_9, s0_10, s0_11, s0_12, s0_13, s0_14, s0_15]
Y = [Y[0], Y[8], Y[4], Y[12], Y[2], Y[10], Y[6], Y[14], Y[1], Y[9], Y[5], Y[13], Y[3], Y[11], Y[7], Y[15]]
# FFT using built-in python function
Y_fft = fft(x, 16)
Y = [abs(ele) for ele in Y]
Y_fft = [abs(ele) for ele in Y_fft]

Y1 = np.subtract(Y, Y_fft)

plt.plot(Y, 'r*') # Our 16 point FFT output plot
plt.title("16 point FFT using decimation in frequency")
plt.show()
plt.plot(Y_fft, 'y^') # Built in 16 point FFT output plot
plt.title("16 point FFT using built-in scipy FFT function")
plt.show()
plt.plot(Y1, 'g') # Difference between two outputs
plt.title("Difference between above two methods")
plt.show()