Faster AAN algorithm for calculating Discrete Cosine Transform

Question

I was implementing the AAN algorithm for computing the 2D Discrete Cosine Transform on an 8*8 array. I perform both the forward and inverse AAN algorithms.

My code is extremely fast compared to the naive DCT algorithm but I feel that my implementation could be faster. How can I improve the performance of my transforms?

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifndef M_PI
#    define M_PI 3.14159265358979323846
#endif
/*
File: dct.c
Summary: 2D AAN implementation of Discrete Cosine transform
     Code for Forward and Inverse DCT using AAN method
Sources: Forward DCT inspired by // https://unix4lyfe.org/dct-1d/
     Inverse DCT inspired by (Everything you need to know about jpeg Youtube Playlist) https://www.youtube.com/watch?v=sb8CQ9knDgI&list=PLpsTn9TA_Q8VMDyOPrDKmSJYt1DLgDZU4

Input : float array of 64 elements representing 8 * 8 block
*Note I use a 1-dimensional array to represent 2d data   
*/


//Compute forward dct
void ForwardDCTComponent(float *component)
{
/*Forward DCT constants*/
const float a1 = 0.707;
const float a2 = 0.541;
const float a3 = 0.707;
const float a4 = 1.307;
const float a5 = 0.383;

const float s0 = 0.353553;
const float s1 = 0.254898;
const float s2 = 0.270598;
const float s3 = 0.300672;
const float s4 = s0;
const float s5 = 0.449988;
const float s6 = 0.653281;
const float s7 = 1.281458;
/*Forward DCT constants*/


    for(int i = 0; i < 8; i++)
    {
        const float b0 = component[0*8 + i] + component[7*8 + i];
        const float b1 = component[1*8 + i] + component[6*8 + i];
        const float b2 = component[2*8 + i] + component[5*8 + i];
        const float b3 = component[3*8 + i] + component[4*8 + i];
        const float b4 =-component[4*8 + i] + component[3*8 + i];
        const float b5 =-component[5*8 + i] + component[2*8 + i];
        const float b6 =-component[6*8 + i] + component[1*8 + i];
        const float b7 =-component[7*8 + i] + component[0*8 + i];
        
        const float c0 = b0 + b3;
        const float c1 = b1 + b2;
        const float c2 =-b2 + b1;
        const float c3 =-b3 + b0;
        const float c4 =-b4 - b5;
        const float c5 = b5 + b6;
        const float c6 = b6 + b7;
        const float c7 = b7;
        
        const float d0 = c0 + c1;
        const float d1 =-c1 + c0;
        const float d2 = c2 + c3;
        const float d3 = c3;
        const float d4 = c4;
        const float d5 = c5;
        const float d6 = c6;
        const float d7 = c7;
        const float d8 = (d4+d6) * a5;
        
        const float e0 = d0;
        const float e1 = d1;
        const float e2 = d2 * a1;
        const float e3 = d3;
        const float e4 = -d4 * a2 - d8;
        const float e5 = d5 * a3;
        const float e6 = d6 * a4 - d8;
        const float e7 = d7;
        
        const float f0 = e0;
        const float f1 = e1;
        const float f2 = e2 + e3;
        const float f3 = e3 - e2;
        const float f4 = e4;
        const float f5 = e5 + e7;
        const float f6 = e6;
        const float f7 = e7 - e5;
        
        const float g0 = f0;
        const float g1 = f1;
        const float g2 = f2;
        const float g3 = f3;
        const float g4 = f4 + f7;
        const float g5 = f5 + f6;
        const float g6 = -f6 + f5;
        const float g7 = f7 - f4;   
        
        component[0*8 + i] = g0 * s0;
        component[4*8 + i] = g1 * s4;
        component[2*8 + i] = g2 * s2;
        component[6*8 + i] = g3 * s6;
        component[5*8 + i] = g4 * s5;
        component[1*8 + i] = g5 * s1;
        component[7*8 + i] = g6 * s7;
        component[3*8 + i] = g7 * s3;       
    }
    for(int i = 0; i < 8; i++)
    {
        const float b0 = component[i*8 + 0]  + component[i*8 + 7];
        const float b1 = component[i*8 + 1]  + component[i*8 + 6];
        const float b2 = component[i*8 + 2]  + component[i*8 + 5];
        const float b3 = component[i*8 + 3]  + component[i*8 + 4];
        const float b4 =-component[i*8 + 4]  + component[i*8 + 3];
        const float b5 =-component[i*8 + 5]  + component[i*8 + 2];
        const float b6 =-component[i*8 + 6]  + component[i*8 + 1];
        const float b7 =-component[i*8 + 7]  + component[i*8 + 0] ;
        
        const float c0 = b0 + b3;
        const float c1 = b1 + b2;
        const float c2 =-b2 + b1;
        const float c3 =-b3 + b0;
        const float c4 =-b4 - b5;
        const float c5 = b5 + b6;
        const float c6 = b6 + b7;
        const float c7 = b7;
        
        const float d0 = c0 + c1;
        const float d1 =-c1 + c0;
        const float d2 = c2 + c3;
        const float d3 = c3;
        const float d4 = c4;
        const float d5 = c5;
        const float d6 = c6;
        const float d7 = c7;
        const float d8 = (d4+d6) * a5;
        
        const float e0 = d0;
        const float e1 = d1;
        const float e2 = d2 * a1;
        const float e3 = d3;
        const float e4 = -d4 * a2 - d8;
        const float e5 = d5 * a3;
        const float e6 = d6 * a4 - d8;
        const float e7 = d7;
        
        const float f0 = e0;
        const float f1 = e1;
        const float f2 = e2 + e3;
        const float f3 = e3 - e2;
        const float f4 = e4;
        const float f5 = e5 + e7;
        const float f6 = e6;
        const float f7 = e7 - e5;
        
        const float g0 = f0;
        const float g1 = f1;
        const float g2 = f2;
        const float g3 = f3;
        const float g4 = f4 + f7;
        const float g5 = f5 + f6;
        const float g6 = -f6 + f5;
        const float g7 = f7 - f4;   
        
        component[i*8 + 0] = g0 * s0;
        component[i*8 + 4] = g1 * s4;
        component[i*8 + 2] = g2 * s2;
        component[i*8 + 6] = g3 * s6;
        component[i*8 + 5] = g4 * s5;
        component[i*8 + 1] = g5 * s1;
        component[i*8 + 7] = g6 * s7;
        component[i*8 + 3] = g7 * s3;   
    }
}

void InverseDCTComponent(float *component)
{
const float m0 = 2.0 * cos(1.0/16.0 * 2.0 * M_PI);
const float m1 = 2.0 * cos(2.0/16.0 * 2.0 * M_PI);
const float m3 = 2.0 * cos(2.0/16.0 * 2.0 * M_PI);
const float m5 = 2.0 * cos(3.0/16.0 * 2.0 * M_PI);
const float m2 = m0-m5;
const float m4 = m0+m5;

const float s0 = cos(0.0/16.0 *M_PI)/sqrt(8);
const float s1 = cos(1.0/16.0 *M_PI)/2.0;
const float s2 = cos(2.0/16.0 *M_PI)/2.0;
const float s3 = cos(3.0/16.0 *M_PI)/2.0;
const float s4 = cos(4.0/16.0 *M_PI)/2.0;
const float s5 = cos(5.0/16.0 *M_PI)/2.0;
const float s6 = cos(6.0/16.0 *M_PI)/2.0;
const float s7 = cos(7.0/16.0 *M_PI)/2.0;

    for(int i = 0; i < 8; i++)
    {
        const float g0 = component[0*8 + i] * s0;
        const float g1 = component[4*8 + i] * s4;
        const float g2 = component[2*8 + i] * s2;
        const float g3 = component[6*8 + i] * s6;
        const float g4 = component[5*8 + i] * s5;
        const float g5 = component[1*8 + i] * s1;
        const float g6 = component[7*8 + i] * s7;
        const float g7 = component[3*8 + i] * s3;
        
        const float f0 = g0;
        const float f1 = g1;
        const float f2 = g2;
        const float f3 = g3;
        const float f4 = g4 - g7;
        const float f5 = g5 + g6;
        const float f6 = g5 - g6;
        const float f7 = g4 + g7;
        
        const float e0 = f0;
        const float e1 = f1;
        const float e2 = f2 - f3;
        const float e3 = f2 + f3;
        const float e4 = f4;
        const float e5 = f5 - f7;
        const float e6 = f6;
        const float e7 = f5 + f7;
        const float e8 = f4 + f6;
        
        const float d0 = e0;
        const float d1 = e1;
        const float d2 = e2 * m1;
        const float d3 = e3;
        const float d4 = e4 * m2;
        const float d5 = e5 * m3;
        const float d6 = e6 * m4;
        const float d7 = e7;
        const float d8 = e8 * m5;
        
        const float c0 = d0 + d1;
        const float c1 = d0 - d1;
        const float c2 = d2 - d3;
        const float c3 = d3;
        const float c4 = d4 + d8;
        const float c5 = d5 + d7;
        const float c6 = d6 - d8;
        const float c7 = d7;
        const float c8 = c5 - c6;
        
        const float b0 = c0 + c3;
        const float b1 = c1 + c2;
        const float b2 = c1 - c2;
        const float b3 = c0 - c3;
        const float b4 = c4 - c8;
        const float b5 = c8;
        const float b6 = c6 - c7;
        const float b7 = c7;
        
        component[0 * 8 + i] = b0 + b7;
        component[1 * 8 + i] = b1 + b6;
        component[2 * 8 + i] = b2 + b5;
        component[3 * 8 + i] = b3 + b4;
        component[4 * 8 + i] = b3 - b4;
        component[5 * 8 + i] = b2 - b5;
        component[6 * 8 + i] = b1 - b6;
        component[7 * 8 + i] = b0 - b7; 
    }
    
    for(int i = 0 ; i < 8; i++)
    {
        const float g0 = component[i*8 + 0] * s0;
        const float g1 = component[i*8 + 4] * s4;
        const float g2 = component[i*8 + 2] * s2;
        const float g3 = component[i*8 + 6] * s6;
        const float g4 = component[i*8 + 5] * s5;
        const float g5 = component[i*8 + 1] * s1;
        const float g6 = component[i*8 + 7] * s7;
        const float g7 = component[i*8 + 3] * s3;
        
        const float f0 = g0;
        const float f1 = g1;
        const float f2 = g2;
        const float f3 = g3;
        const float f4 = g4 - g7;
        const float f5 = g5 + g6;
        const float f6 = g5 - g6;
        const float f7 = g4 + g7;
        
        const float e0 = f0;
        const float e1 = f1;
        const float e2 = f2 - f3;
        const float e3 = f2 + f3;
        const float e4 = f4;
        const float e5 = f5 - f7;
        const float e6 = f6;
        const float e7 = f5 + f7;
        const float e8 = f4 + f6;
        
        const float d0 = e0;
        const float d1 = e1;
        const float d2 = e2 * m1;
        const float d3 = e3;
        const float d4 = e4 * m2;
        const float d5 = e5 * m3;
        const float d6 = e6 * m4;
        const float d7 = e7;
        const float d8 = e8 * m5;
        
        const float c0 = d0 + d1;
        const float c1 = d0 - d1;
        const float c2 = d2 - d3;
        const float c3 = d3;
        const float c4 = d4 + d8;
        const float c5 = d5 + d7;
        const float c6 = d6 - d8;
        const float c7 = d7;
        const float c8 = c5 - c6;
        
        const float b0 = c0 + c3;
        const float b1 = c1 + c2;
        const float b2 = c1 - c2;
        const float b3 = c0 - c3;
        const float b4 = c4 - c8;
        const float b5 = c8;
        const float b6 = c6 - c7;
        const float b7 = c7;
        
        component[i * 8 + 0] = b0 + b7;
        component[i * 8 + 1] = b1 + b6;
        component[i * 8 + 2] = b2 + b5;
        component[i * 8 + 3] = b3 + b4;
        component[i * 8 + 4] = b3 - b4;
        component[i * 8 + 5] = b2 - b5;
        component[i * 8 + 6] = b1 - b6;
        component[i * 8 + 7] = b0 - b7;
    }
}

void Print(float *array, int length)
{
    for(int i = 0; i < 64; i++)
    {
        if(i > 0 && i % 8 == 0)
        {
            putchar('\n');
        }
        printf("%.1f ", array[i]);  
    }
    putchar('\n');
}




int main()
{
    //Create Sample Array 
    float *component = malloc(64 * sizeof(int));
    for(int i = 0; i < 64; i++)
    {
        component[i] = 255;
    }
    //Print Sample Array
    Print(component,64);
    
    //Perform forward DCT
    ForwardDCTComponent(component);
    putchar('\n');
    //Print results of forward dct
    Print(component,64);
    
    //Perform inverse DCT
    InverseDCTComponent(component);
    putchar('\n');
    //Print results of inverse dct
    Print(component,64);
    free(component);
    
    
    return 0;
}

What sort of things are in scope for improving the performance? Maybe SIMD? Fixed-point approximations? — user555045, Commented Jul 30, 2021 at 7:08

chux · Accepted Answer · 2021-07-30 05:25:39Z

Performance

How can I improve the performance of my transforms?

So far, only some small ideas.

float vs. double math

Do not call run-time double functions and constants for float math. A lot of time can be spent up and down converting.

// const float s1 = cos(1.0/16.0 * M_PI)/2.0;
define M_PIf 3.14159265f
const float s1 = cosf(1.0f/16.0f * M_PIf)/2.0f;

On the other hand, even better to pre-calculate with high precision and then ...

Use constants

Some compilers will form compile time constants values involving constant equations with well known functions at compile time, others calculate at run-time. Simple coding the value is fastest. Some issues apply concerng how many digits to code. I reccomend at least 9 for float and OK to have many more. Pedantic code would use 17+ in case float has same precision as double. Your call.

// const float s0 = cos(0.0/16.0 *M_PI)/sqrt(8);
const float s0 = 0.353553391f // cos(0.0/16.0 *M_PI)/sqrt(8);

static

Depending on the compiler, a const object may get re-loaded in a function.

static will certainly prevent that.

// const float s0 = 0.353553391f // cos(0.0/16.0 *M_PI)/sqrt(8);
static const float s0 = 0.353553391f // cos(0.0/16.0 *M_PI)/sqrt(8);

Micro optimization

Might save a tick or two

 //for(int i = 0 ; i < 8; i++)
 //  ...  i*8 ...
 //  ...  i*8 ...
 for(int i = 0 ; i < 64; i+=8)
   ...  i ...
   ...  i ...

Save coding time

This is a saving of your valuable coding time.

void Print(float *array, int length) does not use length. With a well enabled compiler, that would have generated a warning. Yet since it is there, it implies all warnings were not enabled.

Save time, enable all warnings to guide you to better code.

Correctness

Coding float vs. double constants

// const float s0 = 0.353553;
const float s0 = 0.353553f;

See When does appending an 'f' change the value of a floating constant when assigned to a float?

Suspicious values

For a1, a2, ..., why only 3 significant digits constants? Perhaps float a1 = 0.707106781f; is better?

Assessing

To judge the correctness of a float, avoid using "%.1f" as that takes the floating out of floating point. Wee values are all 0.0 and large values have excessive digits.

// printf("%.1f ", array[i]);
printf("%.*g ", FLT_DECIMAL_DIG, array[i]);

Latent Bug

float and int are not specified to have the same size.

Avoid sizing errors. Size to the referenced object, not type. Easier to code right, review and maintain.

// float *component = malloc(64 * sizeof(int)); // int ???
float *component = malloc(sizeof *component * 64);

Stack Exchange Network

Faster AAN algorithm for calculating Discrete Cosine Transform

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
performance
c
image
or ask your own question.

Hot Network Questions

Faster AAN algorithm for calculating Discrete Cosine Transform

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged performancecimage or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
performance
c
image
or ask your own question.