git.sesse.net Git - vlc/blob - src/ac3_decoder/ac3_imdct_c.c

   1 /*****************************************************************************
   2  * ac3_imdct_c.c: ac3 DCT
   3  *****************************************************************************
   4  * Copyright (C) 1999, 2000 VideoLAN
   5  * $Id: ac3_imdct_c.c,v 1.1 2001/04/30 21:04:20 reno 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 #include "defs.h"
  26
  27 #include <math.h>
  28 #include <stdio.h>
  29
  30 #include "config.h"
  31 #include "common.h"
  32 #include "threads.h"
  33 #include "mtime.h"
  34
  35 #include "stream_control.h"
  36 #include "input_ext-dec.h"
  37
  38 #include "ac3_decoder.h"
  39 #include "ac3_internal.h"
  40
  41 void fft_64p_c (complex_t *x);
  42 void fft_128p_c (complex_t *x);
  43 void imdct_do_512_c (imdct_t * p_imdct, float data[], float delay[]);
  44 void imdct_do_512_nol_c (imdct_t * p_imdct, float data[], float delay[]);
  45
  46
  47 static float window[] = {
  48         0.00014, 0.00024, 0.00037, 0.00051, 0.00067, 0.00086, 0.00107, 0.00130,
  49         0.00157, 0.00187, 0.00220, 0.00256, 0.00297, 0.00341, 0.00390, 0.00443,
  50         0.00501, 0.00564, 0.00632, 0.00706, 0.00785, 0.00871, 0.00962, 0.01061,
  51         0.01166, 0.01279, 0.01399, 0.01526, 0.01662, 0.01806, 0.01959, 0.02121,
  52         0.02292, 0.02472, 0.02662, 0.02863, 0.03073, 0.03294, 0.03527, 0.03770,
  53         0.04025, 0.04292, 0.04571, 0.04862, 0.05165, 0.05481, 0.05810, 0.06153,
  54         0.06508, 0.06878, 0.07261, 0.07658, 0.08069, 0.08495, 0.08935, 0.09389,
  55         0.09859, 0.10343, 0.10842, 0.11356, 0.11885, 0.12429, 0.12988, 0.13563,
  56         0.14152, 0.14757, 0.15376, 0.16011, 0.16661, 0.17325, 0.18005, 0.18699,
  57         0.19407, 0.20130, 0.20867, 0.21618, 0.22382, 0.23161, 0.23952, 0.24757,
  58         0.25574, 0.26404, 0.27246, 0.28100, 0.28965, 0.29841, 0.30729, 0.31626,
  59         0.32533, 0.33450, 0.34376, 0.35311, 0.36253, 0.37204, 0.38161, 0.39126,
  60         0.40096, 0.41072, 0.42054, 0.43040, 0.44030, 0.45023, 0.46020, 0.47019,
  61         0.48020, 0.49022, 0.50025, 0.51028, 0.52031, 0.53033, 0.54033, 0.55031,
  62         0.56026, 0.57019, 0.58007, 0.58991, 0.59970, 0.60944, 0.61912, 0.62873,
  63         0.63827, 0.64774, 0.65713, 0.66643, 0.67564, 0.68476, 0.69377, 0.70269,
  64         0.71150, 0.72019, 0.72877, 0.73723, 0.74557, 0.75378, 0.76186, 0.76981,
  65         0.77762, 0.78530, 0.79283, 0.80022, 0.80747, 0.81457, 0.82151, 0.82831,
  66         0.83496, 0.84145, 0.84779, 0.85398, 0.86001, 0.86588, 0.87160, 0.87716,
  67         0.88257, 0.88782, 0.89291, 0.89785, 0.90264, 0.90728, 0.91176, 0.91610,
  68         0.92028, 0.92432, 0.92822, 0.93197, 0.93558, 0.93906, 0.94240, 0.94560,
  69         0.94867, 0.95162, 0.95444, 0.95713, 0.95971, 0.96217, 0.96451, 0.96674,
  70         0.96887, 0.97089, 0.97281, 0.97463, 0.97635, 0.97799, 0.97953, 0.98099,
  71         0.98236, 0.98366, 0.98488, 0.98602, 0.98710, 0.98811, 0.98905, 0.98994,
  72         0.99076, 0.99153, 0.99225, 0.99291, 0.99353, 0.99411, 0.99464, 0.99513,
  73         0.99558, 0.99600, 0.99639, 0.99674, 0.99706, 0.99736, 0.99763, 0.99788,
  74         0.99811, 0.99831, 0.99850, 0.99867, 0.99882, 0.99895, 0.99908, 0.99919,
  75         0.99929, 0.99938, 0.99946, 0.99953, 0.99959, 0.99965, 0.99969, 0.99974,
  76         0.99978, 0.99981, 0.99984, 0.99986, 0.99988, 0.99990, 0.99992, 0.99993,
  77         0.99994, 0.99995, 0.99996, 0.99997, 0.99998, 0.99998, 0.99998, 0.99999,
  78         0.99999, 0.99999, 0.99999, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
  79         1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000
  80 };
  81
  82 static const int pm128[128] =
  83 {
  84         0, 16, 32, 48, 64, 80,  96, 112,  8, 40, 72, 104, 24, 56,  88, 120,
  85         4, 20, 36, 52, 68, 84, 100, 116, 12, 28, 44,  60, 76, 92, 108, 124,
  86         2, 18, 34, 50, 66, 82,  98, 114, 10, 42, 74, 106, 26, 58,  90, 122,
  87         6, 22, 38, 54, 70, 86, 102, 118, 14, 46, 78, 110, 30, 62,  94, 126,
  88         1, 17, 33, 49, 65, 81,  97, 113,  9, 41, 73, 105, 25, 57,  89, 121,
  89         5, 21, 37, 53, 69, 85, 101, 117, 13, 29, 45,  61, 77, 93, 109, 125,
  90         3, 19, 35, 51, 67, 83,  99, 115, 11, 43, 75, 107, 27, 59,  91, 123,
  91         7, 23, 39, 55, 71, 87, 103, 119, 15, 31, 47,  63, 79, 95, 111, 127
  92 };
  93
  94 static const int pm64[64] =
  95 {
  96         0,  8, 16, 24, 32, 40, 48, 56,
  97         4, 20, 36, 52, 12, 28, 44, 60,
  98         2, 10, 18, 26, 34, 42, 50, 58,
  99         6, 14, 22, 30, 38, 46, 54, 62,
 100         1,  9, 17, 25, 33, 41, 49, 57,
 101         5, 21, 37, 53, 13, 29, 45, 61,
 102         3, 11, 19, 27, 35, 43, 51, 59,
 103         7, 23, 39, 55, 15, 31, 47, 63
 104 };
 105
 106 int imdct_init_c (imdct_t * p_imdct)
 107 {
 108         int i;
 109         float scale = 255.99609372;
 110
 111         p_imdct->imdct_do_512 = imdct_do_512_c;
 112         p_imdct->imdct_do_512_nol = imdct_do_512_nol_c;
 113         p_imdct->fft_64p = fft_64p_c;
 114
 115         /* Twiddle factors to turn IFFT into IMDCT */
 116
 117         for (i=0; i < 128; i++) {
 118                 p_imdct->xcos1[i] = cos(2.0f * M_PI * (8*i+1)/(8*N)) * scale;
 119                 p_imdct->xsin1[i] = sin(2.0f * M_PI * (8*i+1)/(8*N)) * scale;
 120         }
 121
 122         return 0;
 123 }
 124
 125 void imdct_do_256 (imdct_t * p_imdct, float data[],float delay[])
 126 {
 127         int i, j, k;
 128         int p, q;
 129
 130         float tmp_a_i;
 131         float tmp_a_r;
 132
 133         float *data_ptr;
 134         float *delay_ptr;
 135         float *window_ptr;
 136
 137         complex_t *buf1, *buf2;
 138
 139         buf1 = &p_imdct->buf[0];
 140         buf2 = &p_imdct->buf[64];
 141
 142     /* Pre IFFT complex multiply plus IFFT complex conjugate */
 143         for (k=0; k<64; k++) {
 144                 /* X1[k] = X[2*k]
 145                  * X2[k] = X[2*k+1]     */
 146
 147                 j = pm64[k];
 148                 p = 2 * (128-2*j-1);
 149                 q = 2 * (2 * j);
 150
 151                 /* Z1[k] = (X1[128-2*k-1] + j * X1[2*k]) * (xcos2[k] + j * xsin2[k]); */
 152                 buf1[k].real =        data[p] * p_imdct->xcos2[j] - data[q] * p_imdct->xsin2[j];
 153                 buf1[k].imag = -1.0f*(data[q] * p_imdct->xcos2[j] + data[p] * p_imdct->xsin2[j]);
 154                 /* Z2[k] = (X2[128-2*k-1] + j * X2[2*k]) * (xcos2[k] + j * xsin2[k]); */
 155                 buf2[k].real =        data[p + 1] * p_imdct->xcos2[j] - data[q + 1] * p_imdct->xsin2[j];
 156                 buf2[k].imag = -1.0f*(data[q + 1] * p_imdct->xcos2[j] + data[p + 1] * p_imdct->xsin2[j]);
 157         }
 158
 159         p_imdct->fft_64p(&buf1[0]);
 160         p_imdct->fft_64p(&buf2[0]);
 161
 162         /* Post IFFT complex multiply */
 163         for( i=0; i < 64; i++) {
 164                 tmp_a_r =  buf1[i].real;
 165                 tmp_a_i = -buf1[i].imag;
 166                 buf1[i].real = (tmp_a_r * p_imdct->xcos2[i]) - (tmp_a_i * p_imdct->xsin2[i]);
 167                 buf1[i].imag = (tmp_a_r * p_imdct->xsin2[i]) + (tmp_a_i * p_imdct->xcos2[i]);
 168                 tmp_a_r =  buf2[i].real;
 169                 tmp_a_i = -buf2[i].imag;
 170                 buf2[i].real = (tmp_a_r * p_imdct->xcos2[i]) - (tmp_a_i * p_imdct->xsin2[i]);
 171                 buf2[i].imag = (tmp_a_r * p_imdct->xsin2[i]) + (tmp_a_i * p_imdct->xcos2[i]);
 172         }
 173
 174         data_ptr = data;
 175         delay_ptr = delay;
 176         window_ptr = window;
 177
 178         /* Window and convert to real valued signal */
 179         for(i=0; i< 64; i++) {
 180                 *data_ptr++ = -buf1[i].imag     * *window_ptr++ + *delay_ptr++;
 181                 *data_ptr++ = buf1[64-i-1].real * *window_ptr++ + *delay_ptr++;
 182         }
 183
 184         for(i=0; i< 64; i++) {
 185                 *data_ptr++ = -buf1[i].real     * *window_ptr++ + *delay_ptr++;
 186                 *data_ptr++ = buf1[64-i-1].imag * *window_ptr++ + *delay_ptr++;
 187         }
 188
 189         delay_ptr = delay;
 190
 191         for(i=0; i< 64; i++) {
 192                 *delay_ptr++ = -buf2[i].real      * *--window_ptr;
 193                 *delay_ptr++ =  buf2[64-i-1].imag * *--window_ptr;
 194         }
 195
 196         for(i=0; i< 64; i++) {
 197                 *delay_ptr++ =  buf2[i].imag      * *--window_ptr;
 198                 *delay_ptr++ = -buf2[64-i-1].real * *--window_ptr;
 199         }
 200 }
 201
 202
 203 void imdct_do_256_nol (imdct_t * p_imdct, float data[], float delay[])
 204 {
 205         int i, j, k;
 206         int p, q;
 207
 208         float tmp_a_i;
 209         float tmp_a_r;
 210
 211         float *data_ptr;
 212         float *delay_ptr;
 213         float *window_ptr;
 214
 215         complex_t *buf1, *buf2;
 216
 217         buf1 = &p_imdct->buf[0];
 218         buf2 = &p_imdct->buf[64];
 219
 220     /* Pre IFFT complex multiply plus IFFT cmplx conjugate */
 221         for(k=0; k<64; k++) {
 222         /* X1[k] = X[2*k]
 223         * X2[k] = X[2*k+1] */
 224         j = pm64[k];
 225         p = 2 * (128-2*j-1);
 226         q = 2 * (2 * j);
 227
 228         /* Z1[k] = (X1[128-2*k-1] + j * X1[2*k]) * (xcos2[k] + j * xsin2[k]); */
 229         buf1[k].real =        data[p] * p_imdct->xcos2[j] - data[q] * p_imdct->xsin2[j];
 230         buf1[k].imag = -1.0f*(data[q] * p_imdct->xcos2[j] + data[p] * p_imdct->xsin2[j]);
 231         /* Z2[k] = (X2[128-2*k-1] + j * X2[2*k]) * (xcos2[k] + j * xsin2[k]); */
 232         buf2[k].real =        data[p + 1] * p_imdct->xcos2[j] - data[q + 1] * p_imdct->xsin2[j];
 233         buf2[k].imag = -1.0f*(data[q + 1] * p_imdct->xcos2[j] + data[p + 1] * p_imdct->xsin2[j]);
 234     }
 235
 236     p_imdct->fft_64p(&buf1[0]);
 237     p_imdct->fft_64p(&buf2[0]);
 238
 239     /* Post IFFT complex multiply */
 240     for( i=0; i < 64; i++) {
 241         /* y1[n] = z1[n] * (xcos2[n] + j * xs in2[n]) ; */
 242         tmp_a_r =  buf1[i].real;
 243         tmp_a_i = -buf1[i].imag;
 244         buf1[i].real =(tmp_a_r * p_imdct->xcos2[i])  -  (tmp_a_i  * p_imdct->xsin2[i]);
 245         buf1[i].imag =(tmp_a_r * p_imdct->xsin2[i])  +  (tmp_a_i  * p_imdct->xcos2[i]);
 246         /* y2[n] = z2[n] * (xcos2[n] + j * xsin2[n]) ; */
 247         tmp_a_r =  buf2[i].real;
 248         tmp_a_i = -buf2[i].imag;
 249         buf2[i].real =(tmp_a_r * p_imdct->xcos2[i])  -  (tmp_a_i  * p_imdct->xsin2[i]);
 250         buf2[i].imag =(tmp_a_r * p_imdct->xsin2[i])  +  (tmp_a_i  * p_imdct->xcos2[i]);
 251     }
 252
 253     data_ptr = data;
 254     delay_ptr = delay;
 255     window_ptr = window;
 256
 257     /* Window and convert to real valued signal, no overlap */
 258     for(i=0; i< 64; i++) {
 259         *data_ptr++ = -buf1[i].imag     * *window_ptr++;
 260         *data_ptr++ = buf1[64-i-1].real * *window_ptr++;
 261     }
 262
 263     for(i=0; i< 64; i++) {
 264         *data_ptr++ = -buf1[i].real     * *window_ptr++ + *delay_ptr++;
 265         *data_ptr++ = buf1[64-i-1].imag * *window_ptr++ + *delay_ptr++;
 266     }
 267
 268     delay_ptr = delay;
 269
 270     for(i=0; i< 64; i++) {
 271         *delay_ptr++ = -buf2[i].real      * *--window_ptr;
 272         *delay_ptr++ =  buf2[64-i-1].imag * *--window_ptr;
 273     }
 274
 275     for(i=0; i< 64; i++) {
 276         *delay_ptr++ =  buf2[i].imag      * *--window_ptr;
 277         *delay_ptr++ = -buf2[64-i-1].real * *--window_ptr;
 278         }
 279 }
 280
 281 void imdct_do_512_c (imdct_t * p_imdct, float data[], float delay[])
 282 {
 283         int i, j;
 284         float tmp_a_r, tmp_a_i;
 285         float *data_ptr;
 286         float *delay_ptr;
 287         float *window_ptr;
 288
 289     /* 512 IMDCT with source and dest data in 'data'
 290      * Pre IFFT complex multiply plus IFFT complex conjugate */
 291
 292     for( i=0; i < 128; i++) {
 293                 j = pm128[i];
 294                 /* a = (data[256-2*j-1] - data[2*j]) * (xcos1[j] + xsin1[j]);
 295                  * c = data[2*j] * xcos1[j];
 296                  * b = data[256-2*j-1] * xsin1[j];
 297                  * buf1[i].real = a - b + c;
 298                  * buf1[i].imag = b + c; */
 299                 p_imdct->buf[i].real = (data[256-2*j-1] * p_imdct->xcos1[j]) - (data[2*j] * p_imdct->xsin1[j]);
 300                 p_imdct->buf[i].imag = -1.0 * (data[2*j] * p_imdct->xcos1[j] + data[256-2*j-1] * p_imdct->xsin1[j]);
 301         }
 302
 303         fft_128p_c (&p_imdct->buf[0]);
 304
 305     /* Post IFFT complex multiply  plus IFFT complex conjugate */
 306         for (i=0; i < 128; i++) {
 307                 tmp_a_r = p_imdct->buf[i].real;
 308                 tmp_a_i = p_imdct->buf[i].imag;
 309                 /* a = (tmp_a_r - tmp_a_i) * (xcos1[j] + xsin1[j]);
 310                  * b = tmp_a_r * xsin1[j];
 311                  * c = tmp_a_i * xcos1[j];
 312                  * buf[j].real = a - b + c;
 313                  * buf[j].imag = b + c; */
 314                 p_imdct->buf[i].real =(tmp_a_r * p_imdct->xcos1[i])  +  (tmp_a_i  * p_imdct->xsin1[i]);
 315                 p_imdct->buf[i].imag =(tmp_a_r * p_imdct->xsin1[i])  -  (tmp_a_i  * p_imdct->xcos1[i]);
 316         }
 317
 318         data_ptr = data;
 319         delay_ptr = delay;
 320         window_ptr = window;
 321
 322     /* Window and convert to real valued signal */
 323         for (i=0; i< 64; i++) {
 324                 *data_ptr++ = -p_imdct->buf[64+i].imag  * *window_ptr++ + *delay_ptr++;
 325                 *data_ptr++ = p_imdct->buf[64-i-1].real * *window_ptr++ + *delay_ptr++;
 326         }
 327
 328         for(i=0; i< 64; i++) {
 329                 *data_ptr++ = -p_imdct->buf[i].real      * *window_ptr++ + *delay_ptr++;
 330                 *data_ptr++ = p_imdct->buf[128-i-1].imag * *window_ptr++ + *delay_ptr++;
 331         }
 332
 333     /* The trailing edge of the window goes into the delay line */
 334         delay_ptr = delay;
 335
 336         for(i=0; i< 64; i++) {
 337                 *delay_ptr++ = -p_imdct->buf[64+i].real   * *--window_ptr;
 338                 *delay_ptr++ =  p_imdct->buf[64-i-1].imag * *--window_ptr;
 339         }
 340
 341         for(i=0; i<64; i++) {
 342                 *delay_ptr++ =  p_imdct->buf[i].imag       * *--window_ptr;
 343                 *delay_ptr++ = -p_imdct->buf[128-i-1].real * *--window_ptr;
 344         }
 345 }
 346
 347
 348 void imdct_do_512_nol_c (imdct_t * p_imdct, float data[], float delay[])
 349 {
 350         int i, j;
 351
 352         float tmp_a_i;
 353         float tmp_a_r;
 354
 355         float *data_ptr;
 356         float *delay_ptr;
 357         float *window_ptr;
 358
 359     /* 512 IMDCT with source and dest data in 'data'
 360          * Pre IFFT complex multiply plus IFFT cmplx conjugate */
 361
 362     for( i=0; i < 128; i++) {
 363         /* z[i] = (X[256-2*i-1] + j * X[2*i]) * (xcos1[i] + j * xsin1[i]) */
 364                 j = pm128[i];
 365         /* a = (data[256-2*j-1] - data[2*j]) * (xcos1[j] + xsin1[j]);
 366          * c = data[2*j] * xcos1[j];
 367          * b = data[256-2*j-1] * xsin1[j];
 368          * buf1[i].real = a - b + c;
 369          * buf1[i].imag = b + c; */
 370                 p_imdct->buf[i].real = (data[256-2*j-1] * p_imdct->xcos1[j]) - (data[2*j] * p_imdct->xsin1[j]);
 371                 p_imdct->buf[i].imag = -1.0 * (data[2*j] * p_imdct->xcos1[j] + data[256-2*j-1] * p_imdct->xsin1[j]);
 372         }
 373
 374         fft_128p_c (&p_imdct->buf[0]);
 375
 376     /* Post IFFT complex multiply  plus IFFT complex conjugate*/
 377         for (i=0; i < 128; i++) {
 378                 /* y[n] = z[n] * (xcos1[n] + j * xsin1[n]) ;
 379                  * int j1 = i; */
 380                 tmp_a_r = p_imdct->buf[i].real;
 381                 tmp_a_i = p_imdct->buf[i].imag;
 382                 /* a = (tmp_a_r - tmp_a_i) * (xcos1[j] + xsin1[j]);
 383                  * b = tmp_a_r * xsin1[j];
 384                  * c = tmp_a_i * xcos1[j];
 385                  * buf[j].real = a - b + c;
 386                  * buf[j].imag = b + c; */
 387                 p_imdct->buf[i].real =(tmp_a_r * p_imdct->xcos1[i]) + (tmp_a_i  * p_imdct->xsin1[i]);
 388                 p_imdct->buf[i].imag =(tmp_a_r * p_imdct->xsin1[i]) - (tmp_a_i  * p_imdct->xcos1[i]);
 389         }
 390
 391         data_ptr = data;
 392         delay_ptr = delay;
 393         window_ptr = window;
 394
 395         /* Window and convert to real valued signal, no overlap here*/
 396         for (i=0; i< 64; i++) {
 397                 *data_ptr++ = -p_imdct->buf[64+i].imag  * *window_ptr++;
 398                 *data_ptr++ = p_imdct->buf[64-i-1].real * *window_ptr++;
 399         }
 400
 401         for(i=0; i< 64; i++) {
 402                 *data_ptr++ = -p_imdct->buf[i].real      * *window_ptr++;
 403                 *data_ptr++ = p_imdct->buf[128-i-1].imag * *window_ptr++;
 404         }
 405
 406         /* The trailing edge of the window goes into the delay line */
 407         delay_ptr = delay;
 408
 409         for(i=0; i< 64; i++) {
 410                 *delay_ptr++ = -p_imdct->buf[64+i].real   * *--window_ptr;
 411                 *delay_ptr++ =  p_imdct->buf[64-i-1].imag * *--window_ptr;
 412         }
 413
 414         for(i=0; i<64; i++) {
 415                 *delay_ptr++ =  p_imdct->buf[i].imag       * *--window_ptr;
 416                 *delay_ptr++ = -p_imdct->buf[128-i-1].real * *--window_ptr;
 417         }
 418 }