Compute the Discrete Cosine Transform

Question

Implement the Discrete Cosine Transform (DCT). This may implemented as either a function or a program and the sequence can be given as either an argument or using standard input. Your program must be able to handle sequences of any length (assuming a hypothetical version of your language which has no limits on things like memory and integer size).

There was a previous challenge for DFT, now lets compute DCT! ;-)

Use the DCT-II definition from Wikipedia:

Your program takes a sequence x_n as input, and must produce the corresponding sequence X_k. In this formula, the cosine is in radians.

Rules

• This is code-golf so the shortest solution wins.
• Builtins that compute the DCT in forward or backward (also known as inverse) directions are not allowed.
• Floating-point inaccuracies will not be counted against you.
• You don't have to use the exact algorithm shown here, as long as you produce the same results.

Answer 1

Mathematica, 79 bytes

Tr/@(s=#;x=Length@s;Table[s[[n+1]]Cos[Pi/x(n+1/2)#],{n,0,x-1}]&/@Range[0,x-1])&

Since there are (currently) no test cases, as it turns out, there is a built-in for this: FourierDCT. However, the function "normalizes" the result by dividing it by a value before returning it. Thankfully, the documentation specifies that this division factor is just the square root of the input list's length (for DCT-II).

So, we can verify our results by multiplying the output of FourierDCT by the square root of the length of our input list. For anyone else who tries this problem, here are some test cases:

FourierDCT@{1,3,2,4}*Sqrt@4

{10, -2.38896, 0, -2.07193}

FourierDCT@{1,1,1,1,1}*Sqrt@5

{5, 0, 0, 0, 0}

Note that I've slightly edited the output of the last one; the last four values output as very small decimals, around 10^-16, due to floating point precision errors around 0. My answer does give 0 for them, if you convert it to a decimal.

This solution gives exact answers when Mathematica can do so. If that's not okay and you want decimals, let me know.

Answer 2

Jelly, 15 bytes

LḶµ+.×þ÷L-*Ḟ³æ.

Try it online!

This uses the identity that cos(pi * x) = real(e^(i * pi * x)) = real((-1)^x)

Explanation

LḶµ+.×þ÷L-*Ḟ³æ.  Input: array A
L                Length
 Ḷ               Lowered range - [0, 1, ..., len(A)-1]
  µ              Begin a new monadic chain
   +.            Add 0.5
     ×þ          Outer product using multiplication
       ÷         Divide by
        L        Length
         -*      Compute (-1)^x
           Ḟ     Real part
             æ.  Dot product with
            ³    The first argument (A)

Answer 3

C (gcc), 88 83 78 77 bytes

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

Try it online!

This assumes the output vector is cleared before use.

Thanks to @Peter Taylor for -5

Answer 4

Clojure, 100 bytes

#(for[c[(count %)]k(range c)](apply +(map(fn[n x](*(Math/cos(*(/ c)Math/PI(+ n 0.5)k))x))(range)%)))

Sometimes it is difficult to be creative.

Answer 5

Octave, 49 46 45 42 bytes

Here we consider the DCT as a matrix multiplication. The matrix is basically the entry-wise cosine of a rank 1 matrix which is very simple to construct.

Thanks @Cowsquack for -1 byte, thanks @ceilingcat for another -3 bytes!

@(x)cos(pi/(n=numel(x))*(0:n-1)'*(.5:n))*x

Try it online!

Answer 6

Mathematica, 51 49 bytes

Re@Total[I^Array[(2##+#2)/n&,{n=Length@#,n},0]#]&

Try it online! (Using Mathics)

Based on my solution to DFT.