# Compute the Discrete Fourier Transform

Implement the Discrete Fourier Transform (DFT) for a sequence of any length. This may implemented as either a function or a program and the sequence can be given as either an argument or using standard input.

The algorithm will compute a result based on standard DFT in the forward direction. The input sequence has length N and consists of [x(0), x(1), ..., x(N-1)]. The output sequence will have the same length and consists of [X(0), X(1), ..., X(N-1)] where each X(k) is defined by the relation below.

$$X_k\overset{\text{def}}{=}\sum_{n=0}^{N-1}x_n\cdot e^{-2\pi ikn/N}$$

## Rules

• This is so the shortest solution wins.
• Builtins that compute the DFT in forward or backward (also known as inverse) directions are not allowed.
• Floating-point inaccuracies will not be counted against you.

## Test Cases

DFT([1, 1, 1, 1]) = [4, 0, 0, 0]
DFT([1, 0, 2, 0, 3, 0, 4, 0]) = [10, -2+2j, -2, -2-2j, 10, -2+2j, -2, -2-2j]
DFT([1, 2, 3, 4, 5]) = [15, -2.5+3.44j, -2.5+0.81j, -2.5-0.81j, -2.5-3.44j]
DFT([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]


### Help

There was a previous challenge for finding the DFT using a FFT algorithm for sequences with lengths equal to a power of 2. You may find some tricks there that might help you here. Keep in mind that this challenge does not limit you to any complexity and also requires your solution to work for sequences of any length.

# Jelly, 16 15 bytes

LR’µ×'÷L-*²³÷S€


Try it online!

### How it works

LR’µ×'÷L-*²³÷S€  Main link. Argument [x(0), ..., x(N-1)].

L                Length; yield N.
R               Range; yield [1, ..., N].
’              Decrement; yield [0, ..., N-1].
µ             Begin a new, monadic chain. Argument: [0, ..., N-1]
×'           Spawned multiply [0, ..., N-1] with itself, yielding the matrix
of all possible products k×n.
÷L         Divide each product by N.
-*       Compute (-1)**(kn÷L) for each kn.
²      Square each result, computing (-1)**(2kn÷L).
³÷    Divide [x(0), ..., x(N-1)] by the results.
S€  Compute the sum for each row, i.e., each X(k).

• ninja'd :( Jun 11 '16 at 3:34

# Python 3, 77 bytes

lambda x,e=enumerate:[sum(t/1j**(4*k*n/len(x))for n,t in e(x))for k,_ in e(x)]


Test it on Ideone.

The code uses the equivalent formula • Wow, humongous figures. It's nice seeing the equivalent formulas that can allow for shorter code. Jun 11 '16 at 5:22

# J, 30 20 bytes

3 bytes thanks to @miles.

Uses the fact that e^ipi = -1.

The formula becomes X_k = sum(x_n / ((-1)^(2nk/N))).

+/@:%_1^2**/~@i.@#%#


### Usage

>> DCT =: +/@:%_1^2**/~@i.@#%#
>> DCT 1 2 3 4 5
<< 15 _2.5j3.44095 _2.5j0.812299 _2.5j_0.812299 _2.5j_3.44095


where >> is STDIN and << is STDOUT.

"Floating-point inaccuracies will not be counted against you."

# APL(Dyalog Unicode), 30 bytes SBCS

{⍵∘{+/⍺÷0j1*(4×⍵×⍳≢⍺)÷≢⍺}¨⍳≢⍵}


Try it on APLgolf!

# APL(Dyalog Unicode), 31 bytes SBCS

{⍵∘{+/⍺×*a÷⍨○0J¯2×⍵×⍳a←≢⍺}¨⍳≢⍵}


Try it on APLgolf!

A dfn submission, assumes ⎕IO←0.

Ungolfed:

X ← {+/⍺×*(○0J¯2×⍵×⍳≢⍺)÷≢⍺}
f ← {⍵∘X¨⍳≢⍵}


APL translated to Python:

from math import *

def f(sequence):
def X(n):
return sum(
sequence[b] # ⍺
* e **      # ×*(...)
(
pi*-2j          # ○0J¯2
* n             # ×⍵
* b             # ×⍳≢⍺
/ len(sequence) # ÷≢⍺
)
for b in range(len(sequence))
)
return [X(i) for i in range(len(sequence))] # ⍵∘X¨⍳≢⍵


First it binds the sequence to a dfn, and runs the resulting function on each element in the range from $$\0\$$ to $$\N\$$ -- {⍵∘{...}¨⍳≢⍵}.

In the inner dfn it sums -- +/ -- the sequence -- ⍺ -- times $$\e\$$ to the power of -- ×*(...) -- $$\\pi\$$ times $$\-2i\$$ -- ○0J¯2 -- times $$\k\$$ -- ⍵ -- times each $$\n\$$ from $$\0\$$ to $$\N\$$ -- ⍳≢⍵ -- divided by $$\N\$$ -- ÷≢⍺.

# MATL, 20 16 bytes

-1yn:qt!Gn/E*^/s


Input is a column vector, i.e. use semicolon as separator:

[1; 1; 1; 1]
[1; 0; 2; 0; 3; 0; 4; 0]
[1; 2; 3; 4; 5]
[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]


This uses the formula in Leaky Nun's answer, based on the facts that exp() = −1, and that MATL's power operaton with non-integer exponent produces (as in most programming languages) the principal branch result.

Try it online!

Due to Octave's weird spacing when displaying complex numbers, the real and imaginary parts in each entry of the output are separated by a space, as are consecutive entries. If that looks too ugly, add a ! at the end of the code (17 bytes) to have each entry of the output on a different line.

### Explanation

-1      % Push -1
y       % Get input implicitly. Push a copy below and one on top of -1
n:q     % Row vector [0 1 ... N-1] where N is implicit input length
t!      % Duplicate and transpose: column vector
Gn      % Push input length
/       % Divide
E       % Multiply by 2
*       % Multiply, element-wise with broadcast. Gives the exponent as a matrix
^       % Power (base is -1), element-wise. Gives a matrix
/       % Divide matrix by input (column vector), element-wise with broadcast
s       % Sum


# Pyth, 30

ms.e*b^.n1****c_2lQk.n0d.j0)QU


Test Suite

Very naive approach, just an implementation of the formula. Runs into various minor floating point issues with values which should be additive inverses adding to result in values that are slightly off of zero.

Oddly .j doesn't seem to work with no arguments, but I'm not sure if I'm using it correctly.

• Congrats on 10k!! Jun 11 '16 at 1:00

# Pyth, 18 bytes

Uses the fact that e^ipi = -1.

The formula becomes X_k = sum(x_n / ((-1)^(2nk/N))).

ms.ecb^_1*cyklQdQU


Test suite.

# Julia, 45 bytes

~=endof
!x=sum(x./im.^(4(r=0:~x-1).*r'/~x),1)


Try it online!

The code uses the equivalent formula ## Python 2, 78 bytes

l=input();p=1
for _ in l:print reduce(lambda a,b:a*p+b,l)*p;p*=1j**(4./len(l))


The polynomial is evaluated for each power p of 1j**(4./len(l)).

The expression reduce(lambda a,b:a*p+b,l) evaluates the polynomial given by l on the value p via Horner's method:

reduce(lambda a,b:a*10+b,[1,2,3,4,5])
=> 12345


Except, out input list is reversed, with constant term at the end. We could reverse it, but because p**len(l)==1 for Fourier coefficients, we can use a hack of inverting p and multiplying the whole result by p.

Some equal-length attempts:

l=input();i=0
for _ in l:print reduce(lambda a,b:(a+b)*1j**i,l,0);i+=4./len(l)

l=input();i=0
for _ in l:print reduce(lambda a,b:a*1j**i+b,l+);i+=4./len(l)


As a function for 1 byte more (79):

lambda l:[reduce(lambda a,b:a*1j**(i*4./len(l))+b,l+)for i in range(len(l))]


An attempt at recursion (80):

f=lambda l,i=0:l[i:]and[reduce(lambda a,b:(a+b)*1j**(i*4./len(l)),l,0)]+f(l,i+1)


Iteratively simulating reduce (80):

l=input();p=1;N=len(l)
exec"s=0\nfor x in l:s=s*p+x\nprint s*p;p*=1j**(4./N);"*N


# C (gcc), 86 78 bytes

k;d(a,b,c)_Complex*a,*b;{for(k=c*c;k--;)b[k/c]+=a[k%c]/cpow(1i,k/c*k%c*4./c);}


Try it online!

This assumes the output vector is zeroed out before use.

# Python 2, 89 bytes

Uses the fact that e^ipi = -1.

The formula becomes X_k = sum(x_n / ((-1)^(2nk/N))).

lambda a:[sum(a[x]/(-1+0j)**(x*y*2./len(a))for x in range(len(a)))for y in range(len(a))]


Ideone it!

• Post that as a separate answer! Jun 11 '16 at 4:09
• Alright, if you say so. Jun 11 '16 at 4:10

# Mathematica, 4948 47 bytes

Total[I^Array[4(+##-##-1)/n&,{n=Length@#,n}]#]&


Based on the formula from @Dennis' solution.

# Axiom, 81 bytes

g(x)==(#x<2=>x;[reduce(+,[x.j/%i^(4*k*(j-1)/#x)for j in 1..#x])for k in 0..#x-1])


using the formula someone post here. Results

(6) -> g([1,1,1,1])
(6)  [4,0,0,0]
Type: List Expression Complex Integer
(7) -> g([1,2,3,4])
(7)  [10,- 2 + 2%i,- 2,- 2 - 2%i]
Type: List Expression Complex Integer
(8) -> g([1,0,2,0,3,0,4,0])
(8)  [10,- 2 + 2%i,- 2,- 2 - 2%i,10,- 2 + 2%i,- 2,- 2 - 2%i]
Type: List Expression Complex Integer
(11) -> g([1,2,3,4,5])
(11)
5+--+4       5+--+3    5+--+2      5+--+
\|%i   + 5%i \|%i   - 4\|%i   - 3%i\|%i  + 2
[15, --------------------------------------------,
5+--+4
\|%i
5+--+4       5+--+3    5+--+2      5+--+
\|%i   + 3%i \|%i   - 5\|%i   - 2%i\|%i  + 4
--------------------------------------------,
5+--+4
\|%i
5+--+4       5+--+3    5+--+2      5+--+
\|%i   + 4%i \|%i   - 2\|%i   - 5%i\|%i  + 3
--------------------------------------------,
5+--+4
\|%i
5+--+4       5+--+3    5+--+2      5+--+
\|%i   + 2%i \|%i   - 3\|%i   - 4%i\|%i  + 5
--------------------------------------------]
5+--+4
\|%i
Type: List Expression Complex Integer
(12) -> g([1,2,3,4,5.])
(12)
[15.0, - 2.5 + 3.4409548011 779338455 %i, - 2.5 + 0.8122992405 822658154 %i,
- 2.5 - 0.8122992405 822658154 %i, - 2.5 - 3.4409548011 779338455 %i]
Type: List Expression Complex Float


# Octave, 43 39 bytes

@(x)j.^(-4*(t=0:(s=rows(x))-1).*t'/s)*x


Try it online!

Multiplies the input vector by the DFT matrix.

# Raku, 40 bytes

{^$_ .map: ((1.roots($_)X**-*)Z*$_).sum}  Try it online! • $_ is the input list.
• ^$_ is the sequence of integers from zero to one less than the size of the list. • .map: maps those indices over the following anonymous function, where the * which follows X**- assumes the value of each successive index. • 1.roots($_) returns the N complex roots of unity, where N is the size of the input list $_. • X** -* cross-exponentiates those roots with the negative of the current index into the input list. • Z*$_ zips-with-multiplication those powers of the roots of unity with the input list.
• .sum sums the zipped list.