Digital Signal Processing/Fast Fourier Transform (FFT) Algorithm

The Fast 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: center|800 px|16-point DIF FFT circuit In the picture, the coefficients are given by: where N is the number of points of the FFT. They are evenly spaced on the lower half of the unit circle and . 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.