4
\$\begingroup\$

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 is a previous challenge for the DFT, now let’s compute the DCT! Use the DCT‑II definition from Wikipedia:

$$ X_{k} = \sum_{n \, = \, 0}^{N - 1} x_{n} \, \cos \left[ \frac{\pi}{N} \, \left( n + \frac{1}{2} \right) \, k \right] \quad\quad k = 0, \, \ldots, N - 1 $$

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

Rules

  • This is 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.
Improve this question
\$\endgroup\$
12
  • 3
    \$\begingroup\$ Mind including an explanation on what DCT is? \$\endgroup\$ Commented May 13, 2017 at 0:32
  • 11
    \$\begingroup\$ Please add some test cases. \$\endgroup\$ Commented May 13, 2017 at 2:33
  • 1
    \$\begingroup\$ very low people even knows what is it, consider it is kind of alchemistry) its almost impossible to compute DCT, as is impossible to get Au from Hg ;-) \$\endgroup\$ Commented May 18, 2017 at 0:57
  • 1
    \$\begingroup\$ please, post test cases as you get ones, from working algorithm, i have none by hand and is too silly for using so hard maths. and we could discuss these cases, if that will be actually required ;-) \$\endgroup\$ Commented May 18, 2017 at 1:00
  • 4
    \$\begingroup\$ @xakepp35 how do we know if our program is correct, if we don't have at least two test cases to test it against? Specifically, as long as you produce the same results as what? \$\endgroup\$ Commented May 18, 2017 at 1:13

8 Answers 8

Reset to default
3
\$\begingroup\$

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.

Improve this answer
\$\endgroup\$
3
\$\begingroup\$

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)
Improve this answer
\$\endgroup\$
3
\$\begingroup\$

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

Improve this answer
\$\endgroup\$
0
2
\$\begingroup\$

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!

Improve this answer
\$\endgroup\$
4
  • \$\begingroup\$ Can you use .5 instead of 1/2? \$\endgroup\$ Commented Aug 9, 2017 at 11:21
  • \$\begingroup\$ D'oh, of course! \$\endgroup\$ Commented Aug 9, 2017 at 11:25
  • \$\begingroup\$ Oh of course, thanks a lot! \$\endgroup\$ Commented Sep 20, 2017 at 16:15
  • \$\begingroup\$ Suggest rows() instead of numel() \$\endgroup\$ Commented Oct 19, 2019 at 18:55
1
\$\begingroup\$

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.

Improve this answer
\$\endgroup\$
1
\$\begingroup\$

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.

Improve this answer
\$\endgroup\$
0
\$\begingroup\$

APL(NARS), 35 chars

{↑+/⍵×2○N÷⍨1p1×,/↑∘.×/1r2 1-⊂⍳N←≢⍵}

It would use the DCT-II definition from Wikipedia showed in the question too. It appears strange numbers in the function because APL start to count from 1 in arrays. One more readable one is this below even if a little more long:

{{+/⍵{X[⍵]×2○N÷⍨1p1×(⍺-1)×⍵-1r2}¨⍳N}¨⍳N←≢X←⍵}

test:

  f←{↑+/⍵×2○N÷⍨1p1×,/↑∘.×/1r2 1-⊂⍳N←≢⍵}
  f 1 3 2 4
10 ¯2.388955165 ¯6.061692394E¯16 ¯2.07192983 
  f 1 ¯1 ¯1 1 
0 ¯6.627740718E¯17 2.828427125 4.800234458E¯16 
  f 1 0 0 0 0 0 0 0 0 0 
1 0.9876883406 0.9510565163 0.8910065242 0.8090169944 0.7071067812 0.5877852523 0.4539904997 0.3090169944 0.156434465
Improve this answer
\$\endgroup\$
0
\$\begingroup\$

Python, 95 bytes

def f(x):N=len(x);R=range(N);return[sum((1j**2)**((n+.5)*k/N)*x[n] for n in R).real for k in R]

Attempt This Online!

Straight implementation of the formula, with the trick that cos(pi * ...) == (-1+0j)**(...) and that (1j**2) == (-1+0j) (not saving characters, I just like the look of the chained exponentiation).

I tried to lambda it, but the saved return comes back with interest for saving N and R as assignment expressions (97 bytes):

lambda x:(N:=len(x),R:=range(N),[sum((1j**2)**((n+.5)*k/N)*x[n] for n in R).real for k in R])[-1]
Improve this answer
\$\endgroup\$