TDAC-MDCTs of a sweep

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

using u16 = unsigned __int16;

#define MAX     (3 << 15)
#define WINSIZE 240

static double  X   [MAX+2048];
static double  Y   [MAX+2048];
static double* Z   [MAX+2048];
static u16     Zlen[MAX+2048];

constexpr double Pi = 3.1415926535897932384626433832795;

static double
w0 (int N, int k /*0...2N-1*/)
{
    double const q = (k - N + 0.5) / N;
    return ::fabs(q) < 0.5 ? 1.0 : 0.0;
}

static double
w1 (int N, int k /*0...2N-1*/)
{
    double const q = (k - N + 0.5) / N;
    double const c = ::cos (Pi / 2 * q);
    return c;
}

static double
w2 (int N, int k /*0...2N-1*/)
{
    double const q = (k - N + 0.5) / N;
    double const c = ::cos (Pi / 2 * q);
    return sin (Pi/2*c*c);
}

#define W(N,k)  w2(N,k)

static double
cs (int n, int N, int k /*0...2N-1*/)
{
    return ::cos (Pi/4 / N * (2*n+1+N) * (2*k+1));
}

static double 
MDCT (int _N0, int N, int k /*0...2N-1*/)
{
    double S = 0.0;
    for (int n = 0; n < 2*N; n++)
        S += X[_N0 + n] * W(N, n) * cs(n, N, k);
    return -S / N * sqrt(2);
}

static double 
iMDCT (int _N0, int N, int n /*0...N-1*/)
{
    double S = 0.0;
    for (int k = 0; k < N; k++)
        S += Z[_N0][k] * cs(n, N, k);
    return -S * W(N, n);
}

static void 
Process ()
{
    for (int N0 = 0; N0 <= MAX; N0 += WINSIZE)
    {
        Z[N0] = (double*)::calloc (sizeof(double), Zlen[N0] = WINSIZE);
        for (int k = 0; k < WINSIZE; k++)
            Z[N0][k] = ::MDCT (N0, WINSIZE, k);
    }
    
    for (int N0 = 0; N0 < MAX; N0 += WINSIZE)
    {
        for (int n = 0; n < 2*WINSIZE; n++)
            Y[N0+n] += iMDCT (N0, WINSIZE, n) * sqrt(2);
    }
}

static void 
Generate ()
{
    for (int i = 0; i <= MAX; i++)
        X[i] = 32000.0 * ::cos (Pi*i*i*i / MAX / MAX / 3);
}

static void 
Print1 (FILE* const fp)
{
    for (int i = 0; i < MAX; i++)
    {
        ::fprintf (fp, "%7u\t%15.6f\t%15.6f", i, X[i], Y[i]);
        for (int j = 0; j < Zlen[i]; j++)
            ::fprintf (fp, "\t%15.6f", Z[i][j]);
        ::fprintf (fp, "\n");
    }
}

static void 
Print1 ()
{
    FILE* const fp = ::fopen ("mdct.csv", "wb");
    if (fp)
    {
        ::Print1 (fp);
        ::fclose (fp);
    }
}

static void 
Print2 (FILE* const fp)
{
    for (int j = 0; j < 1024; j++)
    {
        for (int i = 0; i < MAX; i++)
            if (Zlen[i])
                if (Zlen[i] > j)
                    ::fprintf (fp, "\t%15.6f", Z[i][j]);
                else
                    ::fprintf (fp, "\t");
        ::fprintf (fp, "\n");
    }
}

static void 
Print2 ()
{
    FILE* const fp = ::fopen ("mdct_freq.csv", "wb");
    if (fp)
    {
        ::Print2 (fp);
        ::fclose (fp);
    }
}

static void
Window ()
{
    for (int k = 0; k < WINSIZE*2; k++)
        ::printf ("%4u\t%12.6f\t%12.6f\t%12.6f\n", k, w0(WINSIZE,k), w1(WINSIZE,k), w2(WINSIZE,k));
    ::printf ("\n");
}

int main ()
{
    ::Window();
    ::fprintf (stderr, "Generate %u\n", MAX);
    ::Generate();
    ::fprintf (stderr, "Process\n");
    ::Process();
    ::fprintf (stderr, "Print1\n");
    ::Print1();
    ::fprintf (stderr, "Print2\n");
    ::Print2();
    ::fprintf (stderr, "Ready\n");
    return 0;
}