diff options
Diffstat (limited to 'ac3dec/imdct_c.c')
| -rw-r--r-- | ac3dec/imdct_c.c | 218 |
1 files changed, 218 insertions, 0 deletions
diff --git a/ac3dec/imdct_c.c b/ac3dec/imdct_c.c new file mode 100644 index 00000000..1f2bfe84 --- /dev/null +++ b/ac3dec/imdct_c.c @@ -0,0 +1,218 @@ +/* + * imdct.c + * + * Copyright (C) Aaron Holtzman - May 1999 + * + * This file is part of ac3dec, a free Dolby AC-3 stream decoder. + * + * ac3dec is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * ac3dec 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Make; see the file COPYING. If not, write to + * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + * + * + */ + +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include "ac3.h" +#include "ac3_internal.h" + +#include "downmix.h" +#include "imdct_c.h" +#include "srfft.h" + + +#define N 512 + +extern void (*imdct_do_512) (float data[],float delay[]); +extern void (*imdct_do_512_nol) (float data[], float delay[]); +extern void (*fft_64p) (complex_t *); + +extern const int pm128[]; +extern float window[]; +extern complex_t buf[128]; + +extern void fft_64p_c (complex_t *); +extern void fft_128p_c (complex_t *); + +static void imdct_do_512_c (float data[],float delay[]); +static void imdct_do_512_nol_c (float data[], float delay[]); + +/* Twiddle factors for IMDCT */ +static float xcos1[128] __attribute__((aligned(16))); +static float xsin1[128] __attribute__((aligned(16))); + + +int imdct_init_c (void) +{ + int i; + float scale = 255.99609372; + + imdct_do_512 = imdct_do_512_c; + imdct_do_512_nol = imdct_do_512_nol_c; + fft_64p = fft_64p_c; + + /* Twiddle factors to turn IFFT into IMDCT */ + + for (i=0; i < 128; i++) { + xcos1[i] = cos(2.0f * M_PI * (8*i+1)/(8*N)) * scale; + xsin1[i] = sin(2.0f * M_PI * (8*i+1)/(8*N)) * scale; + } + + return 0; +} + + +static void imdct_do_512_c (float data[], float delay[]) +{ + int i, j; + float tmp_a_r, tmp_a_i; + float *data_ptr; + float *delay_ptr; + float *window_ptr; + +// 512 IMDCT with source and dest data in 'data' +// Pre IFFT complex multiply plus IFFT complex conjugate + + for( i=0; i < 128; i++) { + j = pm128[i]; + //a = (data[256-2*j-1] - data[2*j]) * (xcos1[j] + xsin1[j]); + //c = data[2*j] * xcos1[j]; + //b = data[256-2*j-1] * xsin1[j]; + //buf1[i].re = a - b + c; + //buf1[i].im = b + c; + buf[i].re = (data[256-2*j-1] * xcos1[j]) - (data[2*j] * xsin1[j]); + buf[i].im = -1.0 * (data[2*j] * xcos1[j] + data[256-2*j-1] * xsin1[j]); + } + + fft_128p_c (&buf[0]); + +// Post IFFT complex multiply plus IFFT complex conjugate + for (i=0; i < 128; i++) { + tmp_a_r = buf[i].re; + tmp_a_i = buf[i].im; + //a = (tmp_a_r - tmp_a_i) * (xcos1[j] + xsin1[j]); + //b = tmp_a_r * xsin1[j]; + //c = tmp_a_i * xcos1[j]; + //buf[j].re = a - b + c; + //buf[j].im = b + c; + buf[i].re =(tmp_a_r * xcos1[i]) + (tmp_a_i * xsin1[i]); + buf[i].im =(tmp_a_r * xsin1[i]) - (tmp_a_i * xcos1[i]); + } + + data_ptr = data; + delay_ptr = delay; + window_ptr = window; + +// Window and convert to real valued signal + for (i=0; i< 64; i++) { + *data_ptr++ = -buf[64+i].im * *window_ptr++ + *delay_ptr++; + *data_ptr++ = buf[64-i-1].re * *window_ptr++ + *delay_ptr++; + } + + for(i=0; i< 64; i++) { + *data_ptr++ = -buf[i].re * *window_ptr++ + *delay_ptr++; + *data_ptr++ = buf[128-i-1].im * *window_ptr++ + *delay_ptr++; + } + +// The trailing edge of the window goes into the delay line + delay_ptr = delay; + + for(i=0; i< 64; i++) { + *delay_ptr++ = -buf[64+i].re * *--window_ptr; + *delay_ptr++ = buf[64-i-1].im * *--window_ptr; + } + + for(i=0; i<64; i++) { + *delay_ptr++ = buf[i].im * *--window_ptr; + *delay_ptr++ = -buf[128-i-1].re * *--window_ptr; + } +} + + +static void imdct_do_512_nol_c (float data[], float delay[]) +{ + int i, j; + + float tmp_a_i; + float tmp_a_r; + + float *data_ptr; + float *delay_ptr; + float *window_ptr; + + // + // 512 IMDCT with source and dest data in 'data' + // + + // Pre IFFT complex multiply plus IFFT cmplx conjugate + + for( i=0; i < 128; i++) { + /* z[i] = (X[256-2*i-1] + j * X[2*i]) * (xcos1[i] + j * xsin1[i]) */ + j = pm128[i]; + //a = (data[256-2*j-1] - data[2*j]) * (xcos1[j] + xsin1[j]); + //c = data[2*j] * xcos1[j]; + //b = data[256-2*j-1] * xsin1[j]; + //buf1[i].re = a - b + c; + + //buf1[i].im = b + c; + buf[i].re = (data[256-2*j-1] * xcos1[j]) - (data[2*j] * xsin1[j]); + buf[i].im = -1.0 * (data[2*j] * xcos1[j] + data[256-2*j-1] * xsin1[j]); + } + + fft_128p_c (&buf[0]); + + /* Post IFFT complex multiply plus IFFT complex conjugate*/ + for (i=0; i < 128; i++) { + /* y[n] = z[n] * (xcos1[n] + j * xsin1[n]) ; */ + /* int j1 = i; */ + tmp_a_r = buf[i].re; + tmp_a_i = buf[i].im; + //a = (tmp_a_r - tmp_a_i) * (xcos1[j] + xsin1[j]); + //b = tmp_a_r * xsin1[j]; + //c = tmp_a_i * xcos1[j]; + //buf[j].re = a - b + c; + //buf[j].im = b + c; + buf[i].re =(tmp_a_r * xcos1[i]) + (tmp_a_i * xsin1[i]); + buf[i].im =(tmp_a_r * xsin1[i]) - (tmp_a_i * xcos1[i]); + } + + data_ptr = data; + delay_ptr = delay; + window_ptr = window; + + /* Window and convert to real valued signal, no overlap here*/ + for (i=0; i< 64; i++) { + *data_ptr++ = -buf[64+i].im * *window_ptr++; + *data_ptr++ = buf[64-i-1].re * *window_ptr++; + } + + for(i=0; i< 64; i++) { + *data_ptr++ = -buf[i].re * *window_ptr++; + *data_ptr++ = buf[128-i-1].im * *window_ptr++; + } + + /* The trailing edge of the window goes into the delay line */ + delay_ptr = delay; + + for(i=0; i< 64; i++) { + *delay_ptr++ = -buf[64+i].re * *--window_ptr; + *delay_ptr++ = buf[64-i-1].im * *--window_ptr; + } + + for(i=0; i<64; i++) { + *delay_ptr++ = buf[i].im * *--window_ptr; + *delay_ptr++ = -buf[128-i-1].re * *--window_ptr; + } +} |