Compute the Discrete Fourier Transform

Question

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.

Rules

• This is code-golf 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.

Answer 1

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).

Answer 2

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

Answer 3

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."

Answer 4

MATL, 20 16 bytes

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

Input is a column vector, i.e. replace commas by semicolons:

[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(iπ) = −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 with complex numbers, the real and imaginary parts are separated by a space, as are the different entries of the resulting vector. If that looks too ugly, add a ! at the end (17 bytes) to have each entry of the output in a different row.

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

Answer 5

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.

Answer 6

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.

Answer 7

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.

Answer 8

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

Answer 9

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+[0]);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+[0])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

Answer 10

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

Answer 11

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!

Answer 12

Mathematica, 49 48 47 bytes

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

Based on the formula from @Dennis' solution.

Answer 13

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.