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.

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

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.

Community · Accepted Answer · 2020-06-17 09:04:33Z

8

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

formula

edited Jun 17, 2020 at 9:04

CommunityBot

1

answered Jun 11, 2016 at 4:16

Dennis

210k41 gold badges368 silver badges819 bronze badges

\$\begingroup\$ Wow, humongous figures. It's nice seeing the equivalent formulas that can allow for shorter code. \$\endgroup\$
– miles
Jun 11, 2016 at 5:22

Add a comment |

Dennis · Accepted Answer · 2016-06-11 03:58:12Z

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

\$\begingroup\$ ninja'd :( \$\endgroup\$
– Leaky Nun
Jun 11, 2016 at 3:34 — Leaky Nun, Jun 11, 2016 at 3:34

Community · Accepted Answer · 2017-04-13 12:39:05Z

4

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

edited Apr 13, 2017 at 12:39

CommunityBot

1

answered Jun 11, 2016 at 1:45

Leaky Nun

49.6k6 gold badges107 silver badges281 bronze badges

Add a comment |

Ruan · Accepted Answer · 2021-05-25 15:31:20Z

APL(Dyalog Unicode), 30 bytes ^SBCS

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

Adapted from the formula presented in this Python answer by Dennis

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\$ -- ÷≢⍺.

Luis Mendo · Accepted Answer · 2021-02-17 00:34:41Z

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(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 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

FryAmTheEggman · Accepted Answer · 2016-06-11 00:46:26Z

2

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.

edited Jun 11, 2016 at 0:46

answered Jun 11, 2016 at 0:35

FryAmTheEggman

17.5k3 gold badges39 silver badges97 bronze badges

1

\$\begingroup\$ Congrats on 10k!! \$\endgroup\$
– Luis Mendo
Jun 11, 2016 at 1:00

Add a comment |

Leaky Nun · Accepted Answer · 2016-06-11 03:07:26Z

2

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.

answered Jun 11, 2016 at 3:07

Leaky Nun

49.6k6 gold badges107 silver badges281 bronze badges

Add a comment |

Dennis · Accepted Answer · 2016-06-11 06:39:19Z

2

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

formula

answered Jun 11, 2016 at 6:39

Dennis

210k41 gold badges368 silver badges819 bronze badges

Add a comment |

Community · Accepted Answer · 2020-06-17 09:04:33Z

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

ceilingcat · Accepted Answer · 2017-08-13 09:32:45Z

2

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.

edited Aug 13, 2017 at 9:32

answered Aug 9, 2017 at 7:47

ceilingcat

7,4572 gold badges28 silver badges30 bronze badges

Add a comment |

Community · Accepted Answer · 2020-06-17 09:04:33Z

1

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!

edited Jun 17, 2020 at 9:04

CommunityBot

1

answered Jun 11, 2016 at 3:49

Leaky Nun

49.6k6 gold badges107 silver badges281 bronze badges

\$\begingroup\$ Post that as a separate answer! \$\endgroup\$
– Leaky Nun
Jun 11, 2016 at 4:09
\$\begingroup\$ Alright, if you say so. \$\endgroup\$
– Dennis
Jun 11, 2016 at 4:10

Add a comment |

Community · Accepted Answer · 2017-04-13 12:39:05Z

1

Mathematica, 49 48 47 bytes

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

Based on the formula from @Dennis' solution.

edited Apr 13, 2017 at 12:39

CommunityBot

1

answered Jun 17, 2016 at 10:02

miles

16.9k2 gold badges28 silver badges94 bronze badges

Add a comment |

Community · Accepted Answer · 2020-06-17 09:04:33Z

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

ceilingcat · Accepted Answer · 2019-10-19 17:26:25Z

1

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.

edited Oct 19, 2019 at 17:26

answered Oct 19, 2019 at 17:06

ceilingcat

7,4572 gold badges28 silver badges30 bronze badges

Add a comment |

Sean · Accepted Answer · 2021-05-25 12:13:11Z

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.

Stack Exchange Network

Compute the Discrete Fourier Transform

Rules

Test Cases

Help

15 Answers 15

Python 3, 77 bytes

Jelly, 16 15 bytes

How it works

J, 30 20 bytes

Usage

APL(Dyalog Unicode), 30 bytes ^SBCS

APL(Dyalog Unicode), 31 bytes ^SBCS

MATL, 20 16 bytes

Explanation

Pyth, 30

Pyth, 18 bytes

Julia, 45 bytes

Python 2, 78 bytes

C (gcc), 86 78 bytes

Python 2, 89 bytes

Mathematica, 49 48 47 bytes

Axiom, 81 bytes

Octave, 43 39 bytes

Raku, 40 bytes

Your Answer

Not the answer you're looking for? Browse other questions tagged
code-golf
math
or ask your own question.

Linked

Hot Network Questions

Compute the Discrete Fourier Transform

Rules

Test Cases

Help

15 Answers 15

Python 3, 77 bytes

Jelly, 16 15 bytes

How it works

J, 30 20 bytes

Usage

APL(Dyalog Unicode), 30 bytes SBCS

APL(Dyalog Unicode), 31 bytes SBCS

MATL, 20 16 bytes

Explanation

Pyth, 30

Pyth, 18 bytes

Julia, 45 bytes

Python 2, 78 bytes

C (gcc), 86 78 bytes

Python 2, 89 bytes

Mathematica, 49 48 47 bytes

Axiom, 81 bytes

Octave, 43 39 bytes

Raku, 40 bytes

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged code-golfmath or ask your own question.

Linked

Related

Hot Network Questions

APL(Dyalog Unicode), 30 bytes ^SBCS

APL(Dyalog Unicode), 31 bytes ^SBCS

Not the answer you're looking for? Browse other questions tagged
code-golf
math
or ask your own question.