You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
550 lines
13 KiB
C
550 lines
13 KiB
C
/*
|
|
* Musepack audio compression
|
|
* Copyright (c) 2005-2009, The Musepack Development Team
|
|
* Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf
|
|
*
|
|
* This library is free software; you can redistribute it and/or
|
|
* modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation; either
|
|
* version 2.1 of the License, or (at your option) any later version.
|
|
*
|
|
* This library is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
* Lesser General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU Lesser General Public
|
|
* License along with this library; if not, write to the Free Software
|
|
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
|
*/
|
|
|
|
#include "libmpcpsy.h"
|
|
#include <mpc/mpcmath.h>
|
|
|
|
/* F U N C T I O N S */
|
|
static void makewt ( const int nw, int* ip, float* w );
|
|
static void makect ( const int nc, int* ip, float* c );
|
|
|
|
# define cftmdl(n,l,a,w) cftmdl_i386 ( n, l, a, w )
|
|
|
|
|
|
/* -------- child routines -------- */
|
|
static mpc_inline void
|
|
bitrv2 ( const int n, int* ip, float* a )
|
|
{
|
|
int j, j1, k, k1, l, m, m2;
|
|
float xr, xi, yr, yi;
|
|
|
|
ip[0] = 0;
|
|
l = n;
|
|
m = 1;
|
|
while ( (m << 3) < l ) {
|
|
l >>= 1;
|
|
for ( j = 0; j < m; j++ ) {
|
|
ip[m + j] = ip[j] + l;
|
|
}
|
|
m <<= 1;
|
|
}
|
|
m2 = 2 * m;
|
|
if ( (m << 3) == l ) {
|
|
for ( k = 0; k < m; k++ ) {
|
|
for ( j = 0; j < k; j++ ) {
|
|
j1 = 2 * j + ip[k];
|
|
k1 = 2 * k + ip[j];
|
|
xr = a[j1];
|
|
xi = a[j1 + 1];
|
|
yr = a[k1];
|
|
yi = a[k1 + 1];
|
|
a[j1] = yr;
|
|
a[j1 + 1] = yi;
|
|
a[k1] = xr;
|
|
a[k1 + 1] = xi;
|
|
j1 += m2;
|
|
k1 += 2 * m2;
|
|
xr = a[j1];
|
|
xi = a[j1 + 1];
|
|
yr = a[k1];
|
|
yi = a[k1 + 1];
|
|
a[j1] = yr;
|
|
a[j1 + 1] = yi;
|
|
a[k1] = xr;
|
|
a[k1 + 1] = xi;
|
|
j1 += m2;
|
|
k1 -= m2;
|
|
xr = a[j1];
|
|
xi = a[j1 + 1];
|
|
yr = a[k1];
|
|
yi = a[k1 + 1];
|
|
a[j1] = yr;
|
|
a[j1 + 1] = yi;
|
|
a[k1] = xr;
|
|
a[k1 + 1] = xi;
|
|
j1 += m2;
|
|
k1 += 2 * m2;
|
|
xr = a[j1];
|
|
xi = a[j1 + 1];
|
|
yr = a[k1];
|
|
yi = a[k1 + 1];
|
|
a[j1] = yr;
|
|
a[j1 + 1] = yi;
|
|
a[k1] = xr;
|
|
a[k1 + 1] = xi;
|
|
}
|
|
j1 = 2 * k + m2 + ip[k];
|
|
k1 = j1 + m2;
|
|
xr = a[j1];
|
|
xi = a[j1 + 1];
|
|
yr = a[k1];
|
|
yi = a[k1 + 1];
|
|
a[j1] = yr;
|
|
a[j1 + 1] = yi;
|
|
a[k1] = xr;
|
|
a[k1 + 1] = xi;
|
|
}
|
|
} else {
|
|
for ( k = 1; k < m; k++ ) {
|
|
for ( j = 0; j < k; j++ ) {
|
|
j1 = 2 * j + ip[k];
|
|
k1 = 2 * k + ip[j];
|
|
xr = a[j1];
|
|
xi = a[j1 + 1];
|
|
yr = a[k1];
|
|
yi = a[k1 + 1];
|
|
a[j1] = yr;
|
|
a[j1 + 1] = yi;
|
|
a[k1] = xr;
|
|
a[k1 + 1] = xi;
|
|
j1 += m2;
|
|
k1 += m2;
|
|
xr = a[j1];
|
|
xi = a[j1 + 1];
|
|
yr = a[k1];
|
|
yi = a[k1 + 1];
|
|
a[j1] = yr;
|
|
a[j1 + 1] = yi;
|
|
a[k1] = xr;
|
|
a[k1 + 1] = xi;
|
|
}
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
|
|
static mpc_inline void
|
|
cft1st ( const int n, float* a, float* w )
|
|
{
|
|
int j, k1;
|
|
float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
|
|
float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
|
|
|
|
x0r = a[ 0] + a[ 2];
|
|
x0i = a[ 1] + a[ 3];
|
|
x1r = a[ 0] - a[ 2];
|
|
x1i = a[ 1] - a[ 3];
|
|
x2r = a[ 4] + a[ 6];
|
|
x2i = a[ 5] + a[ 7];
|
|
x3r = a[ 4] - a[ 6];
|
|
x3i = a[ 5] - a[ 7];
|
|
a[ 0] = x0r + x2r;
|
|
a[ 1] = x0i + x2i;
|
|
a[ 4] = x0r - x2r;
|
|
a[ 5] = x0i - x2i;
|
|
a[ 2] = x1r - x3i;
|
|
a[ 3] = x1i + x3r;
|
|
a[ 6] = x1r + x3i;
|
|
a[ 7] = x1i - x3r;
|
|
wk1r = w[ 2];
|
|
x0r = a[ 8] + a[10];
|
|
x0i = a[ 9] + a[11];
|
|
x1r = a[ 8] - a[10];
|
|
x1i = a[ 9] - a[11];
|
|
x2r = a[12] + a[14];
|
|
x2i = a[13] + a[15];
|
|
x3r = a[12] - a[14];
|
|
x3i = a[13] - a[15];
|
|
a[ 8] = x0r + x2r;
|
|
a[ 9] = x0i + x2i;
|
|
a[12] = x2i - x0i;
|
|
a[13] = x0r - x2r;
|
|
x0r = x1r - x3i;
|
|
x0i = x1i + x3r;
|
|
a[10] = wk1r * (x0r - x0i);
|
|
a[11] = wk1r * (x0r + x0i);
|
|
x0r = x3i + x1r;
|
|
x0i = x3r - x1i;
|
|
a[14] = wk1r * (x0i - x0r);
|
|
a[15] = wk1r * (x0i + x0r);
|
|
|
|
k1 = 0;
|
|
j = 16;
|
|
do {
|
|
k1 += 2;
|
|
wk2r = w[k1];
|
|
wk2i = w[k1 + 1];
|
|
wk1r = w[2*k1];
|
|
wk1i = w[2*k1 + 1];
|
|
wk3r = wk1r - 2 * wk2i * wk1i;
|
|
wk3i = 2 * wk2i * wk1r - wk1i;
|
|
x0r = a[j] + a[j + 2];
|
|
x0i = a[j + 1] + a[j + 3];
|
|
x1r = a[j] - a[j + 2];
|
|
x1i = a[j + 1] - a[j + 3];
|
|
x2r = a[j + 4] + a[j + 6];
|
|
x2i = a[j + 5] + a[j + 7];
|
|
x3r = a[j + 4] - a[j + 6];
|
|
x3i = a[j + 5] - a[j + 7];
|
|
a[j] = x0r + x2r;
|
|
a[j + 1] = x0i + x2i;
|
|
x0r -= x2r;
|
|
x0i -= x2i;
|
|
a[j + 4] = wk2r * x0r - wk2i * x0i;
|
|
a[j + 5] = wk2r * x0i + wk2i * x0r;
|
|
x0r = x1r - x3i;
|
|
x0i = x1i + x3r;
|
|
a[j + 2] = wk1r * x0r - wk1i * x0i;
|
|
a[j + 3] = wk1r * x0i + wk1i * x0r;
|
|
x0r = x1r + x3i;
|
|
x0i = x1i - x3r;
|
|
a[j + 6] = wk3r * x0r - wk3i * x0i;
|
|
a[j + 7] = wk3r * x0i + wk3i * x0r;
|
|
wk1r = w[2*k1 + 2];
|
|
wk1i = w[2*k1 + 3];
|
|
wk3r = wk1r - 2 * wk2r * wk1i;
|
|
wk3i = 2 * wk2r * wk1r - wk1i;
|
|
x0r = a[j + 8] + a[j + 10];
|
|
x0i = a[j + 9] + a[j + 11];
|
|
x1r = a[j + 8] - a[j + 10];
|
|
x1i = a[j + 9] - a[j + 11];
|
|
x2r = a[j + 12] + a[j + 14];
|
|
x2i = a[j + 13] + a[j + 15];
|
|
x3r = a[j + 12] - a[j + 14];
|
|
x3i = a[j + 13] - a[j + 15];
|
|
a[j + 8] = x0r + x2r;
|
|
a[j + 9] = x0i + x2i;
|
|
x0r -= x2r;
|
|
x0i -= x2i;
|
|
a[j + 12] = -wk2i * x0r - wk2r * x0i;
|
|
a[j + 13] = -wk2i * x0i + wk2r * x0r;
|
|
x0r = x1r - x3i;
|
|
x0i = x1i + x3r;
|
|
a[j + 10] = wk1r * x0r - wk1i * x0i;
|
|
a[j + 11] = wk1r * x0i + wk1i * x0r;
|
|
x0r = x1r + x3i;
|
|
x0i = x1i - x3r;
|
|
a[j + 14] = wk3r * x0r - wk3i * x0i;
|
|
a[j + 15] = wk3r * x0i + wk3i * x0r;
|
|
} while ( j += 16, j < n );
|
|
return;
|
|
}
|
|
|
|
static mpc_inline void
|
|
cftmdl_i386 ( const int n, const int l, float* a, float* w )
|
|
{
|
|
int j, j1, j2, j3, k, k1, m, m2;
|
|
float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
|
|
float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
|
|
|
|
m = l << 2;
|
|
|
|
for ( j = 0; j < l; j += 2 ) {
|
|
j1 = j + l;
|
|
j2 = j1 + l;
|
|
j3 = j2 + l;
|
|
x0r = a[j] + a[j1];
|
|
x0i = a[j + 1] + a[j1 + 1];
|
|
x1r = a[j] - a[j1];
|
|
x1i = a[j + 1] - a[j1 + 1];
|
|
x2r = a[j2] + a[j3];
|
|
x2i = a[j2 + 1] + a[j3 + 1];
|
|
x3r = a[j2] - a[j3];
|
|
x3i = a[j2 + 1] - a[j3 + 1];
|
|
a[j] = x0r + x2r;
|
|
a[j + 1] = x0i + x2i;
|
|
a[j2] = x0r - x2r;
|
|
a[j2 + 1] = x0i - x2i;
|
|
a[j1] = x1r - x3i;
|
|
a[j1 + 1] = x1i + x3r;
|
|
a[j3] = x1r + x3i;
|
|
a[j3 + 1] = x1i - x3r;
|
|
}
|
|
|
|
wk1r = w[2];
|
|
for ( j = m; j < l + m; j += 2 ) {
|
|
j1 = j + l;
|
|
j2 = j1 + l;
|
|
j3 = j2 + l;
|
|
x0r = a[j] + a[j1];
|
|
x0i = a[j + 1] + a[j1 + 1];
|
|
x1r = a[j] - a[j1];
|
|
x1i = a[j + 1] - a[j1 + 1];
|
|
x2r = a[j2] + a[j3];
|
|
x2i = a[j2 + 1] + a[j3 + 1];
|
|
x3r = a[j2] - a[j3];
|
|
x3i = a[j2 + 1] - a[j3 + 1];
|
|
a[j] = x0r + x2r;
|
|
a[j + 1] = x0i + x2i;
|
|
a[j2] = x2i - x0i;
|
|
a[j2 + 1] = x0r - x2r;
|
|
x0r = x1r - x3i;
|
|
x0i = x1i + x3r;
|
|
a[j1] = wk1r * (x0r - x0i);
|
|
a[j1 + 1] = wk1r * (x0r + x0i);
|
|
x0r = x3i + x1r;
|
|
x0i = x3r - x1i;
|
|
a[j3] = wk1r * (x0i - x0r);
|
|
a[j3 + 1] = wk1r * (x0i + x0r);
|
|
}
|
|
|
|
k1 = 0;
|
|
m2 = 2 * m;
|
|
for ( k = m2; k < n; k += m2 ) {
|
|
k1 += 2;
|
|
wk2r = w[k1];
|
|
wk2i = w[k1 + 1];
|
|
wk1r = w[2*k1];
|
|
wk1i = w[2*k1 + 1];
|
|
wk3r = wk1r - 2 * wk2i * wk1i;
|
|
wk3i = 2 * wk2i * wk1r - wk1i;
|
|
j = k;
|
|
do {
|
|
j1 = j + l;
|
|
j2 = j1 + l;
|
|
j3 = j2 + l;
|
|
x0r = a[j] + a[j1];
|
|
x0i = a[j + 1] + a[j1 + 1];
|
|
x1r = a[j] - a[j1];
|
|
x1i = a[j + 1] - a[j1 + 1];
|
|
x2r = a[j2] + a[j3];
|
|
x2i = a[j2 + 1] + a[j3 + 1];
|
|
x3r = a[j2] - a[j3];
|
|
x3i = a[j2 + 1] - a[j3 + 1];
|
|
a[j] = x0r + x2r;
|
|
a[j + 1] = x0i + x2i;
|
|
x0r -= x2r;
|
|
x0i -= x2i;
|
|
a[j2] = wk2r * x0r - wk2i * x0i;
|
|
a[j2 + 1] = wk2r * x0i + wk2i * x0r;
|
|
x0r = x1r - x3i;
|
|
x0i = x1i + x3r;
|
|
a[j1] = wk1r * x0r - wk1i * x0i;
|
|
a[j1 + 1] = wk1r * x0i + wk1i * x0r;
|
|
x0r = x1r + x3i;
|
|
x0i = x1i - x3r;
|
|
a[j3] = wk3r * x0r - wk3i * x0i;
|
|
a[j3 + 1] = wk3r * x0i + wk3i * x0r;
|
|
} while ( j += 2, j < l + k );
|
|
|
|
wk1r = w[2*k1 + 2];
|
|
wk1i = w[2*k1 + 3];
|
|
wk3r = wk1r - 2 * wk2r * wk1i;
|
|
wk3i = 2 * wk2r * wk1r - wk1i;
|
|
j = k + m;
|
|
do {
|
|
j1 = j + l;
|
|
j2 = j1 + l;
|
|
j3 = j2 + l;
|
|
x0r = a[j] + a[j1];
|
|
x0i = a[j + 1] + a[j1 + 1];
|
|
x1r = a[j] - a[j1];
|
|
x1i = a[j + 1] - a[j1 + 1];
|
|
x2r = a[j2] + a[j3];
|
|
x2i = a[j2 + 1] + a[j3 + 1];
|
|
x3r = a[j2] - a[j3];
|
|
x3i = a[j2 + 1] - a[j3 + 1];
|
|
a[j] = x0r + x2r;
|
|
a[j + 1] = x0i + x2i;
|
|
x0r -= x2r;
|
|
x0i -= x2i;
|
|
a[j2] = -wk2i * x0r - wk2r * x0i;
|
|
a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
|
|
x0r = x1r - x3i;
|
|
x0i = x1i + x3r;
|
|
a[j1] = wk1r * x0r - wk1i * x0i;
|
|
a[j1 + 1] = wk1r * x0i + wk1i * x0r;
|
|
x0r = x1r + x3i;
|
|
x0i = x1i - x3r;
|
|
a[j3] = wk3r * x0r - wk3i * x0i;
|
|
a[j3 + 1] = wk3r * x0i + wk3i * x0r;
|
|
} while ( j += 2, j < l+k+m );
|
|
}
|
|
return;
|
|
}
|
|
|
|
static mpc_inline void
|
|
cftfsub ( const int n, float* a, float* w )
|
|
{
|
|
int j, j1, j2, j3, l;
|
|
float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
|
|
|
|
l = 2;
|
|
if ( n > 8 ) {
|
|
cft1st ( n, a, w );
|
|
l = 8;
|
|
while ( (l << 2) < n ) {
|
|
cftmdl ( n, l, a, w );
|
|
l <<= 2;
|
|
}
|
|
}
|
|
if ( (l << 2) == n ) {
|
|
j = 0;
|
|
do {
|
|
j1 = j + l;
|
|
j2 = j1 + l;
|
|
j3 = j2 + l;
|
|
x0r = a[j] + a[j1];
|
|
x0i = a[j + 1] + a[j1 + 1];
|
|
x1r = a[j] - a[j1];
|
|
x1i = a[j + 1] - a[j1 + 1];
|
|
x2r = a[j2] + a[j3];
|
|
x2i = a[j2 + 1] + a[j3 + 1];
|
|
x3r = a[j2] - a[j3];
|
|
x3i = a[j2 + 1] - a[j3 + 1];
|
|
a[j] = x0r + x2r;
|
|
a[j + 1] = x0i + x2i;
|
|
a[j2] = x0r - x2r;
|
|
a[j2 + 1] = x0i - x2i;
|
|
a[j1] = x1r - x3i;
|
|
a[j1 + 1] = x1i + x3r;
|
|
a[j3] = x1r + x3i;
|
|
a[j3 + 1] = x1i - x3r;
|
|
} while ( j += 2, j < l );
|
|
} else {
|
|
j = 0;
|
|
do {
|
|
j1 = j + l;
|
|
x0r = a[j] - a[j1];
|
|
x0i = a[j + 1] - a[j1 + 1];
|
|
a[j] += a[j1];
|
|
a[j + 1] += a[j1 + 1];
|
|
a[j1] = x0r;
|
|
a[j1 + 1] = x0i;
|
|
} while ( j += 2, j < l );
|
|
}
|
|
return;
|
|
}
|
|
|
|
|
|
static mpc_inline void
|
|
rftfsub ( const int n, float* a, int nc, float* c )
|
|
{
|
|
int j, k, kk, ks, m;
|
|
float wkr, wki, xr, xi, yr, yi;
|
|
|
|
m = n >> 1;
|
|
ks = 2 * nc / m;
|
|
kk = ks;
|
|
j = 2;
|
|
k = n;
|
|
do {
|
|
k -= 2;
|
|
nc -= ks;
|
|
wkr = 0.5f - c[nc];
|
|
wki = c[kk];
|
|
xr = a[j] - a[k];
|
|
xi = a[j + 1] + a[k + 1];
|
|
yr = wkr * xr - wki * xi;
|
|
yi = wkr * xi + wki * xr;
|
|
a[j] -= yr;
|
|
a[j + 1] -= yi;
|
|
a[k] += yr;
|
|
a[k + 1] -= yi;
|
|
kk += ks;
|
|
} while ( j += 2, j < m );
|
|
return;
|
|
}
|
|
|
|
// generates lookup-tables
|
|
void
|
|
Generate_FFT_Tables ( const int n, int* ip, float* w )
|
|
{
|
|
int nw;
|
|
int nc;
|
|
|
|
nw = n >> 2;
|
|
makewt ( nw, ip, w );
|
|
|
|
nc = n >> 2;
|
|
makect ( nc, ip, w + nw );
|
|
}
|
|
|
|
|
|
// patched to only-forward
|
|
void
|
|
rdft ( const int n, float* a, int* ip, float* w )
|
|
{
|
|
float xi;
|
|
|
|
if ( n > 4) {
|
|
bitrv2 ( n, ip + 2, a );
|
|
cftfsub ( n, a, w );
|
|
rftfsub ( n, a, ip[1], w + ip[0] );
|
|
}
|
|
else if ( n == 4 ) {
|
|
cftfsub ( n, a, w );
|
|
}
|
|
xi = a[0] - a[1];
|
|
a[0] += a[1];
|
|
a[1] = xi;
|
|
return;
|
|
}
|
|
|
|
|
|
/* -------- initializing routines -------- */
|
|
static void
|
|
makewt ( const int nw, int* ip, float* w )
|
|
{
|
|
int j;
|
|
int nwh;
|
|
float x;
|
|
float y;
|
|
double delta;
|
|
|
|
ip[0] = nw;
|
|
ip[1] = 1;
|
|
if ( nw > 2 ) {
|
|
nwh = nw >> 1;
|
|
delta = (M_PI/4) / nwh;
|
|
w[0] = 1.;
|
|
w[1] = 0.;
|
|
w[nwh] = COS (delta * nwh);
|
|
w[nwh + 1] = w[nwh];
|
|
if ( nwh > 2 ) {
|
|
for ( j = 2; j < nwh; j += 2 ) {
|
|
x = COS (delta * j);
|
|
y = SIN (delta * j);
|
|
w[j] = x;
|
|
w[j + 1] = y;
|
|
w[nw - j] = y;
|
|
w[nw - j + 1] = x;
|
|
}
|
|
bitrv2 ( nw, ip + 2, w );
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
|
|
static void
|
|
makect ( const int nc, int* ip, float* c )
|
|
{
|
|
int j;
|
|
int nch;
|
|
double delta;
|
|
|
|
ip[1] = nc;
|
|
if ( nc > 1 ) {
|
|
nch = nc >> 1;
|
|
delta = (M_PI/4) / nch;
|
|
c[0] = COS (delta * nch);
|
|
c[nch] = 0.5f * c[0];
|
|
for ( j = 1; j < nch; j++ ) {
|
|
c[j] = 0.5f * COS (delta * j);
|
|
c[nc - j] = 0.5f * SIN (delta * j);
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
/* end of fft4g.c */
|