<p>The link you posted to a 2012 paper (<a href="https://arxiv.org/pdf/1203.3442.pdf" rel="nofollow noreferrer">https://arxiv.org/pdf/1203.3442.pdf</a>) seems to describe a rather interesting DCT algorithm: It has low computational complexity (32*17 + 255 multiplies for a 16*16 block and 16*5 + 63 for size 8*8) but also a very regular structure, which makes it easy to synthesise a double-sized variant.</p>

<p>When implementing these things, one should mostly just focus on the butterfly graph: Read from left to right to implement forward (Type II) DCT and from right to left to implement inverse (Type III) DCT. Read text and formulas only when needed to interpret any special symbols in the graph.</p>

<p>That being said, I tried to implement the 8-point DCT II sub-module using the graph from the paper. In this case, the 8 outputs, starting from top, should be re-defined as X0, X4, X2, X6, X1, X3, X5, X7. The first five seem to be within a constant factor of reference DCT output, but I don't seem to get the bottom 3 right.</p>

<p>Here's my code that tries to calculate the 8-point transform:</p>

<pre><code>void fDCT2bb2(float* data, const float factor) {
    float a = data[0], b = data[1];
    a -= b;
    b *= factor;
    data[0] = a + b;
    data[1] = b - a;
}

void fDCT2bb4(float* data, const float factor) {
    float a[2] = {data[0], data[1]};
    float b[2] = {data[2], data[3]};
    a[0] -= b[1];
    a[1] -= b[0];
    b[0] *= factor;
    b[1] *= factor;
    data[0] = a[0] + b[0];
    data[1] = a[1] + b[1];
    data[2] = b[0] - a[0];
    data[3] = b[1] - a[1];
}

void fDCT8point(const float* input, float* output) {
    float a[4] = {
        input[0] + input[7],
        input[1] + input[6],
        input[2] + input[5],
        input[3] + input[4]
    };

    float c = a[0];
    a[0] += a[3];
    a[3] -= c;

    c = a[1];
    a[1] += a[2];
    a[2] -= c;

    c = a[0];
    a[0] += a[1];
    a[1] -= c;

    c = a[2];
    a[2] = a[3];
    a[3] = c;

    fDCT2bb2(&amp;a[2], 1.41421356f);

    float b[4] = {
        input[7] - input[0],
        input[6] - input[1],
        input[5] - input[2],
        input[4] - input[3]
    };

    fDCT2bb4(b, 1.41421356f);
    fDCT2bb2(b, 1.84775906f);
    fDCT2bb2(&amp;b[2], -0.76536686f);

    output[0] = a[0];
    output[4] = a[1];
    output[2] = a[2];
    output[6] = a[3];
    output[1] = b[0];
    output[7] = b[1];
    output[5] = b[2];
    output[3] = b[3];
}
</code></pre>

<p>Any simple change to the above seems to make the output worse. I may have misinterpreted how to implement a "building block" with 4 inputs and 4 outputs from the rather terse description, but there shouldn't be too many ways to do things as it's only supposed to have 2 multiplies and 6 adds.</p>

<p><strong>Edit</strong>: I got this fixed by using 5pi/8 instead of 3pi/8 in (=-0.765..) unlike the graph and swapping outputs 3 and 7. Apparently this is an 8-point-only thing, so a 16-point transform should be exactly as in the graph.</p>

<p>Anyway, I have also implemented a similar recursive, regular DCT using this paper: <a href="http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.27.3258&amp;rep=rep1&amp;type=pdf" rel="nofollow noreferrer">http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.27.3258&amp;rep=rep1&amp;type=pdf</a></p>

<p>The butterfly graph is on page 8. After implementing the 8-point variant, it's easy to see how to keep doubling the transform size when needed. The 2-D expansion is not very relevant today, because SIMDifying it is difficult.</p>
How to implement a fast 16 point dct
<p>I'm trying to find a way to perform a fast 16 point dct2 and dct3 transform.</p>

<p>I found some articles like <a href="https://arxiv.org/pdf/1203.3442.pdf" rel="nofollow noreferrer">this one</a> talking about how to do this in mathematical theory, but I'm novice when it comes to reading complex math equations, so honestly I can't understand it.</p>

<p>I searched online for an implementation of a fast 16 point dct, and I found this <a href="http://www.spiral.net/" rel="nofollow noreferrer">code generator</a> which outputs code based on your desired DCT parameters.</p>

<p>I asked it to generate a 16 point dct2 and dct3 with double precision, however the outputs were not mirror images as the inputs when ran through both equations. 
This was my input:</p>

<pre><code>// Before DCT
inputArray[ 0] = 12;
inputArray[ 1] = 12;
inputArray[ 2] = 12;
inputArray[ 3] = 14;
inputArray[ 4] = 8;
inputArray[ 5] = 10;
inputArray[ 6] = 12;
inputArray[ 7] = 12;
inputArray[ 8] = 12;
inputArray[ 9] = 12;
inputArray[10] = 12;
inputArray[11] = 12;
inputArray[12] = 12;
inputArray[13] = 12;
inputArray[14] = 12; 
inputArray[15] = 12;
</code></pre>

<p>And this was my output</p>

<pre><code>// After DCT and IDCT
outputArray[ 0] = 184;
outputArray[ 1] = 194;   
outputArray[ 2] = 178;
outputArray[ 3] = 198;
outputArray[ 4] = 155;  
outputArray[ 5] = 141;
outputArray[ 6] = 164;
outputArray[ 7] = 149;
outputArray[ 8] = 138;
outputArray[ 9] = 121; 
outputArray[10] = 107;
outputArray[11] = 90;
outputArray[12] = 74;  
outputArray[13] = 55;
outputArray[14] = 37;
outputArray[15] = 19;
</code></pre>

<p>I realized the first 5 or so indexes do equal the inputs when divided by 16, however this trend doesn't continue as you go down.</p>

<p>Is this the expected behavior? Or is there something else I need to do the get a proper conversion? </p>

<p>Also I did find an <a href="https://www.nayuki.io/page/fast-discrete-cosine-transform-algorithms" rel="nofollow noreferrer">8 point dct</a> that works well and gives the proper results, is there anyway to expand that into a fast 16 point dct?</p>