APL(NARS), 29 chars

{↑+/⍵÷0J1*N÷⍨4×,/∘.×⍨1-⍳N←≢⍵}

In above I would rewrite the formula in the post https://codegolf.stackexchange.com/a/82667/120560 One more readable one is this below even if a little more long:

{{+/⍵{X[⍵]÷0J1*N÷⍨4×(⍵-1)×⍺-1}¨⍳N}¨⍳N←≢X←⍵}

test:

  h←{↑+/⍵÷0J1*N÷⍨4×,/∘.×⍨1-⍳N←≢⍵}
  h 1 2 3 4 5
15J0 ¯2.5J3.440954801 ¯2.5J0.8122992406 ¯2.5J¯0.8122992406 ¯2.5J¯3.440954801 
  h 1 1 1 1
4J0 0J0 0J0 0J0 
  h 1 2 3 4
10J0 ¯2J2 ¯2J0 ¯2J¯2 

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

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.

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

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.

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.

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

Mathematica, 49 48 47 bytes

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

Based on the formula from @Dennis' solution.

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

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

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

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!

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.

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

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.