Digital Signal Processing/Fast Fourier Transform (FFT) Algorithm

The Fourier Transform (FFT) is an algorithm used to calculate a Discrete Fourier Transform (DFT) using a smaller number of multiplications.

Cooley–Tukey algorithm

The Cooley–Tukey algorithm is the most common fast Fourier transform (FFT) algorithm. It breaks down a larger DFT into smaller DFTs. The best known use of the Cooley–Tukey algorithm is to divide a N point transform into two N/2 point transforms, and is therefore limited to power-of-two sizes.

The algorithm can be implemented in two different ways, using the so-called:

• Decimation In Time (DIT)
• Decimation In Frequency (DIF)

Circuit diagram

The following picture represents a 16-point DIF FFT circuit diagram:

16-point DIF FFT circuit

In the picture, the coefficients are given by:

$W^{k}=e^{-j\cdot k/N\cdot 2\pi }$

where N is the number of points of the FFT. They are evenly spaced on the lower half of the unit circle and $W^{0}=1$.

Test script

The following python script allows to test the algorithm illustrated above:

#!/usr/bin/python3
import math
import cmath

# ------------------------------------------------------------------------------
# Parameters
#
coefficientBitNb = 8    # 0 for working with reals

# ------------------------------------------------------------------------------
# Functions
#
# ..............................................................................

def  bit_reversed(n, bit_nb):
    binary_number = bin(n)
    reverse_binary_number = binary_number[-1:1:-1]
    reverse_binary_number = reverse_binary_number + \
        (bit_nb - len(reverse_binary_number))*'0'
    reverse_number = int(reverse_binary_number, bit_nb-1)

    return(reverse_number)

# ..............................................................................

def  fft(x):
    j = complex(0, 1)
                                                                         # sizes
    point_nb = len(x)
    stage_nb = int(math.log2(point_nb))
                                                                       # vectors
    stored = x.copy()
    for index in range(point_nb):
        stored[index] = complex(stored[index], 0)
    calculated = [complex(0, 0)] * point_nb
                                                                  # coefficients
    coefficients = [complex(0, 0)] * (point_nb//2)
    for index in range(len(coefficients)):
        coefficients[index] = cmath.exp(-j * index/point_nb * 2*cmath.pi)
        coefficients[index] = coefficients[index] * 2**coefficientBitNb
                                                                          # loop
    for stage_index in range(stage_nb):
        # print([stored[i].real for i in range(point_nb)])
        index_offset = 2**(stage_nb-stage_index-1)
                                                           # butterfly additions
        for vector_index in range(point_nb):
            isEven = (vector_index // index_offset) % 2 == 0
            if isEven:
                operand_a = stored[vector_index]
                operand_b = stored[vector_index + index_offset]
            else:
                operand_a = - stored[vector_index]
                operand_b = stored[vector_index - index_offset]
            calculated[vector_index] = operand_a + operand_b
                                                   # coefficient multiplications
        for vector_index in range(point_nb):
            isEven = (vector_index // index_offset) % 2 == 0
            if not isEven:
                coefficient_index = (vector_index % index_offset) \
                    * (stage_index+1)
                # print(                                                  \
                #     "(%d, %d) -> %d"                                    \
                #     % (stage_index, vector_index, coefficient_index)    \
                # )
                calculated[vector_index] = \
                    coefficients[coefficient_index] * calculated[vector_index] \
                    / 2**coefficientBitNb
                if coefficientBitNb > 0:
                    calculated[vector_index] = complex(               \
                        math.floor(calculated[vector_index].real),    \
                        math.floor(calculated[vector_index].imag)     \
                    )
                                                                       # storage
        stored = calculated.copy()
                                                               # reorder results
    for index in range(point_nb):
        calculated[bit_reversed(index, stage_nb)] = stored[index]

    return calculated

# ------------------------------------------------------------------------------
# Main program
#
source = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
transformed = fft(source)
amplitude = [0.0] * len(source)
amplitude_print = ''
for index in range(len(transformed)):
    amplitude[index] = abs(transformed[index])
    amplitude_print = amplitude_print + "%5.3f " % amplitude[index]
print()
print(amplitude_print)

It has a special parameter, coefficientBitNb, which allows to determine the calculation results for a circuit using only integer numbers.