# Cythonized version of FDCT (fast discrete cosine transform) function, ported from Java

This is my Cython code for an FDCT function (invoked here):

from __future__ import division

import numpy
cimport numpy
cimport cython

FLOAT64 = numpy.float64

ctypedef numpy.int_t INT_t
ctypedef numpy.float64_t FLOAT64_t

cdef inline FLOAT64_t **cosine = [
[
3.535534e-01, 3.535534e-01, 3.535534e-01, 3.535534e-01,
3.535534e-01, 3.535534e-01, 3.535534e-01, 3.535534e-01
], [
4.903926e-01, 4.157348e-01, 2.777851e-01, 9.754516e-02,
-9.754516e-02, -2.777851e-01, -4.157348e-01, -4.903926e-01
], [
4.619398e-01, 1.913417e-01, -1.913417e-01, -4.619398e-01,
-4.619398e-01, -1.913417e-01, 1.913417e-01, 4.619398e-01
], [
4.157348e-01, -9.754516e-02, -4.903926e-01, -2.777851e-01,
2.777851e-01, 4.903926e-01, 9.754516e-02, -4.157348e-01
], [
3.535534e-01, -3.535534e-01, -3.535534e-01, 3.535534e-01,
3.535534e-01, -3.535534e-01, -3.535534e-01, 3.535534e-01
], [
2.777851e-01, -4.903926e-01, 9.754516e-02, 4.157348e-01,
-4.157348e-01, -9.754516e-02, 4.903926e-01, -2.777851e-01
], [
1.913417e-01, -4.619398e-01, 4.619398e-01, -1.913417e-01,
-1.913417e-01, 4.619398e-01, -4.619398e-01, 1.913417e-01
], [
9.754516e-02, -2.777851e-01, 4.157348e-01, -4.903926e-01,
4.903926e-01, -4.157348e-01, 2.777851e-01, -9.754516e-02
]
]

cdef inline FLOAT64_t ZERO = <FLOAT64_t>0.0
cdef inline INT_t EIGHT = <INT_t>8
cdef inline INT_t ROUND_FACTOR(FLOAT64_t x):
return <INT_t>(<FLOAT64_t>x + <FLOAT64_t>0.499999)

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def FDCT(numpy.ndarray[INT_t, ndim=1] shape_plane not None):
cdef int i, j, k
cdef FLOAT64_t s = ZERO
cdef numpy.ndarray[FLOAT64_t, ndim=1] dct = numpy.zeros_like(
shape_plane, dtype=FLOAT64)

for i in range(EIGHT):
for j in range(EIGHT):
s = ZERO
for k in range(EIGHT):
s += cosine[j][k] * <FLOAT64_t>shape_plane[EIGHT * i + k]
dct[EIGHT * i + j] = s

for i in range(EIGHT):
for j in range(EIGHT):
s = ZERO
for k in range(EIGHT):
s += cosine[i][k] * dct[EIGHT * k + j]
shape_plane[EIGHT * i + j] = ROUND_FACTOR(s)


I'm not really a domain-expert w/r/t this type of signal-processing code. I ported this function from a Java implementation:

private static void Fdct(int[] shapes) {
int i, j, k;
double s;
double[] dct = new double[64];

//calculation of the cos-values of the second sum
for (i = 0; i < 8; i++) {
for (j = 0; j < 8; j++) {
s = 0.0;
for (k = 0; k < 8; k++)
s += arrayCosin[j][k] * shapes[8 * i + k];
dct[8 * i + j] = s;
}
}

for (j = 0; j < 8; j++) {
for (i = 0; i < 8; i++) {
s = 0.0;
for (k = 0; k < 8; k++)
s += arrayCosin[i][k] * dct[8 * k + j];
shapes[8 * i + j] = (int) Math.floor(s + 0.499999);
}
}
}


My goal is to match the output of the original Java code, warts and all (if any) – any feedback is welcome, however.

• Can you explain why your goal is to "match the output of the original Java code, warts and all"? What is the underlying requirement here? What's wrong with scipy.fftpack.dct? – Gareth Rees Jul 1 '14 at 9:16
• @GarethRees I would totally use SciPy’s DCT function (or the one from integer FFTW, or somesuch) but I am trying to make this project produce results that produce something as close as possible to what the Java code yields (that’s what I mean by “warts and all”, incedentally). – fish2000 Jan 10 '15 at 20:53