/*
* imdct.c - naive IDMCT implementation. This can be made O(n log n),
* but as is the rest of the decoder is too slow for this to matter
* we ignore it.
*
* This code is part of the Experimental Haskell MP3 Decoder, version 0.0.1.
* Copyright (c) 2008 Bjorn Edstrom
*
* This software is provided 'as-is', without any express or implied
* warranty. In no event will the authors be held liable for any damages
* arising from the use of this software.
*
* Permission is granted to anyone to use this software for any purpose,
* including commercial applications, and to alter it and redistribute it
* freely, subject to the following restrictions:
*
* 1. The origin of this software must not be misrepresented; you must not
* claim that you wrote the original software. If you use this software
* in a product, an acknowledgment in the product documentation would be
* appreciated but is not required.
*
* 2. Altered source versions must be plainly marked as such, and must not be
* misrepresented as being the original software.
*
* 3. This notice may not be removed or altered from any source
* distribution.
*/
#include
#define PI 3.1415926535897931
void imdct18(double *in, double *out)
{
int k, n;
double sum;
static double lookup[36][18];
static int lookup_init = 0;
if (lookup_init == 0) {
for (n = 0; n < 36; n++)
for (k = 0; k < 18; k++)
lookup[n][k] = cos((PI/18) * (n + 0.5 + 18/2.0) * (k + 0.5));
lookup_init = 1;
}
for (n = 0; n < 36; n++) {
sum = 0.0;
for (k = 0; k < 18; k+=3) {
sum += in[k+0] * lookup[n][k+0];
sum += in[k+1] * lookup[n][k+1];
sum += in[k+2] * lookup[n][k+2];
}
out[n] = sum;
}
}
void imdct(int points, double *in, double *out)
{
int k, n;
double sum;
if (points == 18) {
imdct18(in, out);
return;
}
for (n = 0; n < points*2; n++) {
sum = 0.0;
for (k = 0; k < points; k++)
sum += in[k] * cos((PI/points) * (n + 0.5 + points/2.0) * (k + 0.5));
out[n] = sum;
}
}