git.sesse.net Git - vlc/blob - plugins/imdct/ac3_imdct_c.c

   1 /*****************************************************************************
   2  * ac3_imdct_c.c: ac3 DCT in C
   3  *****************************************************************************
   4  * Copyright (C) 1999, 2000 VideoLAN
   5  * $Id: ac3_imdct_c.c,v 1.1 2001/05/15 16:19:42 sam Exp $
   6  *
   7  * Authors: Renaud Dartus <reno@videolan.org>
   8  *          Aaron Holtzman <aholtzma@engr.uvic.ca>
   9  *
  10  * This program is free software; you can redistribute it and/or modify
  11  * it under the terms of the GNU General Public License as published by
  12  * the Free Software Foundation; either version 2 of the License, or
  13  * (at your option) any later version.
  14  *
  15  * This program is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  * GNU General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU General Public License
  21  * along with this program; if not, write to the Free Software
  22  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111, USA.
  23  *****************************************************************************/
  24
  25 #define MODULE_NAME imdct
  26 #include "modules_inner.h"
  27
  28 /*****************************************************************************
  29  * Preamble
  30  *****************************************************************************/
  31 #include "defs.h"
  32
  33 #include <string.h>                                              /* memcpy() */
  34
  35 #include <math.h>
  36 #include <stdio.h>
  37
  38 #include "config.h"
  39 #include "common.h"
  40 #include "threads.h"
  41 #include "mtime.h"
  42
  43 #include "ac3_imdct.h"
  44 #include "ac3_imdct_common.h"
  45
  46 #ifndef M_PI
  47 #   define M_PI 3.14159265358979323846
  48 #endif
  49
  50 void _M( fft_64p )  ( complex_t *x );
  51 void _M( fft_128p ) ( complex_t *x );
  52
  53 static float window[] = {
  54     0.00014, 0.00024, 0.00037, 0.00051, 0.00067, 0.00086, 0.00107, 0.00130,
  55     0.00157, 0.00187, 0.00220, 0.00256, 0.00297, 0.00341, 0.00390, 0.00443,
  56     0.00501, 0.00564, 0.00632, 0.00706, 0.00785, 0.00871, 0.00962, 0.01061,
  57     0.01166, 0.01279, 0.01399, 0.01526, 0.01662, 0.01806, 0.01959, 0.02121,
  58     0.02292, 0.02472, 0.02662, 0.02863, 0.03073, 0.03294, 0.03527, 0.03770,
  59     0.04025, 0.04292, 0.04571, 0.04862, 0.05165, 0.05481, 0.05810, 0.06153,
  60     0.06508, 0.06878, 0.07261, 0.07658, 0.08069, 0.08495, 0.08935, 0.09389,
  61     0.09859, 0.10343, 0.10842, 0.11356, 0.11885, 0.12429, 0.12988, 0.13563,
  62     0.14152, 0.14757, 0.15376, 0.16011, 0.16661, 0.17325, 0.18005, 0.18699,
  63     0.19407, 0.20130, 0.20867, 0.21618, 0.22382, 0.23161, 0.23952, 0.24757,
  64     0.25574, 0.26404, 0.27246, 0.28100, 0.28965, 0.29841, 0.30729, 0.31626,
  65     0.32533, 0.33450, 0.34376, 0.35311, 0.36253, 0.37204, 0.38161, 0.39126,
  66     0.40096, 0.41072, 0.42054, 0.43040, 0.44030, 0.45023, 0.46020, 0.47019,
  67     0.48020, 0.49022, 0.50025, 0.51028, 0.52031, 0.53033, 0.54033, 0.55031,
  68     0.56026, 0.57019, 0.58007, 0.58991, 0.59970, 0.60944, 0.61912, 0.62873,
  69     0.63827, 0.64774, 0.65713, 0.66643, 0.67564, 0.68476, 0.69377, 0.70269,
  70     0.71150, 0.72019, 0.72877, 0.73723, 0.74557, 0.75378, 0.76186, 0.76981,
  71     0.77762, 0.78530, 0.79283, 0.80022, 0.80747, 0.81457, 0.82151, 0.82831,
  72     0.83496, 0.84145, 0.84779, 0.85398, 0.86001, 0.86588, 0.87160, 0.87716,
  73     0.88257, 0.88782, 0.89291, 0.89785, 0.90264, 0.90728, 0.91176, 0.91610,
  74     0.92028, 0.92432, 0.92822, 0.93197, 0.93558, 0.93906, 0.94240, 0.94560,
  75     0.94867, 0.95162, 0.95444, 0.95713, 0.95971, 0.96217, 0.96451, 0.96674,
  76     0.96887, 0.97089, 0.97281, 0.97463, 0.97635, 0.97799, 0.97953, 0.98099,
  77     0.98236, 0.98366, 0.98488, 0.98602, 0.98710, 0.98811, 0.98905, 0.98994,
  78     0.99076, 0.99153, 0.99225, 0.99291, 0.99353, 0.99411, 0.99464, 0.99513,
  79     0.99558, 0.99600, 0.99639, 0.99674, 0.99706, 0.99736, 0.99763, 0.99788,
  80     0.99811, 0.99831, 0.99850, 0.99867, 0.99882, 0.99895, 0.99908, 0.99919,
  81     0.99929, 0.99938, 0.99946, 0.99953, 0.99959, 0.99965, 0.99969, 0.99974,
  82     0.99978, 0.99981, 0.99984, 0.99986, 0.99988, 0.99990, 0.99992, 0.99993,
  83     0.99994, 0.99995, 0.99996, 0.99997, 0.99998, 0.99998, 0.99998, 0.99999,
  84     0.99999, 0.99999, 0.99999, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
  85     1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000
  86 };
  87
  88 static const int pm128[128] =
  89 {
  90     0, 16, 32, 48, 64, 80,  96, 112,  8, 40, 72, 104, 24, 56,  88, 120,
  91     4, 20, 36, 52, 68, 84, 100, 116, 12, 28, 44,  60, 76, 92, 108, 124,
  92     2, 18, 34, 50, 66, 82,  98, 114, 10, 42, 74, 106, 26, 58,  90, 122,
  93     6, 22, 38, 54, 70, 86, 102, 118, 14, 46, 78, 110, 30, 62,  94, 126,
  94     1, 17, 33, 49, 65, 81,  97, 113,  9, 41, 73, 105, 25, 57,  89, 121,
  95     5, 21, 37, 53, 69, 85, 101, 117, 13, 29, 45,  61, 77, 93, 109, 125,
  96     3, 19, 35, 51, 67, 83,  99, 115, 11, 43, 75, 107, 27, 59,  91, 123,
  97     7, 23, 39, 55, 71, 87, 103, 119, 15, 31, 47,  63, 79, 95, 111, 127
  98 };
  99
 100 static const int pm64[64] =
 101 {
 102     0,  8, 16, 24, 32, 40, 48, 56,
 103     4, 20, 36, 52, 12, 28, 44, 60,
 104     2, 10, 18, 26, 34, 42, 50, 58,
 105     6, 14, 22, 30, 38, 46, 54, 62,
 106     1,  9, 17, 25, 33, 41, 49, 57,
 107     5, 21, 37, 53, 13, 29, 45, 61,
 108     3, 11, 19, 27, 35, 43, 51, 59,
 109     7, 23, 39, 55, 15, 31, 47, 63
 110 };
 111
 112 void _M( imdct_init ) (imdct_t * p_imdct)
 113 {
 114     int i;
 115     float scale = 181.019;
 116
 117     /* Twiddle factors to turn IFFT into IMDCT */
 118     for (i=0; i < 128; i++) {
 119         p_imdct->xcos1[i] = cos(2.0f * M_PI * (8*i+1)/(8*N)) * scale;
 120         p_imdct->xsin1[i] = sin(2.0f * M_PI * (8*i+1)/(8*N)) * scale;
 121     }
 122 }
 123
 124 void _M( imdct_do_512 ) (imdct_t * p_imdct, float data[], float delay[])
 125 {
 126     int i, j;
 127     float tmp_a_r, tmp_a_i;
 128     float *data_ptr;
 129     float *delay_ptr;
 130     float *window_ptr;
 131
 132     /* 512 IMDCT with source and dest data in 'data'
 133      * Pre IFFT complex multiply plus IFFT complex conjugate */
 134
 135     for( i=0; i < 128; i++) {
 136         j = pm128[i];
 137         /* a = (data[256-2*j-1] - data[2*j]) * (xcos1[j] + xsin1[j]);
 138          * c = data[2*j] * xcos1[j];
 139          * b = data[256-2*j-1] * xsin1[j];
 140          * buf1[i].real = a - b + c;
 141          * buf1[i].imag = b + c; */
 142         p_imdct->buf[i].real = (data[256-2*j-1] * p_imdct->xcos1[j]) - (data[2*j] * p_imdct->xsin1[j]);
 143         p_imdct->buf[i].imag = -1.0 * (data[2*j] * p_imdct->xcos1[j] + data[256-2*j-1] * p_imdct->xsin1[j]);
 144     }
 145
 146     _M( fft_128p ) ( &p_imdct->buf[0] );
 147
 148     /* Post IFFT complex multiply  plus IFFT complex conjugate */
 149     for (i=0; i < 128; i++) {
 150         tmp_a_r = p_imdct->buf[i].real;
 151         tmp_a_i = p_imdct->buf[i].imag;
 152         /* a = (tmp_a_r - tmp_a_i) * (xcos1[j] + xsin1[j]);
 153          * b = tmp_a_r * xsin1[j];
 154          * c = tmp_a_i * xcos1[j];
 155          * buf[j].real = a - b + c;
 156          * buf[j].imag = b + c; */
 157         p_imdct->buf[i].real =(tmp_a_r * p_imdct->xcos1[i])  +  (tmp_a_i  * p_imdct->xsin1[i]);
 158         p_imdct->buf[i].imag =(tmp_a_r * p_imdct->xsin1[i])  -  (tmp_a_i  * p_imdct->xcos1[i]);
 159     }
 160
 161     data_ptr = data;
 162     delay_ptr = delay;
 163     window_ptr = window;
 164
 165     /* Window and convert to real valued signal */
 166     for (i=0; i< 64; i++) {
 167         *data_ptr++ = -p_imdct->buf[64+i].imag  * *window_ptr++ + *delay_ptr++;
 168         *data_ptr++ = p_imdct->buf[64-i-1].real * *window_ptr++ + *delay_ptr++;
 169     }
 170
 171     for(i=0; i< 64; i++) {
 172         *data_ptr++ = -p_imdct->buf[i].real      * *window_ptr++ + *delay_ptr++;
 173         *data_ptr++ = p_imdct->buf[128-i-1].imag * *window_ptr++ + *delay_ptr++;
 174     }
 175
 176     /* The trailing edge of the window goes into the delay line */
 177     delay_ptr = delay;
 178
 179     for(i=0; i< 64; i++) {
 180         *delay_ptr++ = -p_imdct->buf[64+i].real   * *--window_ptr;
 181         *delay_ptr++ =  p_imdct->buf[64-i-1].imag * *--window_ptr;
 182     }
 183
 184     for(i=0; i<64; i++) {
 185         *delay_ptr++ =  p_imdct->buf[i].imag       * *--window_ptr;
 186         *delay_ptr++ = -p_imdct->buf[128-i-1].real * *--window_ptr;
 187     }
 188 }
 189
 190
 191 void _M( imdct_do_512_nol ) (imdct_t * p_imdct, float data[], float delay[])
 192 {
 193     int i, j;
 194
 195     float tmp_a_i;
 196     float tmp_a_r;
 197
 198     float *data_ptr;
 199     float *delay_ptr;
 200     float *window_ptr;
 201
 202     /* 512 IMDCT with source and dest data in 'data'
 203      * Pre IFFT complex multiply plus IFFT cmplx conjugate */
 204
 205     for( i=0; i < 128; i++) {
 206         /* z[i] = (X[256-2*i-1] + j * X[2*i]) * (xcos1[i] + j * xsin1[i]) */
 207         j = pm128[i];
 208         /* a = (data[256-2*j-1] - data[2*j]) * (xcos1[j] + xsin1[j]);
 209          * c = data[2*j] * xcos1[j];
 210          * b = data[256-2*j-1] * xsin1[j];
 211          * buf1[i].real = a - b + c;
 212          * buf1[i].imag = b + c; */
 213         p_imdct->buf[i].real = (data[256-2*j-1] * p_imdct->xcos1[j]) - (data[2*j] * p_imdct->xsin1[j]);
 214         p_imdct->buf[i].imag = -1.0 * (data[2*j] * p_imdct->xcos1[j] + data[256-2*j-1] * p_imdct->xsin1[j]);
 215     }
 216
 217     _M( fft_128p ) ( &p_imdct->buf[0] );
 218
 219     /* Post IFFT complex multiply  plus IFFT complex conjugate*/
 220     for (i=0; i < 128; i++) {
 221         /* y[n] = z[n] * (xcos1[n] + j * xsin1[n]) ;
 222          * int j1 = i; */
 223         tmp_a_r = p_imdct->buf[i].real;
 224         tmp_a_i = p_imdct->buf[i].imag;
 225         /* a = (tmp_a_r - tmp_a_i) * (xcos1[j] + xsin1[j]);
 226          * b = tmp_a_r * xsin1[j];
 227          * c = tmp_a_i * xcos1[j];
 228          * buf[j].real = a - b + c;
 229          * buf[j].imag = b + c; */
 230         p_imdct->buf[i].real =(tmp_a_r * p_imdct->xcos1[i]) + (tmp_a_i  * p_imdct->xsin1[i]);
 231         p_imdct->buf[i].imag =(tmp_a_r * p_imdct->xsin1[i]) - (tmp_a_i  * p_imdct->xcos1[i]);
 232     }
 233
 234     data_ptr = data;
 235     delay_ptr = delay;
 236     window_ptr = window;
 237
 238     /* Window and convert to real valued signal, no overlap here*/
 239     for (i=0; i< 64; i++) {
 240         *data_ptr++ = -p_imdct->buf[64+i].imag  * *window_ptr++;
 241         *data_ptr++ = p_imdct->buf[64-i-1].real * *window_ptr++;
 242     }
 243
 244     for(i=0; i< 64; i++) {
 245         *data_ptr++ = -p_imdct->buf[i].real      * *window_ptr++;
 246         *data_ptr++ = p_imdct->buf[128-i-1].imag * *window_ptr++;
 247     }
 248
 249     /* The trailing edge of the window goes into the delay line */
 250     delay_ptr = delay;
 251
 252     for(i=0; i< 64; i++) {
 253         *delay_ptr++ = -p_imdct->buf[64+i].real   * *--window_ptr;
 254         *delay_ptr++ =  p_imdct->buf[64-i-1].imag * *--window_ptr;
 255     }
 256
 257     for(i=0; i<64; i++) {
 258         *delay_ptr++ =  p_imdct->buf[i].imag       * *--window_ptr;
 259         *delay_ptr++ = -p_imdct->buf[128-i-1].real * *--window_ptr;
 260     }
 261 }
 262