lib/bladeenc/subs.c

   1 /*
   2                         (c) Copyright 1998-2001 - Tord Jansson
   3                         ======================================
   4
   5                 This file is part of the BladeEnc MP3 Encoder, based on
   6                 ISO's reference code for MPEG Layer 3 compression, and might
   7                 contain smaller or larger sections that are directly taken
   8                 from ISO's reference code.
   9
  10                 All changes to the ISO reference code herein are either
  11                 copyrighted by Tord Jansson (tord.jansson@swipnet.se)
  12                 or sublicensed to Tord Jansson by a third party.
  13
  14         BladeEnc is free software; you can redistribute this file
  15         and/or modify it under the terms of the GNU Lesser General Public
  16         License as published by the Free Software Foundation; either
  17         version 2.1 of the License, or (at your option) any later version.
  18
  19
  20
  21         ------------    Changes    ------------
  22
  23         2000-11-04  Andre Piotrowski
  24
  25         -       speed up: don«t calculate the unneeded phi values!
  26
  27         2000-11-22      ap
  28
  29         -       speed up: slightly improved way to handle the 'special case'
  30
  31         2000-12-05  ap
  32
  33         -       speed up: implemented prepacking
  34
  35         2000-12-11  ap
  36
  37         -       speed up: avoid data reordering
  38
  39         2000-12-12  ap
  40
  41         -       moved fft configuration switches to "l3psy.h"
  42
  43         2001-01-12  ap
  44
  45         -       use some explicit type casting to avoid compiler warnings
  46 */
  47
  48 /*****************************************************************************
  49  * FFT computes fast fourier transform of BLKSIZE samples of data            *
  50  *   based on the decimation-in-frequency algorithm described in "Digital    *
  51  *   Signal Processing" by Oppenheim and Schafer, refer to pages 304         *
  52  *   (flow graph) and 330-332 (Fortran program in problem 5).                *
  53  *                                                                           *
  54  *   required constants:                                                     *
  55  *         PI          3.14159265358979                                      *
  56  *         BLKSIZE     should be 2^^M for a positive integer M               *
  57  *                                                                           *
  58  *****************************************************************************/
  59
  60 #include        "common.h"
  61 #include        "encoder.h"
  62 #include        "l3psy.h"
  63
  64
  65
  66
  67
  68 /* The switches have been moved to "l3psy.h" */
  69
  70
  71
  72
  73
  74 int                                             fInit_fft;
  75
  76
  77
  78
  79
  80 void fft (FLOAT x_real[], FLOAT x_imag[], FLOAT energy[], FLOAT phi[], int N)
  81 {
  82         int                                             i, j, k, off;
  83 #if USED_VALUES_ONLY
  84         int                                             min_i, max_i;
  85 #endif
  86         int                                             ip, le, le1;
  87         static double                   w_real[BLKSIZE/2], w_imag[BLKSIZE/2];
  88         double                                  t_real, t_imag, u_real, u_imag;
  89         int                                             N_ORG;
  90
  91 #if !REORDER_DATA
  92         static  int                             swap_l[BLKSIZE/2+1];
  93         static  int                             swap_s[BLKSIZE_s/2+1];
  94         int                                             *pSwap, a, b;
  95 #endif
  96
  97         double  t1, t2, t3, t4, t5, t6;
  98
  99         if (fInit_fft == 0)
 100         {
 101                 for (i = 0;  i < BLKSIZE/2;  i++)
 102                 {
 103                         w_real[i] =  cos(PI*i/(BLKSIZE/2));
 104 #if NORMAL_FFT
 105                         w_imag[i] = -sin(PI*i/(BLKSIZE/2));
 106 #else
 107                         w_imag[i] =  sin(PI*i/(BLKSIZE/2));
 108 #endif
 109                 }
 110
 111 #if !REORDER_DATA
 112                 j = 0;
 113                 for (i = 0;  i < BLKSIZE/2-1;  i++)
 114                 {
 115                         swap_l[i] = j;  k = BLKSIZE/4;  while (k <= j) {j -= k;  k >>= 1;}  j += k;
 116                 }
 117                 swap_l[i] = i;  swap_l[i+1] = i+1;
 118
 119                 j = 0;
 120                 for (i = 0;  i < BLKSIZE_s/2-1;  i++)
 121                 {
 122                         swap_s[i] = j;  k = BLKSIZE_s/4;  while (k <= j) {j -= k;  k >>= 1;}  j += k;
 123                 }
 124                 swap_s[i] = i;  swap_s[i+1] = i+1;
 125 #endif
 126
 127                 fInit_fft++;
 128         }
 129
 130
 131 #if REAL_SEQUENCE
 132         N_ORG = N;
 133         N >>= 1;
 134 #if !PREPACKED
 135         /* packing the sequence to the half length */
 136         for (i = 0;  i < N;  i++)
 137         {
 138                 x_real[i] = x_real[2*i];
 139                 x_imag[i] = x_real[2*i+1];
 140         }
 141 #endif
 142 #endif
 143
 144
 145         off = BLKSIZE/N;
 146         for (le = N;  le > 1;  le = le1)
 147         {
 148                 le1 = le >> 1;
 149
 150                         /* special case: k=0; u_real=1.0; u_imag=0.0; j=0 */
 151                         for (i = 0;  i < N;  i += le)
 152                         {
 153                                 ip = i + le1;
 154                                 t_real = x_real[i] - x_real[ip];
 155                                 t_imag = x_imag[i] - x_imag[ip];
 156                                 x_real[i] += x_real[ip];
 157                                 x_imag[i] += x_imag[ip];
 158                                 x_real[ip] = (FLOAT) t_real;
 159                                 x_imag[ip] = (FLOAT) t_imag;
 160                         }
 161
 162                 k = off;
 163                 for (j = 1;  j < le1;  j++)
 164                 {
 165                         u_real = w_real[k];
 166                         u_imag = w_imag[k];
 167                         for (i = j;  i < N;  i += le)
 168                         {
 169                                 ip = i + le1;
 170                                 t_real = x_real[i] - x_real[ip];
 171                                 t_imag = x_imag[i] - x_imag[ip];
 172                                 x_real[i] += x_real[ip];
 173                                 x_imag[i] += x_imag[ip];
 174                                 x_real[ip] = (FLOAT) (t_real*u_real - t_imag*u_imag);
 175                                 x_imag[ip] = (FLOAT) (t_imag*u_real + t_real*u_imag);
 176                         }
 177                         k += off;
 178                 }
 179
 180                 off += off;
 181         }
 182
 183
 184 #if REORDER_DATA
 185         /* this section reorders the data to the correct ordering */
 186         j = 0;
 187         for (i = 0;  i < N-1;  i++)
 188         {
 189                 if (i < j)
 190                 {
 191                         t_real = x_real[j];
 192                         t_imag = x_imag[j];
 193                         x_real[j] = x_real[i];
 194                         x_imag[j] = x_imag[i];
 195                         x_real[i] = (FLOAT) t_real;
 196                         x_imag[i] = (FLOAT) t_imag;
 197                 }
 198                 k=N/2;
 199                 while (k <= j)
 200                 {
 201                         j -= k;
 202                         k >>= 1;
 203                 }
 204                 j += k;
 205         }
 206 #else
 207         /*
 208                 We don«t reorder the data to the correct ordering,
 209                 but access the data by the bit reverse order index array.
 210         */
 211         pSwap = (N_ORG == BLKSIZE) ? swap_l : swap_s;
 212 #endif
 213
 214
 215 #if REAL_SEQUENCE
 216         /* unpacking the sequence */
 217         t_real = x_real[0];
 218         t_imag = x_imag[0];
 219         x_real[0] = (FLOAT) (t_real+t_imag);
 220         x_imag[0] = 0.0;
 221         x_real[N] = (FLOAT) (t_real-t_imag);
 222         x_imag[N] = 0.0;
 223
 224         k = off = BLKSIZE/N_ORG;
 225         for (i = 1;  i < N/2;  i++)
 226         {
 227 #if REORDER_DATA
 228 #define         a       i
 229 #define         b       (N-i)
 230 #else
 231                 a = pSwap[i];
 232                 b = pSwap[N-i];
 233 #endif
 234                 t1 = x_real[a] + x_real[b];
 235                 t2 = x_real[a] - x_real[b];
 236                 t3 = x_imag[a] + x_imag[b];
 237                 t4 = x_imag[a] - x_imag[b];
 238                 t5 = t2*w_imag[k] + t3*w_real[k];
 239                 t6 = t3*w_imag[k] - t2*w_real[k];
 240
 241                 x_real[a] = (FLOAT) (t1+t5) / 2.0;
 242                 x_imag[a] = (FLOAT) (t6+t4) / 2.0;
 243                 x_real[b] = (FLOAT) (t1-t5) / 2.0;
 244                 x_imag[b] = (FLOAT) (t6-t4) / 2.0;
 245
 246                 k += off;
 247         }
 248         /* x_real[N/2] doesn«t change */
 249         /* x_real[N/2] changes the sign in case of a normal fft */
 250 #if (NORMAL_FFT)
 251 #if REORDER_DATA
 252         x_imag[i] = -x_imag[i];
 253 #else
 254         x_imag[pSwap[i]] *= -1.0;
 255 #endif   /* REORDER_DATA */
 256 #endif   /* NORMAL_FFT */
 257         N = N_ORG;
 258 #endif   /* REAL_SEQUENCE */
 259
 260
 261         /* calculating the energy and phase, phi */
 262 #if USED_VALUES_ONLY
 263         if (N == BLKSIZE)
 264         {
 265                 min_i = LONG_FFT_MIN_IDX;
 266                 max_i = LONG_FFT_MAX_IDX;
 267         }
 268         else
 269         {
 270                 min_i = SHORT_FFT_MIN_IDX;
 271                 max_i = SHORT_FFT_MAX_IDX;
 272         }
 273 #endif
 274 #if REAL_SEQUENCE
 275         for (i = 0;  i <= N/2;  i++)
 276 #else
 277         for (i = 0;  i < N;  i++)
 278 #endif
 279         {
 280 #if REORDER_DATA
 281 #define         a       i
 282 #else
 283                 a = pSwap[i];
 284 #endif
 285                 energy[i] = x_real[a]*x_real[a] + x_imag[a]*x_imag[a];
 286                 if(energy[i] <= 0.0005)
 287                 {
 288                         energy[i] = (FLOAT) 0.0005;        /* keep the identity */
 289                         x_real[a] = (FLOAT) sqrt(0.0005);  /* energy[i] * cos(phi[i]) == x_real[i] */
 290                         x_imag[a] = 0.0;                   /* energy[i] * sin(phi[i]) == x_imag[i] */
 291                 }
 292 #if USED_VALUES_ONLY
 293                 if (i >= min_i  &&  i <= max_i)
 294 #endif
 295                 phi[i] = (FLOAT) atan2((double) x_imag[a], (double) x_real[a]);
 296         }
 297 }
 298
 299
 300