3
\$\begingroup\$

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.

\$\endgroup\$
  • 3
    \$\begingroup\$ 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? \$\endgroup\$ – Gareth Rees Jul 1 '14 at 9:16
  • \$\begingroup\$ @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). \$\endgroup\$ – fish2000 Jan 10 '15 at 20:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.