sysy-data/hidden_functional_c/sy/37_dct.sy

#include "sylib.h"

#include <string.h>
// reference:
// https://stackoverflow.com/questions/8310749/discrete-cosine-transform-dct-implementation-c

// const int MAX_DIM_X = 8, MAX_DIM_Y = 8;

#define MAX_DIM_X 8
#define MAX_DIM_Y 8

float test_block[MAX_DIM_X][MAX_DIM_Y];
float test_dct[MAX_DIM_X][MAX_DIM_Y];
float test_idct[MAX_DIM_X][MAX_DIM_Y];

const float PI = 3.14159265359, TWO_PI = 6.28318530718, EPSILON = 1e-6;

float my_fabs(float x)
{
    if (x > 0)
        return x;
    return -x;
}

float p(float x) { return 3 * x - 4 * x * x * x; }

float my_sin_impl(float x)
{
    if (my_fabs(x) <= EPSILON)
        return x;
    return p(my_sin_impl(x / 3.0));
}

float my_sin(float x)
{
    if (x > TWO_PI || x < -TWO_PI) {
        int xx = x / TWO_PI;
        x = x - xx * TWO_PI;
    }
    if (x > PI)
        x = x - TWO_PI;
    if (x < -PI)
        x = x + TWO_PI;
    return my_sin_impl(x);
}

float my_cos(float x) { return my_sin(x + PI / 2); }

void write_mat(float mat[][MAX_DIM_Y], int n, int m)
{
    int i = 0;
    while (i < n) {
        putfloat(mat[i][0]);
        int j = 1;
        while (j < m) {
            putch(32);
            putfloat(mat[i][j]);
            j = j + 1;
        }
        putch(10);
        i = i + 1;
    }
    putch(10);
}

void dct(float dct_mat[][MAX_DIM_Y], float mat[][MAX_DIM_Y], int n, int m)
{
    int u = 0;
    while (u < n) {
        int v = 0;
        while (v < m) {
            dct_mat[u][v] = 0;
            int i = 0;
            while (i < n) {
                int j = 0;
                while (j < m) {
                    dct_mat[u][v] = dct_mat[u][v] + mat[i][j] * my_cos(PI / n * (i + 1. / 2.) * u) * my_cos(PI / m * (j + 1. / 2.) * v);
                    j = j + 1;
                }
                i = i + 1;
            }
            v = v + 1;
        }
        u = u + 1;
    }
}

void idct(float mat[][MAX_DIM_Y], float dct_mat[][MAX_DIM_Y], int n, int m)
{
    int u = 0;
    while (u < n) {
        int v = 0;
        while (v < m) {
            mat[u][v] = 1 / 4. * dct_mat[0][0];
            int i, j;

            i = 1;
            while (i < n) {
                mat[u][v] = mat[u][v] + 1 / 2. * dct_mat[i][0];
                i = i + 1;
            }
            j = 1;
            while (j < m) {
                mat[u][v] = mat[u][v] + 1 / 2. * dct_mat[0][j];
                j = j + 1;
            }

            i = 1;
            while (i < n) {
                j = 1;
                while (j < m) {
                    mat[u][v] = mat[u][v] + dct_mat[i][j] * my_cos(PI / n * (u + 1. / 2.) * i) * my_cos(PI / m * (v + 1. / 2.) * j);
                    j = j + 1;
                }
                i = i + 1;
            }
            mat[u][v] = mat[u][v] * 2. / n * 2. / m;
            v = v + 1;
        }
        u = u + 1;
    }
}

int main()
{
    int dim_x = getint(), dim_y = getint();

    memset(test_block, 0, MAX_DIM_X * MAX_DIM_Y);
    memset(test_dct, 0, MAX_DIM_X * MAX_DIM_Y);
    memset(test_idct, 0, MAX_DIM_X * MAX_DIM_Y);

    int i = 0;
    while (i < dim_x) {
        int j = 0;
        while (j < dim_y) {
            test_block[i][j] = getfloat();
            j = j + 1;
        }
        i = i + 1;
    }

    dct(test_dct, test_block, dim_x, dim_y);
    write_mat(test_dct, dim_x, dim_y);

    idct(test_idct, test_dct, dim_x, dim_y);
    write_mat(test_idct, dim_x, dim_y);

    return 0;
}