The Fast and The Fourier
code-golfrestricted-complexity
Implement the Discrete Fourier Transform (DFT) for a sequence of any length using a Fast Fourier Transform algorithm (FFT). This may implemented as either a function or a program and the sequence can be given as either an argument or using standard input. A DFT has time complexity of O(n2) whereas a FFT has time complexity of O(n log n).
The algorithm will compute a result based on standard DFT in the forward direction. The input sequence has length n and consists of the complex values {x0, x1, ..., xn-1}. The output sequence will have the same length and consists of {y0, y1, ..., yn-1} is defined by the relation below.

Bluestein's algorithm
One algorithm that meets these requirements Bluestein's algorithm. It is a special case of the Chirp-Z transform and is able to compute the FFT for a sequence of any length n by transforming it in order to solve it as a cyclic convolution which can be solved with a time complexity of O(n log n).
Keep in mind that it is not required that you only use this algorithm in your implementation. If you know a better way, feel free to use it.
First, an identity is used to rewrite the initial DFT in a form where a convolution can easily be recognized.

You can obtain two sequences from this new form

which allow you to write the DFT as a convolution of two sequences.

Sample
Get the input
x = [1, 2, 3, 4, 5]
Get the length of the input
n = 5
Compute the 'a' sequence
a = [1, 1.618 - 1.176j, -2.427 - 1.763j, 3.236 + 2.351j, -4.045 + 2.939j]
Compute the 'b' sequence
b = [1, 0.809 + 0.588j, -0.809 + 0.588j, 0.809 - 0.588j, -0.809 - 0.588j]
Compute the convolution of 'a' and 'b' (using summation)
y[0] = (a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + a[3]*b[3] + a[4]*b[4]) / b[0]
= 15 / 1 = 15
y[1] = (a[1]*b[0] + a[0]*b[1] + a[2]*b[1] + a[3]*b[2] + a[4]*b[3]) / b[1]
= (-4.045 + 1.314j) / (0.809 + 0.588j) = -2.5 + 3.441j
y[2] = (a[2]*b[0] + a[1]*b[1] + a[3]*b[1] + a[0]*b[2] + a[4]*b[2]) / b[2]
= (1.545 - 2.127j) / (-0.809 + 0.588j) = -2.5 + 0.813j
y[3] = (a[3]*b[0] + a[2]*b[1] + a[4]*b[1] + a[1]*b[2] + a[0]*b[3]) / b[3]
= (-2.5 + 0.813j) / (0.809 - 0.588j) = -2.5 - 0.813j
y[4] = (a[4]*b[1] + a[3]*b[1] + a[2]*b[2] + a[1]*b[3] + a[0]*b[4]) / b[4]
= 4.253j / (-0.809 - 0.588j) = -2.5 - 3.441j
The Fourier tranform of x
y = [15, -2.5 + 3.441j, -2.5 + 0.813j, -2.5 - 0.813j, -2.5 - 3.441j]
Rules
- This is code-golf so the shortest solution wins.
- Builtins that compute FFT in forward or backward (also known as inverse) directions are not allowed.
- Builtins that compute the convolution are not allowed. (Most will have not been allowed by the previous rule as they use FFT internally.)
- Your solution must have time complexity of O(n log n) where n is the length of the input sequence.
- Floating-point inaccuracies will not be counted against you.
Test Cases
FFT([1, 1, 1, 1]) = [4, 0, 0, 0]
FFT([1, 0, 2, 0, 3, 0, 4, 0]) = [10, -2+2j, -2, -2-2j, 10, -2+2j, -2, -2-2j]
FFT([1, 2, 3, 4, 5]) = [15, -2.5+3.44j, -2.5+0.81j, -2.5-0.81j, -2.5-3.44j]
FFT([5-3.28571j, -0.816474-0.837162j, 0.523306-0.303902j, 0.806172-3.69346j, -4.41953+2.59494j, -0.360252+2.59411j, 1.26678+2.93119j] = [2, -3j, 5, -7j, 11, -13j, 17]
Related
- Compute the Discrete Fourier Transform - This contains some implementations for the standard DFT algorithm which has time complexity O(n2). You'll want to understand how to implement this before trying FFT.
- Too Fast, Too Fourier: FFT Code Golf - This previous challenge is the precursor to the current challenge here. Before, you only had to consider sequences where the length n was a power of 2 which allowed for simpler recursive implementations. The difference here is that you now have to implement an FFT algorithm that will work for sequences with any length.