src/spectral/ooura_fft8g.c

   1 // modifications made for aubio:
   2 //  - replace all 'double' with 'smpl_t'
   3 //  - include "aubio_priv.h" (for config.h and types.h)
   4 //  - add missing prototypes
   5 //  - use COS, SIN, and ATAN macros
   6 //  - add cast to (smpl_t) to avoid float conversion warnings
   7 //  - declare initialization as static
   8 //  - prefix public function with aubio_ooura_
   9
  10 #include "aubio_priv.h"
  11
  12 void aubio_ooura_cdft(int n, int isgn, smpl_t *a, int *ip, smpl_t *w);
  13 void aubio_ooura_rdft(int n, int isgn, smpl_t *a, int *ip, smpl_t *w);
  14 void aubio_ooura_ddct(int n, int isgn, smpl_t *a, int *ip, smpl_t *w);
  15 void aubio_ooura_ddst(int n, int isgn, smpl_t *a, int *ip, smpl_t *w);
  16 void aubio_ooura_dfct(int n, smpl_t *a, smpl_t *t, int *ip, smpl_t *w);
  17 void aubio_ooura_dfst(int n, smpl_t *a, smpl_t *t, int *ip, smpl_t *w);
  18 static void makewt(int nw, int *ip, smpl_t *w);
  19 static void makect(int nc, int *ip, smpl_t *c);
  20 static void bitrv2(int n, int *ip, smpl_t *a);
  21 static void bitrv2conj(int n, int *ip, smpl_t *a);
  22 static void cftfsub(int n, smpl_t *a, smpl_t *w);
  23 static void cftbsub(int n, smpl_t *a, smpl_t *w);
  24 static void cft1st(int n, smpl_t *a, smpl_t *w);
  25 static void cftmdl(int n, int l, smpl_t *a, smpl_t *w);
  26 static void rftfsub(int n, smpl_t *a, int nc, smpl_t *c);
  27 static void rftbsub(int n, smpl_t *a, int nc, smpl_t *c);
  28 static void dctsub(int n, smpl_t *a, int nc, smpl_t *c);
  29 static void dstsub(int n, smpl_t *a, int nc, smpl_t *c);
  30
  31 /*
  32 Fast Fourier/Cosine/Sine Transform
  33     dimension   :one
  34     data length :power of 2
  35     decimation  :frequency
  36     radix       :8, 4, 2
  37     data        :inplace
  38     table       :use
  39 functions
  40     cdft: Complex Discrete Fourier Transform
  41     rdft: Real Discrete Fourier Transform
  42     ddct: Discrete Cosine Transform
  43     ddst: Discrete Sine Transform
  44     dfct: Cosine Transform of RDFT (Real Symmetric DFT)
  45     dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
  46 function prototypes
  47     void cdft(int, int, smpl_t *, int *, smpl_t *);
  48     void rdft(int, int, smpl_t *, int *, smpl_t *);
  49     void ddct(int, int, smpl_t *, int *, smpl_t *);
  50     void ddst(int, int, smpl_t *, int *, smpl_t *);
  51     void dfct(int, smpl_t *, smpl_t *, int *, smpl_t *);
  52     void dfst(int, smpl_t *, smpl_t *, int *, smpl_t *);
  53
  54
  55 -------- Complex DFT (Discrete Fourier Transform) --------
  56     [definition]
  57         <case1>
  58             X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
  59         <case2>
  60             X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
  61         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
  62     [usage]
  63         <case1>
  64             ip[0] = 0; // first time only
  65             cdft(2*n, 1, a, ip, w);
  66         <case2>
  67             ip[0] = 0; // first time only
  68             cdft(2*n, -1, a, ip, w);
  69     [parameters]
  70         2*n            :data length (int)
  71                         n >= 1, n = power of 2
  72         a[0...2*n-1]   :input/output data (smpl_t *)
  73                         input data
  74                             a[2*j] = Re(x[j]),
  75                             a[2*j+1] = Im(x[j]), 0<=j<n
  76                         output data
  77                             a[2*k] = Re(X[k]),
  78                             a[2*k+1] = Im(X[k]), 0<=k<n
  79         ip[0...*]      :work area for bit reversal (int *)
  80                         length of ip >= 2+sqrt(n)
  81                         strictly,
  82                         length of ip >=
  83                             2+(1<<(int)(log(n+0.5)/log(2))/2).
  84                         ip[0],ip[1] are pointers of the cos/sin table.
  85         w[0...n/2-1]   :cos/sin table (smpl_t *)
  86                         w[],ip[] are initialized if ip[0] == 0.
  87     [remark]
  88         Inverse of
  89             cdft(2*n, -1, a, ip, w);
  90         is
  91             cdft(2*n, 1, a, ip, w);
  92             for (j = 0; j <= 2 * n - 1; j++) {
  93                 a[j] *= 1.0 / n;
  94             }
  95         .
  96
  97
  98 -------- Real DFT / Inverse of Real DFT --------
  99     [definition]
 100         <case1> RDFT
 101             R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
 102             I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
 103         <case2> IRDFT (excluding scale)
 104             a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
 105                    sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
 106                    sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
 107     [usage]
 108         <case1>
 109             ip[0] = 0; // first time only
 110             rdft(n, 1, a, ip, w);
 111         <case2>
 112             ip[0] = 0; // first time only
 113             rdft(n, -1, a, ip, w);
 114     [parameters]
 115         n              :data length (int)
 116                         n >= 2, n = power of 2
 117         a[0...n-1]     :input/output data (smpl_t *)
 118                         <case1>
 119                             output data
 120                                 a[2*k] = R[k], 0<=k<n/2
 121                                 a[2*k+1] = I[k], 0<k<n/2
 122                                 a[1] = R[n/2]
 123                         <case2>
 124                             input data
 125                                 a[2*j] = R[j], 0<=j<n/2
 126                                 a[2*j+1] = I[j], 0<j<n/2
 127                                 a[1] = R[n/2]
 128         ip[0...*]      :work area for bit reversal (int *)
 129                         length of ip >= 2+sqrt(n/2)
 130                         strictly,
 131                         length of ip >=
 132                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
 133                         ip[0],ip[1] are pointers of the cos/sin table.
 134         w[0...n/2-1]   :cos/sin table (smpl_t *)
 135                         w[],ip[] are initialized if ip[0] == 0.
 136     [remark]
 137         Inverse of
 138             rdft(n, 1, a, ip, w);
 139         is
 140             rdft(n, -1, a, ip, w);
 141             for (j = 0; j <= n - 1; j++) {
 142                 a[j] *= 2.0 / n;
 143             }
 144         .
 145
 146
 147 -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
 148     [definition]
 149         <case1> IDCT (excluding scale)
 150             C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
 151         <case2> DCT
 152             C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
 153     [usage]
 154         <case1>
 155             ip[0] = 0; // first time only
 156             ddct(n, 1, a, ip, w);
 157         <case2>
 158             ip[0] = 0; // first time only
 159             ddct(n, -1, a, ip, w);
 160     [parameters]
 161         n              :data length (int)
 162                         n >= 2, n = power of 2
 163         a[0...n-1]     :input/output data (smpl_t *)
 164                         output data
 165                             a[k] = C[k], 0<=k<n
 166         ip[0...*]      :work area for bit reversal (int *)
 167                         length of ip >= 2+sqrt(n/2)
 168                         strictly,
 169                         length of ip >=
 170                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
 171                         ip[0],ip[1] are pointers of the cos/sin table.
 172         w[0...n*5/4-1] :cos/sin table (smpl_t *)
 173                         w[],ip[] are initialized if ip[0] == 0.
 174     [remark]
 175         Inverse of
 176             ddct(n, -1, a, ip, w);
 177         is
 178             a[0] *= 0.5;
 179             ddct(n, 1, a, ip, w);
 180             for (j = 0; j <= n - 1; j++) {
 181                 a[j] *= 2.0 / n;
 182             }
 183         .
 184
 185
 186 -------- DST (Discrete Sine Transform) / Inverse of DST --------
 187     [definition]
 188         <case1> IDST (excluding scale)
 189             S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
 190         <case2> DST
 191             S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
 192     [usage]
 193         <case1>
 194             ip[0] = 0; // first time only
 195             ddst(n, 1, a, ip, w);
 196         <case2>
 197             ip[0] = 0; // first time only
 198             ddst(n, -1, a, ip, w);
 199     [parameters]
 200         n              :data length (int)
 201                         n >= 2, n = power of 2
 202         a[0...n-1]     :input/output data (smpl_t *)
 203                         <case1>
 204                             input data
 205                                 a[j] = A[j], 0<j<n
 206                                 a[0] = A[n]
 207                             output data
 208                                 a[k] = S[k], 0<=k<n
 209                         <case2>
 210                             output data
 211                                 a[k] = S[k], 0<k<n
 212                                 a[0] = S[n]
 213         ip[0...*]      :work area for bit reversal (int *)
 214                         length of ip >= 2+sqrt(n/2)
 215                         strictly,
 216                         length of ip >=
 217                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
 218                         ip[0],ip[1] are pointers of the cos/sin table.
 219         w[0...n*5/4-1] :cos/sin table (smpl_t *)
 220                         w[],ip[] are initialized if ip[0] == 0.
 221     [remark]
 222         Inverse of
 223             ddst(n, -1, a, ip, w);
 224         is
 225             a[0] *= 0.5;
 226             ddst(n, 1, a, ip, w);
 227             for (j = 0; j <= n - 1; j++) {
 228                 a[j] *= 2.0 / n;
 229             }
 230         .
 231
 232
 233 -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
 234     [definition]
 235         C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
 236     [usage]
 237         ip[0] = 0; // first time only
 238         dfct(n, a, t, ip, w);
 239     [parameters]
 240         n              :data length - 1 (int)
 241                         n >= 2, n = power of 2
 242         a[0...n]       :input/output data (smpl_t *)
 243                         output data
 244                             a[k] = C[k], 0<=k<=n
 245         t[0...n/2]     :work area (smpl_t *)
 246         ip[0...*]      :work area for bit reversal (int *)
 247                         length of ip >= 2+sqrt(n/4)
 248                         strictly,
 249                         length of ip >=
 250                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
 251                         ip[0],ip[1] are pointers of the cos/sin table.
 252         w[0...n*5/8-1] :cos/sin table (smpl_t *)
 253                         w[],ip[] are initialized if ip[0] == 0.
 254     [remark]
 255         Inverse of
 256             a[0] *= 0.5;
 257             a[n] *= 0.5;
 258             dfct(n, a, t, ip, w);
 259         is
 260             a[0] *= 0.5;
 261             a[n] *= 0.5;
 262             dfct(n, a, t, ip, w);
 263             for (j = 0; j <= n; j++) {
 264                 a[j] *= 2.0 / n;
 265             }
 266         .
 267
 268
 269 -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
 270     [definition]
 271         S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
 272     [usage]
 273         ip[0] = 0; // first time only
 274         dfst(n, a, t, ip, w);
 275     [parameters]
 276         n              :data length + 1 (int)
 277                         n >= 2, n = power of 2
 278         a[0...n-1]     :input/output data (smpl_t *)
 279                         output data
 280                             a[k] = S[k], 0<k<n
 281                         (a[0] is used for work area)
 282         t[0...n/2-1]   :work area (smpl_t *)
 283         ip[0...*]      :work area for bit reversal (int *)
 284                         length of ip >= 2+sqrt(n/4)
 285                         strictly,
 286                         length of ip >=
 287                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
 288                         ip[0],ip[1] are pointers of the cos/sin table.
 289         w[0...n*5/8-1] :cos/sin table (smpl_t *)
 290                         w[],ip[] are initialized if ip[0] == 0.
 291     [remark]
 292         Inverse of
 293             dfst(n, a, t, ip, w);
 294         is
 295             dfst(n, a, t, ip, w);
 296             for (j = 1; j <= n - 1; j++) {
 297                 a[j] *= 2.0 / n;
 298             }
 299         .
 300
 301
 302 Appendix :
 303     The cos/sin table is recalculated when the larger table required.
 304     w[] and ip[] are compatible with all routines.
 305 */
 306
 307
 308 void aubio_ooura_cdft(int n, int isgn, smpl_t *a, int *ip, smpl_t *w)
 309 {
 310     void makewt(int nw, int *ip, smpl_t *w);
 311     void bitrv2(int n, int *ip, smpl_t *a);
 312     void bitrv2conj(int n, int *ip, smpl_t *a);
 313     void cftfsub(int n, smpl_t *a, smpl_t *w);
 314     void cftbsub(int n, smpl_t *a, smpl_t *w);
 315
 316     if (n > (ip[0] << 2)) {
 317         makewt(n >> 2, ip, w);
 318     }
 319     if (n > 4) {
 320         if (isgn >= 0) {
 321             bitrv2(n, ip + 2, a);
 322             cftfsub(n, a, w);
 323         } else {
 324             bitrv2conj(n, ip + 2, a);
 325             cftbsub(n, a, w);
 326         }
 327     } else if (n == 4) {
 328         cftfsub(n, a, w);
 329     }
 330 }
 331
 332
 333 void aubio_ooura_rdft(int n, int isgn, smpl_t *a, int *ip, smpl_t *w)
 334 {
 335     void makewt(int nw, int *ip, smpl_t *w);
 336     void makect(int nc, int *ip, smpl_t *c);
 337     void bitrv2(int n, int *ip, smpl_t *a);
 338     void cftfsub(int n, smpl_t *a, smpl_t *w);
 339     void cftbsub(int n, smpl_t *a, smpl_t *w);
 340     void rftfsub(int n, smpl_t *a, int nc, smpl_t *c);
 341     void rftbsub(int n, smpl_t *a, int nc, smpl_t *c);
 342     int nw, nc;
 343     smpl_t xi;
 344
 345     nw = ip[0];
 346     if (n > (nw << 2)) {
 347         nw = n >> 2;
 348         makewt(nw, ip, w);
 349     }
 350     nc = ip[1];
 351     if (n > (nc << 2)) {
 352         nc = n >> 2;
 353         makect(nc, ip, w + nw);
 354     }
 355     if (isgn >= 0) {
 356         if (n > 4) {
 357             bitrv2(n, ip + 2, a);
 358             cftfsub(n, a, w);
 359             rftfsub(n, a, nc, w + nw);
 360         } else if (n == 4) {
 361             cftfsub(n, a, w);
 362         }
 363         xi = a[0] - a[1];
 364         a[0] += a[1];
 365         a[1] = xi;
 366     } else {
 367         a[1] = (smpl_t)0.5 * (a[0] - a[1]);
 368         a[0] -= a[1];
 369         if (n > 4) {
 370             rftbsub(n, a, nc, w + nw);
 371             bitrv2(n, ip + 2, a);
 372             cftbsub(n, a, w);
 373         } else if (n == 4) {
 374             cftfsub(n, a, w);
 375         }
 376     }
 377 }
 378
 379
 380 void aubio_ooura_ddct(int n, int isgn, smpl_t *a, int *ip, smpl_t *w)
 381 {
 382     void makewt(int nw, int *ip, smpl_t *w);
 383     void makect(int nc, int *ip, smpl_t *c);
 384     void bitrv2(int n, int *ip, smpl_t *a);
 385     void cftfsub(int n, smpl_t *a, smpl_t *w);
 386     void cftbsub(int n, smpl_t *a, smpl_t *w);
 387     void rftfsub(int n, smpl_t *a, int nc, smpl_t *c);
 388     void rftbsub(int n, smpl_t *a, int nc, smpl_t *c);
 389     void dctsub(int n, smpl_t *a, int nc, smpl_t *c);
 390     int j, nw, nc;
 391     smpl_t xr;
 392
 393     nw = ip[0];
 394     if (n > (nw << 2)) {
 395         nw = n >> 2;
 396         makewt(nw, ip, w);
 397     }
 398     nc = ip[1];
 399     if (n > nc) {
 400         nc = n;
 401         makect(nc, ip, w + nw);
 402     }
 403     if (isgn < 0) {
 404         xr = a[n - 1];
 405         for (j = n - 2; j >= 2; j -= 2) {
 406             a[j + 1] = a[j] - a[j - 1];
 407             a[j] += a[j - 1];
 408         }
 409         a[1] = a[0] - xr;
 410         a[0] += xr;
 411         if (n > 4) {
 412             rftbsub(n, a, nc, w + nw);
 413             bitrv2(n, ip + 2, a);
 414             cftbsub(n, a, w);
 415         } else if (n == 4) {
 416             cftfsub(n, a, w);
 417         }
 418     }
 419     dctsub(n, a, nc, w + nw);
 420     if (isgn >= 0) {
 421         if (n > 4) {
 422             bitrv2(n, ip + 2, a);
 423             cftfsub(n, a, w);
 424             rftfsub(n, a, nc, w + nw);
 425         } else if (n == 4) {
 426             cftfsub(n, a, w);
 427         }
 428         xr = a[0] - a[1];
 429         a[0] += a[1];
 430         for (j = 2; j < n; j += 2) {
 431             a[j - 1] = a[j] - a[j + 1];
 432             a[j] += a[j + 1];
 433         }
 434         a[n - 1] = xr;
 435     }
 436 }
 437
 438
 439 void aubio_ooura_ddst(int n, int isgn, smpl_t *a, int *ip, smpl_t *w)
 440 {
 441     void makewt(int nw, int *ip, smpl_t *w);
 442     void makect(int nc, int *ip, smpl_t *c);
 443     void bitrv2(int n, int *ip, smpl_t *a);
 444     void cftfsub(int n, smpl_t *a, smpl_t *w);
 445     void cftbsub(int n, smpl_t *a, smpl_t *w);
 446     void rftfsub(int n, smpl_t *a, int nc, smpl_t *c);
 447     void rftbsub(int n, smpl_t *a, int nc, smpl_t *c);
 448     void dstsub(int n, smpl_t *a, int nc, smpl_t *c);
 449     int j, nw, nc;
 450     smpl_t xr;
 451
 452     nw = ip[0];
 453     if (n > (nw << 2)) {
 454         nw = n >> 2;
 455         makewt(nw, ip, w);
 456     }
 457     nc = ip[1];
 458     if (n > nc) {
 459         nc = n;
 460         makect(nc, ip, w + nw);
 461     }
 462     if (isgn < 0) {
 463         xr = a[n - 1];
 464         for (j = n - 2; j >= 2; j -= 2) {
 465             a[j + 1] = -a[j] - a[j - 1];
 466             a[j] -= a[j - 1];
 467         }
 468         a[1] = a[0] + xr;
 469         a[0] -= xr;
 470         if (n > 4) {
 471             rftbsub(n, a, nc, w + nw);
 472             bitrv2(n, ip + 2, a);
 473             cftbsub(n, a, w);
 474         } else if (n == 4) {
 475             cftfsub(n, a, w);
 476         }
 477     }
 478     dstsub(n, a, nc, w + nw);
 479     if (isgn >= 0) {
 480         if (n > 4) {
 481             bitrv2(n, ip + 2, a);
 482             cftfsub(n, a, w);
 483             rftfsub(n, a, nc, w + nw);
 484         } else if (n == 4) {
 485             cftfsub(n, a, w);
 486         }
 487         xr = a[0] - a[1];
 488         a[0] += a[1];
 489         for (j = 2; j < n; j += 2) {
 490             a[j - 1] = -a[j] - a[j + 1];
 491             a[j] -= a[j + 1];
 492         }
 493         a[n - 1] = -xr;
 494     }
 495 }
 496
 497
 498 void aubio_ooura_dfct(int n, smpl_t *a, smpl_t *t, int *ip, smpl_t *w)
 499 {
 500     void makewt(int nw, int *ip, smpl_t *w);
 501     void makect(int nc, int *ip, smpl_t *c);
 502     void bitrv2(int n, int *ip, smpl_t *a);
 503     void cftfsub(int n, smpl_t *a, smpl_t *w);
 504     void rftfsub(int n, smpl_t *a, int nc, smpl_t *c);
 505     void dctsub(int n, smpl_t *a, int nc, smpl_t *c);
 506     int j, k, l, m, mh, nw, nc;
 507     smpl_t xr, xi, yr, yi;
 508
 509     nw = ip[0];
 510     if (n > (nw << 3)) {
 511         nw = n >> 3;
 512         makewt(nw, ip, w);
 513     }
 514     nc = ip[1];
 515     if (n > (nc << 1)) {
 516         nc = n >> 1;
 517         makect(nc, ip, w + nw);
 518     }
 519     m = n >> 1;
 520     yi = a[m];
 521     xi = a[0] + a[n];
 522     a[0] -= a[n];
 523     t[0] = xi - yi;
 524     t[m] = xi + yi;
 525     if (n > 2) {
 526         mh = m >> 1;
 527         for (j = 1; j < mh; j++) {
 528             k = m - j;
 529             xr = a[j] - a[n - j];
 530             xi = a[j] + a[n - j];
 531             yr = a[k] - a[n - k];
 532             yi = a[k] + a[n - k];
 533             a[j] = xr;
 534             a[k] = yr;
 535             t[j] = xi - yi;
 536             t[k] = xi + yi;
 537         }
 538         t[mh] = a[mh] + a[n - mh];
 539         a[mh] -= a[n - mh];
 540         dctsub(m, a, nc, w + nw);
 541         if (m > 4) {
 542             bitrv2(m, ip + 2, a);
 543             cftfsub(m, a, w);
 544             rftfsub(m, a, nc, w + nw);
 545         } else if (m == 4) {
 546             cftfsub(m, a, w);
 547         }
 548         a[n - 1] = a[0] - a[1];
 549         a[1] = a[0] + a[1];
 550         for (j = m - 2; j >= 2; j -= 2) {
 551             a[2 * j + 1] = a[j] + a[j + 1];
 552             a[2 * j - 1] = a[j] - a[j + 1];
 553         }
 554         l = 2;
 555         m = mh;
 556         while (m >= 2) {
 557             dctsub(m, t, nc, w + nw);
 558             if (m > 4) {
 559                 bitrv2(m, ip + 2, t);
 560                 cftfsub(m, t, w);
 561                 rftfsub(m, t, nc, w + nw);
 562             } else if (m == 4) {
 563                 cftfsub(m, t, w);
 564             }
 565             a[n - l] = t[0] - t[1];
 566             a[l] = t[0] + t[1];
 567             k = 0;
 568             for (j = 2; j < m; j += 2) {
 569                 k += l << 2;
 570                 a[k - l] = t[j] - t[j + 1];
 571                 a[k + l] = t[j] + t[j + 1];
 572             }
 573             l <<= 1;
 574             mh = m >> 1;
 575             for (j = 0; j < mh; j++) {
 576                 k = m - j;
 577                 t[j] = t[m + k] - t[m + j];
 578                 t[k] = t[m + k] + t[m + j];
 579             }
 580             t[mh] = t[m + mh];
 581             m = mh;
 582         }
 583         a[l] = t[0];
 584         a[n] = t[2] - t[1];
 585         a[0] = t[2] + t[1];
 586     } else {
 587         a[1] = a[0];
 588         a[2] = t[0];
 589         a[0] = t[1];
 590     }
 591 }
 592
 593
 594 void aubio_ooura_dfst(int n, smpl_t *a, smpl_t *t, int *ip, smpl_t *w)
 595 {
 596     void makewt(int nw, int *ip, smpl_t *w);
 597     void makect(int nc, int *ip, smpl_t *c);
 598     void bitrv2(int n, int *ip, smpl_t *a);
 599     void cftfsub(int n, smpl_t *a, smpl_t *w);
 600     void rftfsub(int n, smpl_t *a, int nc, smpl_t *c);
 601     void dstsub(int n, smpl_t *a, int nc, smpl_t *c);
 602     int j, k, l, m, mh, nw, nc;
 603     smpl_t xr, xi, yr, yi;
 604
 605     nw = ip[0];
 606     if (n > (nw << 3)) {
 607         nw = n >> 3;
 608         makewt(nw, ip, w);
 609     }
 610     nc = ip[1];
 611     if (n > (nc << 1)) {
 612         nc = n >> 1;
 613         makect(nc, ip, w + nw);
 614     }
 615     if (n > 2) {
 616         m = n >> 1;
 617         mh = m >> 1;
 618         for (j = 1; j < mh; j++) {
 619             k = m - j;
 620             xr = a[j] + a[n - j];
 621             xi = a[j] - a[n - j];
 622             yr = a[k] + a[n - k];
 623             yi = a[k] - a[n - k];
 624             a[j] = xr;
 625             a[k] = yr;
 626             t[j] = xi + yi;
 627             t[k] = xi - yi;
 628         }
 629         t[0] = a[mh] - a[n - mh];
 630         a[mh] += a[n - mh];
 631         a[0] = a[m];
 632         dstsub(m, a, nc, w + nw);
 633         if (m > 4) {
 634             bitrv2(m, ip + 2, a);
 635             cftfsub(m, a, w);
 636             rftfsub(m, a, nc, w + nw);
 637         } else if (m == 4) {
 638             cftfsub(m, a, w);
 639         }
 640         a[n - 1] = a[1] - a[0];
 641         a[1] = a[0] + a[1];
 642         for (j = m - 2; j >= 2; j -= 2) {
 643             a[2 * j + 1] = a[j] - a[j + 1];
 644             a[2 * j - 1] = -a[j] - a[j + 1];
 645         }
 646         l = 2;
 647         m = mh;
 648         while (m >= 2) {
 649             dstsub(m, t, nc, w + nw);
 650             if (m > 4) {
 651                 bitrv2(m, ip + 2, t);
 652                 cftfsub(m, t, w);
 653                 rftfsub(m, t, nc, w + nw);
 654             } else if (m == 4) {
 655                 cftfsub(m, t, w);
 656             }
 657             a[n - l] = t[1] - t[0];
 658             a[l] = t[0] + t[1];
 659             k = 0;
 660             for (j = 2; j < m; j += 2) {
 661                 k += l << 2;
 662                 a[k - l] = -t[j] - t[j + 1];
 663                 a[k + l] = t[j] - t[j + 1];
 664             }
 665             l <<= 1;
 666             mh = m >> 1;
 667             for (j = 1; j < mh; j++) {
 668                 k = m - j;
 669                 t[j] = t[m + k] + t[m + j];
 670                 t[k] = t[m + k] - t[m + j];
 671             }
 672             t[0] = t[m + mh];
 673             m = mh;
 674         }
 675         a[l] = t[0];
 676     }
 677     a[0] = 0;
 678 }
 679
 680
 681 /* -------- initializing routines -------- */
 682
 683
 684 #include <math.h>
 685
 686 void makewt(int nw, int *ip, smpl_t *w)
 687 {
 688     void bitrv2(int n, int *ip, smpl_t *a);
 689     int j, nwh;
 690     smpl_t delta, x, y;
 691
 692     ip[0] = nw;
 693     ip[1] = 1;
 694     if (nw > 2) {
 695         nwh = nw >> 1;
 696         delta = ATAN(1.0) / nwh;
 697         w[0] = 1;
 698         w[1] = 0;
 699         w[nwh] = COS(delta * nwh);
 700         w[nwh + 1] = w[nwh];
 701         if (nwh > 2) {
 702             for (j = 2; j < nwh; j += 2) {
 703                 x = COS(delta * j);
 704                 y = SIN(delta * j);
 705                 w[j] = x;
 706                 w[j + 1] = y;
 707                 w[nw - j] = y;
 708                 w[nw - j + 1] = x;
 709             }
 710             for (j = nwh - 2; j >= 2; j -= 2) {
 711                 x = w[2 * j];
 712                 y = w[2 * j + 1];
 713                 w[nwh + j] = x;
 714                 w[nwh + j + 1] = y;
 715             }
 716             bitrv2(nw, ip + 2, w);
 717         }
 718     }
 719 }
 720
 721
 722 void makect(int nc, int *ip, smpl_t *c)
 723 {
 724     int j, nch;
 725     smpl_t delta;
 726
 727     ip[1] = nc;
 728     if (nc > 1) {
 729         nch = nc >> 1;
 730         delta = ATAN(1.0) / nch;
 731         c[0] = COS(delta * nch);
 732         c[nch] = (smpl_t)0.5 * c[0];
 733         for (j = 1; j < nch; j++) {
 734             c[j] = (smpl_t)0.5 * COS(delta * j);
 735             c[nc - j] = (smpl_t)0.5 * SIN(delta * j);
 736         }
 737     }
 738 }
 739
 740
 741 /* -------- child routines -------- */
 742
 743
 744 void bitrv2(int n, int *ip, smpl_t *a)
 745 {
 746     int j, j1, k, k1, l, m, m2;
 747     smpl_t xr, xi, yr, yi;
 748
 749     ip[0] = 0;
 750     l = n;
 751     m = 1;
 752     while ((m << 3) < l) {
 753         l >>= 1;
 754         for (j = 0; j < m; j++) {
 755             ip[m + j] = ip[j] + l;
 756         }
 757         m <<= 1;
 758     }
 759     m2 = 2 * m;
 760     if ((m << 3) == l) {
 761         for (k = 0; k < m; k++) {
 762             for (j = 0; j < k; j++) {
 763                 j1 = 2 * j + ip[k];
 764                 k1 = 2 * k + ip[j];
 765                 xr = a[j1];
 766                 xi = a[j1 + 1];
 767                 yr = a[k1];
 768                 yi = a[k1 + 1];
 769                 a[j1] = yr;
 770                 a[j1 + 1] = yi;
 771                 a[k1] = xr;
 772                 a[k1 + 1] = xi;
 773                 j1 += m2;
 774                 k1 += 2 * m2;
 775                 xr = a[j1];
 776                 xi = a[j1 + 1];
 777                 yr = a[k1];
 778                 yi = a[k1 + 1];
 779                 a[j1] = yr;
 780                 a[j1 + 1] = yi;
 781                 a[k1] = xr;
 782                 a[k1 + 1] = xi;
 783                 j1 += m2;
 784                 k1 -= m2;
 785                 xr = a[j1];
 786                 xi = a[j1 + 1];
 787                 yr = a[k1];
 788                 yi = a[k1 + 1];
 789                 a[j1] = yr;
 790                 a[j1 + 1] = yi;
 791                 a[k1] = xr;
 792                 a[k1 + 1] = xi;
 793                 j1 += m2;
 794                 k1 += 2 * m2;
 795                 xr = a[j1];
 796                 xi = a[j1 + 1];
 797                 yr = a[k1];
 798                 yi = a[k1 + 1];
 799                 a[j1] = yr;
 800                 a[j1 + 1] = yi;
 801                 a[k1] = xr;
 802                 a[k1 + 1] = xi;
 803             }
 804             j1 = 2 * k + m2 + ip[k];
 805             k1 = j1 + m2;
 806             xr = a[j1];
 807             xi = a[j1 + 1];
 808             yr = a[k1];
 809             yi = a[k1 + 1];
 810             a[j1] = yr;
 811             a[j1 + 1] = yi;
 812             a[k1] = xr;
 813             a[k1 + 1] = xi;
 814         }
 815     } else {
 816         for (k = 1; k < m; k++) {
 817             for (j = 0; j < k; j++) {
 818                 j1 = 2 * j + ip[k];
 819                 k1 = 2 * k + ip[j];
 820                 xr = a[j1];
 821                 xi = a[j1 + 1];
 822                 yr = a[k1];
 823                 yi = a[k1 + 1];
 824                 a[j1] = yr;
 825                 a[j1 + 1] = yi;
 826                 a[k1] = xr;
 827                 a[k1 + 1] = xi;
 828                 j1 += m2;
 829                 k1 += m2;
 830                 xr = a[j1];
 831                 xi = a[j1 + 1];
 832                 yr = a[k1];
 833                 yi = a[k1 + 1];
 834                 a[j1] = yr;
 835                 a[j1 + 1] = yi;
 836                 a[k1] = xr;
 837                 a[k1 + 1] = xi;
 838             }
 839         }
 840     }
 841 }
 842
 843
 844 void bitrv2conj(int n, int *ip, smpl_t *a)
 845 {
 846     int j, j1, k, k1, l, m, m2;
 847     smpl_t xr, xi, yr, yi;
 848
 849     ip[0] = 0;
 850     l = n;
 851     m = 1;
 852     while ((m << 3) < l) {
 853         l >>= 1;
 854         for (j = 0; j < m; j++) {
 855             ip[m + j] = ip[j] + l;
 856         }
 857         m <<= 1;
 858     }
 859     m2 = 2 * m;
 860     if ((m << 3) == l) {
 861         for (k = 0; k < m; k++) {
 862             for (j = 0; j < k; j++) {
 863                 j1 = 2 * j + ip[k];
 864                 k1 = 2 * k + ip[j];
 865                 xr = a[j1];
 866                 xi = -a[j1 + 1];
 867                 yr = a[k1];
 868                 yi = -a[k1 + 1];
 869                 a[j1] = yr;
 870                 a[j1 + 1] = yi;
 871                 a[k1] = xr;
 872                 a[k1 + 1] = xi;
 873                 j1 += m2;
 874                 k1 += 2 * m2;
 875                 xr = a[j1];
 876                 xi = -a[j1 + 1];
 877                 yr = a[k1];
 878                 yi = -a[k1 + 1];
 879                 a[j1] = yr;
 880                 a[j1 + 1] = yi;
 881                 a[k1] = xr;
 882                 a[k1 + 1] = xi;
 883                 j1 += m2;
 884                 k1 -= m2;
 885                 xr = a[j1];
 886                 xi = -a[j1 + 1];
 887                 yr = a[k1];
 888                 yi = -a[k1 + 1];
 889                 a[j1] = yr;
 890                 a[j1 + 1] = yi;
 891                 a[k1] = xr;
 892                 a[k1 + 1] = xi;
 893                 j1 += m2;
 894                 k1 += 2 * m2;
 895                 xr = a[j1];
 896                 xi = -a[j1 + 1];
 897                 yr = a[k1];
 898                 yi = -a[k1 + 1];
 899                 a[j1] = yr;
 900                 a[j1 + 1] = yi;
 901                 a[k1] = xr;
 902                 a[k1 + 1] = xi;
 903             }
 904             k1 = 2 * k + ip[k];
 905             a[k1 + 1] = -a[k1 + 1];
 906             j1 = k1 + m2;
 907             k1 = j1 + m2;
 908             xr = a[j1];
 909             xi = -a[j1 + 1];
 910             yr = a[k1];
 911             yi = -a[k1 + 1];
 912             a[j1] = yr;
 913             a[j1 + 1] = yi;
 914             a[k1] = xr;
 915             a[k1 + 1] = xi;
 916             k1 += m2;
 917             a[k1 + 1] = -a[k1 + 1];
 918         }
 919     } else {
 920         a[1] = -a[1];
 921         a[m2 + 1] = -a[m2 + 1];
 922         for (k = 1; k < m; k++) {
 923             for (j = 0; j < k; j++) {
 924                 j1 = 2 * j + ip[k];
 925                 k1 = 2 * k + ip[j];
 926                 xr = a[j1];
 927                 xi = -a[j1 + 1];
 928                 yr = a[k1];
 929                 yi = -a[k1 + 1];
 930                 a[j1] = yr;
 931                 a[j1 + 1] = yi;
 932                 a[k1] = xr;
 933                 a[k1 + 1] = xi;
 934                 j1 += m2;
 935                 k1 += m2;
 936                 xr = a[j1];
 937                 xi = -a[j1 + 1];
 938                 yr = a[k1];
 939                 yi = -a[k1 + 1];
 940                 a[j1] = yr;
 941                 a[j1 + 1] = yi;
 942                 a[k1] = xr;
 943                 a[k1 + 1] = xi;
 944             }
 945             k1 = 2 * k + ip[k];
 946             a[k1 + 1] = -a[k1 + 1];
 947             a[k1 + m2 + 1] = -a[k1 + m2 + 1];
 948         }
 949     }
 950 }
 951
 952
 953 void cftfsub(int n, smpl_t *a, smpl_t *w)
 954 {
 955     void cft1st(int n, smpl_t *a, smpl_t *w);
 956     void cftmdl(int n, int l, smpl_t *a, smpl_t *w);
 957     int j, j1, j2, j3, l;
 958     smpl_t x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
 959
 960     l = 2;
 961     if (n >= 16) {
 962         cft1st(n, a, w);
 963         l = 16;
 964         while ((l << 3) <= n) {
 965             cftmdl(n, l, a, w);
 966             l <<= 3;
 967         }
 968     }
 969     if ((l << 1) < n) {
 970         for (j = 0; j < l; j += 2) {
 971             j1 = j + l;
 972             j2 = j1 + l;
 973             j3 = j2 + l;
 974             x0r = a[j] + a[j1];
 975             x0i = a[j + 1] + a[j1 + 1];
 976             x1r = a[j] - a[j1];
 977             x1i = a[j + 1] - a[j1 + 1];
 978             x2r = a[j2] + a[j3];
 979             x2i = a[j2 + 1] + a[j3 + 1];
 980             x3r = a[j2] - a[j3];
 981             x3i = a[j2 + 1] - a[j3 + 1];
 982             a[j] = x0r + x2r;
 983             a[j + 1] = x0i + x2i;
 984             a[j2] = x0r - x2r;
 985             a[j2 + 1] = x0i - x2i;
 986             a[j1] = x1r - x3i;
 987             a[j1 + 1] = x1i + x3r;
 988             a[j3] = x1r + x3i;
 989             a[j3 + 1] = x1i - x3r;
 990         }
 991     } else if ((l << 1) == n) {
 992         for (j = 0; j < l; j += 2) {
 993             j1 = j + l;
 994             x0r = a[j] - a[j1];
 995             x0i = a[j + 1] - a[j1 + 1];
 996             a[j] += a[j1];
 997             a[j + 1] += a[j1 + 1];
 998             a[j1] = x0r;
 999             a[j1 + 1] = x0i;
1000         }
1001     }
1002 }
1003
1004
1005 void cftbsub(int n, smpl_t *a, smpl_t *w)
1006 {
1007     void cft1st(int n, smpl_t *a, smpl_t *w);
1008     void cftmdl(int n, int l, smpl_t *a, smpl_t *w);
1009     int j, j1, j2, j3, j4, j5, j6, j7, l;
1010     smpl_t wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1011         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
1012         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
1013
1014     l = 2;
1015     if (n > 16) {
1016         cft1st(n, a, w);
1017         l = 16;
1018         while ((l << 3) < n) {
1019             cftmdl(n, l, a, w);
1020             l <<= 3;
1021         }
1022     }
1023     if ((l << 2) < n) {
1024         wn4r = w[2];
1025         for (j = 0; j < l; j += 2) {
1026             j1 = j + l;
1027             j2 = j1 + l;
1028             j3 = j2 + l;
1029             j4 = j3 + l;
1030             j5 = j4 + l;
1031             j6 = j5 + l;
1032             j7 = j6 + l;
1033             x0r = a[j] + a[j1];
1034             x0i = -a[j + 1] - a[j1 + 1];
1035             x1r = a[j] - a[j1];
1036             x1i = -a[j + 1] + a[j1 + 1];
1037             x2r = a[j2] + a[j3];
1038             x2i = a[j2 + 1] + a[j3 + 1];
1039             x3r = a[j2] - a[j3];
1040             x3i = a[j2 + 1] - a[j3 + 1];
1041             y0r = x0r + x2r;
1042             y0i = x0i - x2i;
1043             y2r = x0r - x2r;
1044             y2i = x0i + x2i;
1045             y1r = x1r - x3i;
1046             y1i = x1i - x3r;
1047             y3r = x1r + x3i;
1048             y3i = x1i + x3r;
1049             x0r = a[j4] + a[j5];
1050             x0i = a[j4 + 1] + a[j5 + 1];
1051             x1r = a[j4] - a[j5];
1052             x1i = a[j4 + 1] - a[j5 + 1];
1053             x2r = a[j6] + a[j7];
1054             x2i = a[j6 + 1] + a[j7 + 1];
1055             x3r = a[j6] - a[j7];
1056             x3i = a[j6 + 1] - a[j7 + 1];
1057             y4r = x0r + x2r;
1058             y4i = x0i + x2i;
1059             y6r = x0r - x2r;
1060             y6i = x0i - x2i;
1061             x0r = x1r - x3i;
1062             x0i = x1i + x3r;
1063             x2r = x1r + x3i;
1064             x2i = x1i - x3r;
1065             y5r = wn4r * (x0r - x0i);
1066             y5i = wn4r * (x0r + x0i);
1067             y7r = wn4r * (x2r - x2i);
1068             y7i = wn4r * (x2r + x2i);
1069             a[j1] = y1r + y5r;
1070             a[j1 + 1] = y1i - y5i;
1071             a[j5] = y1r - y5r;
1072             a[j5 + 1] = y1i + y5i;
1073             a[j3] = y3r - y7i;
1074             a[j3 + 1] = y3i - y7r;
1075             a[j7] = y3r + y7i;
1076             a[j7 + 1] = y3i + y7r;
1077             a[j] = y0r + y4r;
1078             a[j + 1] = y0i - y4i;
1079             a[j4] = y0r - y4r;
1080             a[j4 + 1] = y0i + y4i;
1081             a[j2] = y2r - y6i;
1082             a[j2 + 1] = y2i - y6r;
1083             a[j6] = y2r + y6i;
1084             a[j6 + 1] = y2i + y6r;
1085         }
1086     } else if ((l << 2) == n) {
1087         for (j = 0; j < l; j += 2) {
1088             j1 = j + l;
1089             j2 = j1 + l;
1090             j3 = j2 + l;
1091             x0r = a[j] + a[j1];
1092             x0i = -a[j + 1] - a[j1 + 1];
1093             x1r = a[j] - a[j1];
1094             x1i = -a[j + 1] + a[j1 + 1];
1095             x2r = a[j2] + a[j3];
1096             x2i = a[j2 + 1] + a[j3 + 1];
1097             x3r = a[j2] - a[j3];
1098             x3i = a[j2 + 1] - a[j3 + 1];
1099             a[j] = x0r + x2r;
1100             a[j + 1] = x0i - x2i;
1101             a[j2] = x0r - x2r;
1102             a[j2 + 1] = x0i + x2i;
1103             a[j1] = x1r - x3i;
1104             a[j1 + 1] = x1i - x3r;
1105             a[j3] = x1r + x3i;
1106             a[j3 + 1] = x1i + x3r;
1107         }
1108     } else {
1109         for (j = 0; j < l; j += 2) {
1110             j1 = j + l;
1111             x0r = a[j] - a[j1];
1112             x0i = -a[j + 1] + a[j1 + 1];
1113             a[j] += a[j1];
1114             a[j + 1] = -a[j + 1] - a[j1 + 1];
1115             a[j1] = x0r;
1116             a[j1 + 1] = x0i;
1117         }
1118     }
1119 }
1120
1121
1122 void cft1st(int n, smpl_t *a, smpl_t *w)
1123 {
1124     int j, k1;
1125     smpl_t wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
1126         wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i;
1127     smpl_t x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1128         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
1129         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
1130
1131     wn4r = w[2];
1132     x0r = a[0] + a[2];
1133     x0i = a[1] + a[3];
1134     x1r = a[0] - a[2];
1135     x1i = a[1] - a[3];
1136     x2r = a[4] + a[6];
1137     x2i = a[5] + a[7];
1138     x3r = a[4] - a[6];
1139     x3i = a[5] - a[7];
1140     y0r = x0r + x2r;
1141     y0i = x0i + x2i;
1142     y2r = x0r - x2r;
1143     y2i = x0i - x2i;
1144     y1r = x1r - x3i;
1145     y1i = x1i + x3r;
1146     y3r = x1r + x3i;
1147     y3i = x1i - x3r;
1148     x0r = a[8] + a[10];
1149     x0i = a[9] + a[11];
1150     x1r = a[8] - a[10];
1151     x1i = a[9] - a[11];
1152     x2r = a[12] + a[14];
1153     x2i = a[13] + a[15];
1154     x3r = a[12] - a[14];
1155     x3i = a[13] - a[15];
1156     y4r = x0r + x2r;
1157     y4i = x0i + x2i;
1158     y6r = x0r - x2r;
1159     y6i = x0i - x2i;
1160     x0r = x1r - x3i;
1161     x0i = x1i + x3r;
1162     x2r = x1r + x3i;
1163     x2i = x1i - x3r;
1164     y5r = wn4r * (x0r - x0i);
1165     y5i = wn4r * (x0r + x0i);
1166     y7r = wn4r * (x2r - x2i);
1167     y7i = wn4r * (x2r + x2i);
1168     a[2] = y1r + y5r;
1169     a[3] = y1i + y5i;
1170     a[10] = y1r - y5r;
1171     a[11] = y1i - y5i;
1172     a[6] = y3r - y7i;
1173     a[7] = y3i + y7r;
1174     a[14] = y3r + y7i;
1175     a[15] = y3i - y7r;
1176     a[0] = y0r + y4r;
1177     a[1] = y0i + y4i;
1178     a[8] = y0r - y4r;
1179     a[9] = y0i - y4i;
1180     a[4] = y2r - y6i;
1181     a[5] = y2i + y6r;
1182     a[12] = y2r + y6i;
1183     a[13] = y2i - y6r;
1184     if (n > 16) {
1185         wk1r = w[4];
1186         wk1i = w[5];
1187         x0r = a[16] + a[18];
1188         x0i = a[17] + a[19];
1189         x1r = a[16] - a[18];
1190         x1i = a[17] - a[19];
1191         x2r = a[20] + a[22];
1192         x2i = a[21] + a[23];
1193         x3r = a[20] - a[22];
1194         x3i = a[21] - a[23];
1195         y0r = x0r + x2r;
1196         y0i = x0i + x2i;
1197         y2r = x0r - x2r;
1198         y2i = x0i - x2i;
1199         y1r = x1r - x3i;
1200         y1i = x1i + x3r;
1201         y3r = x1r + x3i;
1202         y3i = x1i - x3r;
1203         x0r = a[24] + a[26];
1204         x0i = a[25] + a[27];
1205         x1r = a[24] - a[26];
1206         x1i = a[25] - a[27];
1207         x2r = a[28] + a[30];
1208         x2i = a[29] + a[31];
1209         x3r = a[28] - a[30];
1210         x3i = a[29] - a[31];
1211         y4r = x0r + x2r;
1212         y4i = x0i + x2i;
1213         y6r = x0r - x2r;
1214         y6i = x0i - x2i;
1215         x0r = x1r - x3i;
1216         x0i = x1i + x3r;
1217         x2r = x1r + x3i;
1218         x2i = x3r - x1i;
1219         y5r = wk1i * x0r - wk1r * x0i;
1220         y5i = wk1i * x0i + wk1r * x0r;
1221         y7r = wk1r * x2r + wk1i * x2i;
1222         y7i = wk1r * x2i - wk1i * x2r;
1223         x0r = wk1r * y1r - wk1i * y1i;
1224         x0i = wk1r * y1i + wk1i * y1r;
1225         a[18] = x0r + y5r;
1226         a[19] = x0i + y5i;
1227         a[26] = y5i - x0i;
1228         a[27] = x0r - y5r;
1229         x0r = wk1i * y3r - wk1r * y3i;
1230         x0i = wk1i * y3i + wk1r * y3r;
1231         a[22] = x0r - y7r;
1232         a[23] = x0i + y7i;
1233         a[30] = y7i - x0i;
1234         a[31] = x0r + y7r;
1235         a[16] = y0r + y4r;
1236         a[17] = y0i + y4i;
1237         a[24] = y4i - y0i;
1238         a[25] = y0r - y4r;
1239         x0r = y2r - y6i;
1240         x0i = y2i + y6r;
1241         a[20] = wn4r * (x0r - x0i);
1242         a[21] = wn4r * (x0i + x0r);
1243         x0r = y6r - y2i;
1244         x0i = y2r + y6i;
1245         a[28] = wn4r * (x0r - x0i);
1246         a[29] = wn4r * (x0i + x0r);
1247         k1 = 4;
1248         for (j = 32; j < n; j += 16) {
1249             k1 += 4;
1250             wk1r = w[k1];
1251             wk1i = w[k1 + 1];
1252             wk2r = w[k1 + 2];
1253             wk2i = w[k1 + 3];
1254             wtmp = 2 * wk2i;
1255             wk3r = wk1r - wtmp * wk1i;
1256             wk3i = wtmp * wk1r - wk1i;
1257             wk4r = 1 - wtmp * wk2i;
1258             wk4i = wtmp * wk2r;
1259             wtmp = 2 * wk4i;
1260             wk5r = wk3r - wtmp * wk1i;
1261             wk5i = wtmp * wk1r - wk3i;
1262             wk6r = wk2r - wtmp * wk2i;
1263             wk6i = wtmp * wk2r - wk2i;
1264             wk7r = wk1r - wtmp * wk3i;
1265             wk7i = wtmp * wk3r - wk1i;
1266             x0r = a[j] + a[j + 2];
1267             x0i = a[j + 1] + a[j + 3];
1268             x1r = a[j] - a[j + 2];
1269             x1i = a[j + 1] - a[j + 3];
1270             x2r = a[j + 4] + a[j + 6];
1271             x2i = a[j + 5] + a[j + 7];
1272             x3r = a[j + 4] - a[j + 6];
1273             x3i = a[j + 5] - a[j + 7];
1274             y0r = x0r + x2r;
1275             y0i = x0i + x2i;
1276             y2r = x0r - x2r;
1277             y2i = x0i - x2i;
1278             y1r = x1r - x3i;
1279             y1i = x1i + x3r;
1280             y3r = x1r + x3i;
1281             y3i = x1i - x3r;
1282             x0r = a[j + 8] + a[j + 10];
1283             x0i = a[j + 9] + a[j + 11];
1284             x1r = a[j + 8] - a[j + 10];
1285             x1i = a[j + 9] - a[j + 11];
1286             x2r = a[j + 12] + a[j + 14];
1287             x2i = a[j + 13] + a[j + 15];
1288             x3r = a[j + 12] - a[j + 14];
1289             x3i = a[j + 13] - a[j + 15];
1290             y4r = x0r + x2r;
1291             y4i = x0i + x2i;
1292             y6r = x0r - x2r;
1293             y6i = x0i - x2i;
1294             x0r = x1r - x3i;
1295             x0i = x1i + x3r;
1296             x2r = x1r + x3i;
1297             x2i = x1i - x3r;
1298             y5r = wn4r * (x0r - x0i);
1299             y5i = wn4r * (x0r + x0i);
1300             y7r = wn4r * (x2r - x2i);
1301             y7i = wn4r * (x2r + x2i);
1302             x0r = y1r + y5r;
1303             x0i = y1i + y5i;
1304             a[j + 2] = wk1r * x0r - wk1i * x0i;
1305             a[j + 3] = wk1r * x0i + wk1i * x0r;
1306             x0r = y1r - y5r;
1307             x0i = y1i - y5i;
1308             a[j + 10] = wk5r * x0r - wk5i * x0i;
1309             a[j + 11] = wk5r * x0i + wk5i * x0r;
1310             x0r = y3r - y7i;
1311             x0i = y3i + y7r;
1312             a[j + 6] = wk3r * x0r - wk3i * x0i;
1313             a[j + 7] = wk3r * x0i + wk3i * x0r;
1314             x0r = y3r + y7i;
1315             x0i = y3i - y7r;
1316             a[j + 14] = wk7r * x0r - wk7i * x0i;
1317             a[j + 15] = wk7r * x0i + wk7i * x0r;
1318             a[j] = y0r + y4r;
1319             a[j + 1] = y0i + y4i;
1320             x0r = y0r - y4r;
1321             x0i = y0i - y4i;
1322             a[j + 8] = wk4r * x0r - wk4i * x0i;
1323             a[j + 9] = wk4r * x0i + wk4i * x0r;
1324             x0r = y2r - y6i;
1325             x0i = y2i + y6r;
1326             a[j + 4] = wk2r * x0r - wk2i * x0i;
1327             a[j + 5] = wk2r * x0i + wk2i * x0r;
1328             x0r = y2r + y6i;
1329             x0i = y2i - y6r;
1330             a[j + 12] = wk6r * x0r - wk6i * x0i;
1331             a[j + 13] = wk6r * x0i + wk6i * x0r;
1332         }
1333     }
1334 }
1335
1336
1337 void cftmdl(int n, int l, smpl_t *a, smpl_t *w)
1338 {
1339     int j, j1, j2, j3, j4, j5, j6, j7, k, k1, m;
1340     smpl_t wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
1341         wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i;
1342     smpl_t x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1343         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
1344         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
1345
1346     m = l << 3;
1347     wn4r = w[2];
1348     for (j = 0; j < l; j += 2) {
1349         j1 = j + l;
1350         j2 = j1 + l;
1351         j3 = j2 + l;
1352         j4 = j3 + l;
1353         j5 = j4 + l;
1354         j6 = j5 + l;
1355         j7 = j6 + l;
1356         x0r = a[j] + a[j1];
1357         x0i = a[j + 1] + a[j1 + 1];
1358         x1r = a[j] - a[j1];
1359         x1i = a[j + 1] - a[j1 + 1];
1360         x2r = a[j2] + a[j3];
1361         x2i = a[j2 + 1] + a[j3 + 1];
1362         x3r = a[j2] - a[j3];
1363         x3i = a[j2 + 1] - a[j3 + 1];
1364         y0r = x0r + x2r;
1365         y0i = x0i + x2i;
1366         y2r = x0r - x2r;
1367         y2i = x0i - x2i;
1368         y1r = x1r - x3i;
1369         y1i = x1i + x3r;
1370         y3r = x1r + x3i;
1371         y3i = x1i - x3r;
1372         x0r = a[j4] + a[j5];
1373         x0i = a[j4 + 1] + a[j5 + 1];
1374         x1r = a[j4] - a[j5];
1375         x1i = a[j4 + 1] - a[j5 + 1];
1376         x2r = a[j6] + a[j7];
1377         x2i = a[j6 + 1] + a[j7 + 1];
1378         x3r = a[j6] - a[j7];
1379         x3i = a[j6 + 1] - a[j7 + 1];
1380         y4r = x0r + x2r;
1381         y4i = x0i + x2i;
1382         y6r = x0r - x2r;
1383         y6i = x0i - x2i;
1384         x0r = x1r - x3i;
1385         x0i = x1i + x3r;
1386         x2r = x1r + x3i;
1387         x2i = x1i - x3r;
1388         y5r = wn4r * (x0r - x0i);
1389         y5i = wn4r * (x0r + x0i);
1390         y7r = wn4r * (x2r - x2i);
1391         y7i = wn4r * (x2r + x2i);
1392         a[j1] = y1r + y5r;
1393         a[j1 + 1] = y1i + y5i;
1394         a[j5] = y1r - y5r;
1395         a[j5 + 1] = y1i - y5i;
1396         a[j3] = y3r - y7i;
1397         a[j3 + 1] = y3i + y7r;
1398         a[j7] = y3r + y7i;
1399         a[j7 + 1] = y3i - y7r;
1400         a[j] = y0r + y4r;
1401         a[j + 1] = y0i + y4i;
1402         a[j4] = y0r - y4r;
1403         a[j4 + 1] = y0i - y4i;
1404         a[j2] = y2r - y6i;
1405         a[j2 + 1] = y2i + y6r;
1406         a[j6] = y2r + y6i;
1407         a[j6 + 1] = y2i - y6r;
1408     }
1409     if (m < n) {
1410         wk1r = w[4];
1411         wk1i = w[5];
1412         for (j = m; j < l + m; j += 2) {
1413             j1 = j + l;
1414             j2 = j1 + l;
1415             j3 = j2 + l;
1416             j4 = j3 + l;
1417             j5 = j4 + l;
1418             j6 = j5 + l;
1419             j7 = j6 + l;
1420             x0r = a[j] + a[j1];
1421             x0i = a[j + 1] + a[j1 + 1];
1422             x1r = a[j] - a[j1];
1423             x1i = a[j + 1] - a[j1 + 1];
1424             x2r = a[j2] + a[j3];
1425             x2i = a[j2 + 1] + a[j3 + 1];
1426             x3r = a[j2] - a[j3];
1427             x3i = a[j2 + 1] - a[j3 + 1];
1428             y0r = x0r + x2r;
1429             y0i = x0i + x2i;
1430             y2r = x0r - x2r;
1431             y2i = x0i - x2i;
1432             y1r = x1r - x3i;
1433             y1i = x1i + x3r;
1434             y3r = x1r + x3i;
1435             y3i = x1i - x3r;
1436             x0r = a[j4] + a[j5];
1437             x0i = a[j4 + 1] + a[j5 + 1];
1438             x1r = a[j4] - a[j5];
1439             x1i = a[j4 + 1] - a[j5 + 1];
1440             x2r = a[j6] + a[j7];
1441             x2i = a[j6 + 1] + a[j7 + 1];
1442             x3r = a[j6] - a[j7];
1443             x3i = a[j6 + 1] - a[j7 + 1];
1444             y4r = x0r + x2r;
1445             y4i = x0i + x2i;
1446             y6r = x0r - x2r;
1447             y6i = x0i - x2i;
1448             x0r = x1r - x3i;
1449             x0i = x1i + x3r;
1450             x2r = x1r + x3i;
1451             x2i = x3r - x1i;
1452             y5r = wk1i * x0r - wk1r * x0i;
1453             y5i = wk1i * x0i + wk1r * x0r;
1454             y7r = wk1r * x2r + wk1i * x2i;
1455             y7i = wk1r * x2i - wk1i * x2r;
1456             x0r = wk1r * y1r - wk1i * y1i;
1457             x0i = wk1r * y1i + wk1i * y1r;
1458             a[j1] = x0r + y5r;
1459             a[j1 + 1] = x0i + y5i;
1460             a[j5] = y5i - x0i;
1461             a[j5 + 1] = x0r - y5r;
1462             x0r = wk1i * y3r - wk1r * y3i;
1463             x0i = wk1i * y3i + wk1r * y3r;
1464             a[j3] = x0r - y7r;
1465             a[j3 + 1] = x0i + y7i;
1466             a[j7] = y7i - x0i;
1467             a[j7 + 1] = x0r + y7r;
1468             a[j] = y0r + y4r;
1469             a[j + 1] = y0i + y4i;
1470             a[j4] = y4i - y0i;
1471             a[j4 + 1] = y0r - y4r;
1472             x0r = y2r - y6i;
1473             x0i = y2i + y6r;
1474             a[j2] = wn4r * (x0r - x0i);
1475             a[j2 + 1] = wn4r * (x0i + x0r);
1476             x0r = y6r - y2i;
1477             x0i = y2r + y6i;
1478             a[j6] = wn4r * (x0r - x0i);
1479             a[j6 + 1] = wn4r * (x0i + x0r);
1480         }
1481         k1 = 4;
1482         for (k = 2 * m; k < n; k += m) {
1483             k1 += 4;
1484             wk1r = w[k1];
1485             wk1i = w[k1 + 1];
1486             wk2r = w[k1 + 2];
1487             wk2i = w[k1 + 3];
1488             wtmp = 2 * wk2i;
1489             wk3r = wk1r - wtmp * wk1i;
1490             wk3i = wtmp * wk1r - wk1i;
1491             wk4r = 1 - wtmp * wk2i;
1492             wk4i = wtmp * wk2r;
1493             wtmp = 2 * wk4i;
1494             wk5r = wk3r - wtmp * wk1i;
1495             wk5i = wtmp * wk1r - wk3i;
1496             wk6r = wk2r - wtmp * wk2i;
1497             wk6i = wtmp * wk2r - wk2i;
1498             wk7r = wk1r - wtmp * wk3i;
1499             wk7i = wtmp * wk3r - wk1i;
1500             for (j = k; j < l + k; j += 2) {
1501                 j1 = j + l;
1502                 j2 = j1 + l;
1503                 j3 = j2 + l;
1504                 j4 = j3 + l;
1505                 j5 = j4 + l;
1506                 j6 = j5 + l;
1507                 j7 = j6 + l;
1508                 x0r = a[j] + a[j1];
1509                 x0i = a[j + 1] + a[j1 + 1];
1510                 x1r = a[j] - a[j1];
1511                 x1i = a[j + 1] - a[j1 + 1];
1512                 x2r = a[j2] + a[j3];
1513                 x2i = a[j2 + 1] + a[j3 + 1];
1514                 x3r = a[j2] - a[j3];
1515                 x3i = a[j2 + 1] - a[j3 + 1];
1516                 y0r = x0r + x2r;
1517                 y0i = x0i + x2i;
1518                 y2r = x0r - x2r;
1519                 y2i = x0i - x2i;
1520                 y1r = x1r - x3i;
1521                 y1i = x1i + x3r;
1522                 y3r = x1r + x3i;
1523                 y3i = x1i - x3r;
1524                 x0r = a[j4] + a[j5];
1525                 x0i = a[j4 + 1] + a[j5 + 1];
1526                 x1r = a[j4] - a[j5];
1527                 x1i = a[j4 + 1] - a[j5 + 1];
1528                 x2r = a[j6] + a[j7];
1529                 x2i = a[j6 + 1] + a[j7 + 1];
1530                 x3r = a[j6] - a[j7];
1531                 x3i = a[j6 + 1] - a[j7 + 1];
1532                 y4r = x0r + x2r;
1533                 y4i = x0i + x2i;
1534                 y6r = x0r - x2r;
1535                 y6i = x0i - x2i;
1536                 x0r = x1r - x3i;
1537                 x0i = x1i + x3r;
1538                 x2r = x1r + x3i;
1539                 x2i = x1i - x3r;
1540                 y5r = wn4r * (x0r - x0i);
1541                 y5i = wn4r * (x0r + x0i);
1542                 y7r = wn4r * (x2r - x2i);
1543                 y7i = wn4r * (x2r + x2i);
1544                 x0r = y1r + y5r;
1545                 x0i = y1i + y5i;
1546                 a[j1] = wk1r * x0r - wk1i * x0i;
1547                 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
1548                 x0r = y1r - y5r;
1549                 x0i = y1i - y5i;
1550                 a[j5] = wk5r * x0r - wk5i * x0i;
1551                 a[j5 + 1] = wk5r * x0i + wk5i * x0r;
1552                 x0r = y3r - y7i;
1553                 x0i = y3i + y7r;
1554                 a[j3] = wk3r * x0r - wk3i * x0i;
1555                 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
1556                 x0r = y3r + y7i;
1557                 x0i = y3i - y7r;
1558                 a[j7] = wk7r * x0r - wk7i * x0i;
1559                 a[j7 + 1] = wk7r * x0i + wk7i * x0r;
1560                 a[j] = y0r + y4r;
1561                 a[j + 1] = y0i + y4i;
1562                 x0r = y0r - y4r;
1563                 x0i = y0i - y4i;
1564                 a[j4] = wk4r * x0r - wk4i * x0i;
1565                 a[j4 + 1] = wk4r * x0i + wk4i * x0r;
1566                 x0r = y2r - y6i;
1567                 x0i = y2i + y6r;
1568                 a[j2] = wk2r * x0r - wk2i * x0i;
1569                 a[j2 + 1] = wk2r * x0i + wk2i * x0r;
1570                 x0r = y2r + y6i;
1571                 x0i = y2i - y6r;
1572                 a[j6] = wk6r * x0r - wk6i * x0i;
1573                 a[j6 + 1] = wk6r * x0i + wk6i * x0r;
1574             }
1575         }
1576     }
1577 }
1578
1579
1580 void rftfsub(int n, smpl_t *a, int nc, smpl_t *c)
1581 {
1582     int j, k, kk, ks, m;
1583     smpl_t wkr, wki, xr, xi, yr, yi;
1584
1585     m = n >> 1;
1586     ks = 2 * nc / m;
1587     kk = 0;
1588     for (j = 2; j < m; j += 2) {
1589         k = n - j;
1590         kk += ks;
1591         wkr = (smpl_t)0.5 - c[nc - kk];
1592         wki = c[kk];
1593         xr = a[j] - a[k];
1594         xi = a[j + 1] + a[k + 1];
1595         yr = wkr * xr - wki * xi;
1596         yi = wkr * xi + wki * xr;
1597         a[j] -= yr;
1598         a[j + 1] -= yi;
1599         a[k] += yr;
1600         a[k + 1] -= yi;
1601     }
1602 }
1603
1604
1605 void rftbsub(int n, smpl_t *a, int nc, smpl_t *c)
1606 {
1607     int j, k, kk, ks, m;
1608     smpl_t wkr, wki, xr, xi, yr, yi;
1609
1610     a[1] = -a[1];
1611     m = n >> 1;
1612     ks = 2 * nc / m;
1613     kk = 0;
1614     for (j = 2; j < m; j += 2) {
1615         k = n - j;
1616         kk += ks;
1617         wkr = (smpl_t)0.5 - c[nc - kk];
1618         wki = c[kk];
1619         xr = a[j] - a[k];
1620         xi = a[j + 1] + a[k + 1];
1621         yr = wkr * xr + wki * xi;
1622         yi = wkr * xi - wki * xr;
1623         a[j] -= yr;
1624         a[j + 1] = yi - a[j + 1];
1625         a[k] += yr;
1626         a[k + 1] = yi - a[k + 1];
1627     }
1628     a[m + 1] = -a[m + 1];
1629 }
1630
1631
1632 void dctsub(int n, smpl_t *a, int nc, smpl_t *c)
1633 {
1634     int j, k, kk, ks, m;
1635     smpl_t wkr, wki, xr;
1636
1637     m = n >> 1;
1638     ks = nc / n;
1639     kk = 0;
1640     for (j = 1; j < m; j++) {
1641         k = n - j;
1642         kk += ks;
1643         wkr = c[kk] - c[nc - kk];
1644         wki = c[kk] + c[nc - kk];
1645         xr = wki * a[j] - wkr * a[k];
1646         a[j] = wkr * a[j] + wki * a[k];
1647         a[k] = xr;
1648     }
1649     a[m] *= c[0];
1650 }
1651
1652
1653 void dstsub(int n, smpl_t *a, int nc, smpl_t *c)
1654 {
1655     int j, k, kk, ks, m;
1656     smpl_t wkr, wki, xr;
1657
1658     m = n >> 1;
1659     ks = nc / n;
1660     kk = 0;
1661     for (j = 1; j < m; j++) {
1662         k = n - j;
1663         kk += ks;
1664         wkr = c[kk] - c[nc - kk];
1665         wki = c[kk] + c[nc - kk];
1666         xr = wki * a[k] - wkr * a[j];
1667         a[k] = wkr * a[k] + wki * a[j];
1668         a[j] = xr;
1669     }
1670     a[m] *= c[0];
1671 }
1672