dbf-halloween2015

annotate libs/vorbis/mdct.c @ 1:c3f5c32cb210

find changesets by author, revision, files, or words in the commit message
barfed all the libraries in the source tree to make porting easier
author John Tsiombikas <nuclear@member.fsf.org>
date Sun, 01 Nov 2015 00:36:56 +0200
parents
children
rev   line source
nuclear@1 1 /********************************************************************
nuclear@1 2 * *
nuclear@1 3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
nuclear@1 4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
nuclear@1 5 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
nuclear@1 6 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
nuclear@1 7 * *
nuclear@1 8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 *
nuclear@1 9 * by the Xiph.Org Foundation http://www.xiph.org/ *
nuclear@1 10 * *
nuclear@1 11 ********************************************************************
nuclear@1 12
nuclear@1 13 function: normalized modified discrete cosine transform
nuclear@1 14 power of two length transform only [64 <= n ]
nuclear@1 15 last mod: $Id: mdct.c 16227 2009-07-08 06:58:46Z xiphmont $
nuclear@1 16
nuclear@1 17 Original algorithm adapted long ago from _The use of multirate filter
nuclear@1 18 banks for coding of high quality digital audio_, by T. Sporer,
nuclear@1 19 K. Brandenburg and B. Edler, collection of the European Signal
nuclear@1 20 Processing Conference (EUSIPCO), Amsterdam, June 1992, Vol.1, pp
nuclear@1 21 211-214
nuclear@1 22
nuclear@1 23 The below code implements an algorithm that no longer looks much like
nuclear@1 24 that presented in the paper, but the basic structure remains if you
nuclear@1 25 dig deep enough to see it.
nuclear@1 26
nuclear@1 27 This module DOES NOT INCLUDE code to generate/apply the window
nuclear@1 28 function. Everybody has their own weird favorite including me... I
nuclear@1 29 happen to like the properties of y=sin(.5PI*sin^2(x)), but others may
nuclear@1 30 vehemently disagree.
nuclear@1 31
nuclear@1 32 ********************************************************************/
nuclear@1 33
nuclear@1 34 /* this can also be run as an integer transform by uncommenting a
nuclear@1 35 define in mdct.h; the integerization is a first pass and although
nuclear@1 36 it's likely stable for Vorbis, the dynamic range is constrained and
nuclear@1 37 roundoff isn't done (so it's noisy). Consider it functional, but
nuclear@1 38 only a starting point. There's no point on a machine with an FPU */
nuclear@1 39
nuclear@1 40 #include <stdio.h>
nuclear@1 41 #include <stdlib.h>
nuclear@1 42 #include <string.h>
nuclear@1 43 #include <math.h>
nuclear@1 44 #include "vorbis/codec.h"
nuclear@1 45 #include "mdct.h"
nuclear@1 46 #include "os.h"
nuclear@1 47 #include "misc.h"
nuclear@1 48
nuclear@1 49 /* build lookups for trig functions; also pre-figure scaling and
nuclear@1 50 some window function algebra. */
nuclear@1 51
nuclear@1 52 void mdct_init(mdct_lookup *lookup,int n){
nuclear@1 53 int *bitrev=_ogg_malloc(sizeof(*bitrev)*(n/4));
nuclear@1 54 DATA_TYPE *T=_ogg_malloc(sizeof(*T)*(n+n/4));
nuclear@1 55
nuclear@1 56 int i;
nuclear@1 57 int n2=n>>1;
nuclear@1 58 int log2n=lookup->log2n=rint(log((float)n)/log(2.f));
nuclear@1 59 lookup->n=n;
nuclear@1 60 lookup->trig=T;
nuclear@1 61 lookup->bitrev=bitrev;
nuclear@1 62
nuclear@1 63 /* trig lookups... */
nuclear@1 64
nuclear@1 65 for(i=0;i<n/4;i++){
nuclear@1 66 T[i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i)));
nuclear@1 67 T[i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i)));
nuclear@1 68 T[n2+i*2]=FLOAT_CONV(cos((M_PI/(2*n))*(2*i+1)));
nuclear@1 69 T[n2+i*2+1]=FLOAT_CONV(sin((M_PI/(2*n))*(2*i+1)));
nuclear@1 70 }
nuclear@1 71 for(i=0;i<n/8;i++){
nuclear@1 72 T[n+i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i+2))*.5);
nuclear@1 73 T[n+i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i+2))*.5);
nuclear@1 74 }
nuclear@1 75
nuclear@1 76 /* bitreverse lookup... */
nuclear@1 77
nuclear@1 78 {
nuclear@1 79 int mask=(1<<(log2n-1))-1,i,j;
nuclear@1 80 int msb=1<<(log2n-2);
nuclear@1 81 for(i=0;i<n/8;i++){
nuclear@1 82 int acc=0;
nuclear@1 83 for(j=0;msb>>j;j++)
nuclear@1 84 if((msb>>j)&i)acc|=1<<j;
nuclear@1 85 bitrev[i*2]=((~acc)&mask)-1;
nuclear@1 86 bitrev[i*2+1]=acc;
nuclear@1 87
nuclear@1 88 }
nuclear@1 89 }
nuclear@1 90 lookup->scale=FLOAT_CONV(4.f/n);
nuclear@1 91 }
nuclear@1 92
nuclear@1 93 /* 8 point butterfly (in place, 4 register) */
nuclear@1 94 STIN void mdct_butterfly_8(DATA_TYPE *x){
nuclear@1 95 REG_TYPE r0 = x[6] + x[2];
nuclear@1 96 REG_TYPE r1 = x[6] - x[2];
nuclear@1 97 REG_TYPE r2 = x[4] + x[0];
nuclear@1 98 REG_TYPE r3 = x[4] - x[0];
nuclear@1 99
nuclear@1 100 x[6] = r0 + r2;
nuclear@1 101 x[4] = r0 - r2;
nuclear@1 102
nuclear@1 103 r0 = x[5] - x[1];
nuclear@1 104 r2 = x[7] - x[3];
nuclear@1 105 x[0] = r1 + r0;
nuclear@1 106 x[2] = r1 - r0;
nuclear@1 107
nuclear@1 108 r0 = x[5] + x[1];
nuclear@1 109 r1 = x[7] + x[3];
nuclear@1 110 x[3] = r2 + r3;
nuclear@1 111 x[1] = r2 - r3;
nuclear@1 112 x[7] = r1 + r0;
nuclear@1 113 x[5] = r1 - r0;
nuclear@1 114
nuclear@1 115 }
nuclear@1 116
nuclear@1 117 /* 16 point butterfly (in place, 4 register) */
nuclear@1 118 STIN void mdct_butterfly_16(DATA_TYPE *x){
nuclear@1 119 REG_TYPE r0 = x[1] - x[9];
nuclear@1 120 REG_TYPE r1 = x[0] - x[8];
nuclear@1 121
nuclear@1 122 x[8] += x[0];
nuclear@1 123 x[9] += x[1];
nuclear@1 124 x[0] = MULT_NORM((r0 + r1) * cPI2_8);
nuclear@1 125 x[1] = MULT_NORM((r0 - r1) * cPI2_8);
nuclear@1 126
nuclear@1 127 r0 = x[3] - x[11];
nuclear@1 128 r1 = x[10] - x[2];
nuclear@1 129 x[10] += x[2];
nuclear@1 130 x[11] += x[3];
nuclear@1 131 x[2] = r0;
nuclear@1 132 x[3] = r1;
nuclear@1 133
nuclear@1 134 r0 = x[12] - x[4];
nuclear@1 135 r1 = x[13] - x[5];
nuclear@1 136 x[12] += x[4];
nuclear@1 137 x[13] += x[5];
nuclear@1 138 x[4] = MULT_NORM((r0 - r1) * cPI2_8);
nuclear@1 139 x[5] = MULT_NORM((r0 + r1) * cPI2_8);
nuclear@1 140
nuclear@1 141 r0 = x[14] - x[6];
nuclear@1 142 r1 = x[15] - x[7];
nuclear@1 143 x[14] += x[6];
nuclear@1 144 x[15] += x[7];
nuclear@1 145 x[6] = r0;
nuclear@1 146 x[7] = r1;
nuclear@1 147
nuclear@1 148 mdct_butterfly_8(x);
nuclear@1 149 mdct_butterfly_8(x+8);
nuclear@1 150 }
nuclear@1 151
nuclear@1 152 /* 32 point butterfly (in place, 4 register) */
nuclear@1 153 STIN void mdct_butterfly_32(DATA_TYPE *x){
nuclear@1 154 REG_TYPE r0 = x[30] - x[14];
nuclear@1 155 REG_TYPE r1 = x[31] - x[15];
nuclear@1 156
nuclear@1 157 x[30] += x[14];
nuclear@1 158 x[31] += x[15];
nuclear@1 159 x[14] = r0;
nuclear@1 160 x[15] = r1;
nuclear@1 161
nuclear@1 162 r0 = x[28] - x[12];
nuclear@1 163 r1 = x[29] - x[13];
nuclear@1 164 x[28] += x[12];
nuclear@1 165 x[29] += x[13];
nuclear@1 166 x[12] = MULT_NORM( r0 * cPI1_8 - r1 * cPI3_8 );
nuclear@1 167 x[13] = MULT_NORM( r0 * cPI3_8 + r1 * cPI1_8 );
nuclear@1 168
nuclear@1 169 r0 = x[26] - x[10];
nuclear@1 170 r1 = x[27] - x[11];
nuclear@1 171 x[26] += x[10];
nuclear@1 172 x[27] += x[11];
nuclear@1 173 x[10] = MULT_NORM(( r0 - r1 ) * cPI2_8);
nuclear@1 174 x[11] = MULT_NORM(( r0 + r1 ) * cPI2_8);
nuclear@1 175
nuclear@1 176 r0 = x[24] - x[8];
nuclear@1 177 r1 = x[25] - x[9];
nuclear@1 178 x[24] += x[8];
nuclear@1 179 x[25] += x[9];
nuclear@1 180 x[8] = MULT_NORM( r0 * cPI3_8 - r1 * cPI1_8 );
nuclear@1 181 x[9] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 );
nuclear@1 182
nuclear@1 183 r0 = x[22] - x[6];
nuclear@1 184 r1 = x[7] - x[23];
nuclear@1 185 x[22] += x[6];
nuclear@1 186 x[23] += x[7];
nuclear@1 187 x[6] = r1;
nuclear@1 188 x[7] = r0;
nuclear@1 189
nuclear@1 190 r0 = x[4] - x[20];
nuclear@1 191 r1 = x[5] - x[21];
nuclear@1 192 x[20] += x[4];
nuclear@1 193 x[21] += x[5];
nuclear@1 194 x[4] = MULT_NORM( r1 * cPI1_8 + r0 * cPI3_8 );
nuclear@1 195 x[5] = MULT_NORM( r1 * cPI3_8 - r0 * cPI1_8 );
nuclear@1 196
nuclear@1 197 r0 = x[2] - x[18];
nuclear@1 198 r1 = x[3] - x[19];
nuclear@1 199 x[18] += x[2];
nuclear@1 200 x[19] += x[3];
nuclear@1 201 x[2] = MULT_NORM(( r1 + r0 ) * cPI2_8);
nuclear@1 202 x[3] = MULT_NORM(( r1 - r0 ) * cPI2_8);
nuclear@1 203
nuclear@1 204 r0 = x[0] - x[16];
nuclear@1 205 r1 = x[1] - x[17];
nuclear@1 206 x[16] += x[0];
nuclear@1 207 x[17] += x[1];
nuclear@1 208 x[0] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 );
nuclear@1 209 x[1] = MULT_NORM( r1 * cPI1_8 - r0 * cPI3_8 );
nuclear@1 210
nuclear@1 211 mdct_butterfly_16(x);
nuclear@1 212 mdct_butterfly_16(x+16);
nuclear@1 213
nuclear@1 214 }
nuclear@1 215
nuclear@1 216 /* N point first stage butterfly (in place, 2 register) */
nuclear@1 217 STIN void mdct_butterfly_first(DATA_TYPE *T,
nuclear@1 218 DATA_TYPE *x,
nuclear@1 219 int points){
nuclear@1 220
nuclear@1 221 DATA_TYPE *x1 = x + points - 8;
nuclear@1 222 DATA_TYPE *x2 = x + (points>>1) - 8;
nuclear@1 223 REG_TYPE r0;
nuclear@1 224 REG_TYPE r1;
nuclear@1 225
nuclear@1 226 do{
nuclear@1 227
nuclear@1 228 r0 = x1[6] - x2[6];
nuclear@1 229 r1 = x1[7] - x2[7];
nuclear@1 230 x1[6] += x2[6];
nuclear@1 231 x1[7] += x2[7];
nuclear@1 232 x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]);
nuclear@1 233 x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]);
nuclear@1 234
nuclear@1 235 r0 = x1[4] - x2[4];
nuclear@1 236 r1 = x1[5] - x2[5];
nuclear@1 237 x1[4] += x2[4];
nuclear@1 238 x1[5] += x2[5];
nuclear@1 239 x2[4] = MULT_NORM(r1 * T[5] + r0 * T[4]);
nuclear@1 240 x2[5] = MULT_NORM(r1 * T[4] - r0 * T[5]);
nuclear@1 241
nuclear@1 242 r0 = x1[2] - x2[2];
nuclear@1 243 r1 = x1[3] - x2[3];
nuclear@1 244 x1[2] += x2[2];
nuclear@1 245 x1[3] += x2[3];
nuclear@1 246 x2[2] = MULT_NORM(r1 * T[9] + r0 * T[8]);
nuclear@1 247 x2[3] = MULT_NORM(r1 * T[8] - r0 * T[9]);
nuclear@1 248
nuclear@1 249 r0 = x1[0] - x2[0];
nuclear@1 250 r1 = x1[1] - x2[1];
nuclear@1 251 x1[0] += x2[0];
nuclear@1 252 x1[1] += x2[1];
nuclear@1 253 x2[0] = MULT_NORM(r1 * T[13] + r0 * T[12]);
nuclear@1 254 x2[1] = MULT_NORM(r1 * T[12] - r0 * T[13]);
nuclear@1 255
nuclear@1 256 x1-=8;
nuclear@1 257 x2-=8;
nuclear@1 258 T+=16;
nuclear@1 259
nuclear@1 260 }while(x2>=x);
nuclear@1 261 }
nuclear@1 262
nuclear@1 263 /* N/stage point generic N stage butterfly (in place, 2 register) */
nuclear@1 264 STIN void mdct_butterfly_generic(DATA_TYPE *T,
nuclear@1 265 DATA_TYPE *x,
nuclear@1 266 int points,
nuclear@1 267 int trigint){
nuclear@1 268
nuclear@1 269 DATA_TYPE *x1 = x + points - 8;
nuclear@1 270 DATA_TYPE *x2 = x + (points>>1) - 8;
nuclear@1 271 REG_TYPE r0;
nuclear@1 272 REG_TYPE r1;
nuclear@1 273
nuclear@1 274 do{
nuclear@1 275
nuclear@1 276 r0 = x1[6] - x2[6];
nuclear@1 277 r1 = x1[7] - x2[7];
nuclear@1 278 x1[6] += x2[6];
nuclear@1 279 x1[7] += x2[7];
nuclear@1 280 x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]);
nuclear@1 281 x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]);
nuclear@1 282
nuclear@1 283 T+=trigint;
nuclear@1 284
nuclear@1 285 r0 = x1[4] - x2[4];
nuclear@1 286 r1 = x1[5] - x2[5];
nuclear@1 287 x1[4] += x2[4];
nuclear@1 288 x1[5] += x2[5];
nuclear@1 289 x2[4] = MULT_NORM(r1 * T[1] + r0 * T[0]);
nuclear@1 290 x2[5] = MULT_NORM(r1 * T[0] - r0 * T[1]);
nuclear@1 291
nuclear@1 292 T+=trigint;
nuclear@1 293
nuclear@1 294 r0 = x1[2] - x2[2];
nuclear@1 295 r1 = x1[3] - x2[3];
nuclear@1 296 x1[2] += x2[2];
nuclear@1 297 x1[3] += x2[3];
nuclear@1 298 x2[2] = MULT_NORM(r1 * T[1] + r0 * T[0]);
nuclear@1 299 x2[3] = MULT_NORM(r1 * T[0] - r0 * T[1]);
nuclear@1 300
nuclear@1 301 T+=trigint;
nuclear@1 302
nuclear@1 303 r0 = x1[0] - x2[0];
nuclear@1 304 r1 = x1[1] - x2[1];
nuclear@1 305 x1[0] += x2[0];
nuclear@1 306 x1[1] += x2[1];
nuclear@1 307 x2[0] = MULT_NORM(r1 * T[1] + r0 * T[0]);
nuclear@1 308 x2[1] = MULT_NORM(r1 * T[0] - r0 * T[1]);
nuclear@1 309
nuclear@1 310 T+=trigint;
nuclear@1 311 x1-=8;
nuclear@1 312 x2-=8;
nuclear@1 313
nuclear@1 314 }while(x2>=x);
nuclear@1 315 }
nuclear@1 316
nuclear@1 317 STIN void mdct_butterflies(mdct_lookup *init,
nuclear@1 318 DATA_TYPE *x,
nuclear@1 319 int points){
nuclear@1 320
nuclear@1 321 DATA_TYPE *T=init->trig;
nuclear@1 322 int stages=init->log2n-5;
nuclear@1 323 int i,j;
nuclear@1 324
nuclear@1 325 if(--stages>0){
nuclear@1 326 mdct_butterfly_first(T,x,points);
nuclear@1 327 }
nuclear@1 328
nuclear@1 329 for(i=1;--stages>0;i++){
nuclear@1 330 for(j=0;j<(1<<i);j++)
nuclear@1 331 mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i);
nuclear@1 332 }
nuclear@1 333
nuclear@1 334 for(j=0;j<points;j+=32)
nuclear@1 335 mdct_butterfly_32(x+j);
nuclear@1 336
nuclear@1 337 }
nuclear@1 338
nuclear@1 339 void mdct_clear(mdct_lookup *l){
nuclear@1 340 if(l){
nuclear@1 341 if(l->trig)_ogg_free(l->trig);
nuclear@1 342 if(l->bitrev)_ogg_free(l->bitrev);
nuclear@1 343 memset(l,0,sizeof(*l));
nuclear@1 344 }
nuclear@1 345 }
nuclear@1 346
nuclear@1 347 STIN void mdct_bitreverse(mdct_lookup *init,
nuclear@1 348 DATA_TYPE *x){
nuclear@1 349 int n = init->n;
nuclear@1 350 int *bit = init->bitrev;
nuclear@1 351 DATA_TYPE *w0 = x;
nuclear@1 352 DATA_TYPE *w1 = x = w0+(n>>1);
nuclear@1 353 DATA_TYPE *T = init->trig+n;
nuclear@1 354
nuclear@1 355 do{
nuclear@1 356 DATA_TYPE *x0 = x+bit[0];
nuclear@1 357 DATA_TYPE *x1 = x+bit[1];
nuclear@1 358
nuclear@1 359 REG_TYPE r0 = x0[1] - x1[1];
nuclear@1 360 REG_TYPE r1 = x0[0] + x1[0];
nuclear@1 361 REG_TYPE r2 = MULT_NORM(r1 * T[0] + r0 * T[1]);
nuclear@1 362 REG_TYPE r3 = MULT_NORM(r1 * T[1] - r0 * T[0]);
nuclear@1 363
nuclear@1 364 w1 -= 4;
nuclear@1 365
nuclear@1 366 r0 = HALVE(x0[1] + x1[1]);
nuclear@1 367 r1 = HALVE(x0[0] - x1[0]);
nuclear@1 368
nuclear@1 369 w0[0] = r0 + r2;
nuclear@1 370 w1[2] = r0 - r2;
nuclear@1 371 w0[1] = r1 + r3;
nuclear@1 372 w1[3] = r3 - r1;
nuclear@1 373
nuclear@1 374 x0 = x+bit[2];
nuclear@1 375 x1 = x+bit[3];
nuclear@1 376
nuclear@1 377 r0 = x0[1] - x1[1];
nuclear@1 378 r1 = x0[0] + x1[0];
nuclear@1 379 r2 = MULT_NORM(r1 * T[2] + r0 * T[3]);
nuclear@1 380 r3 = MULT_NORM(r1 * T[3] - r0 * T[2]);
nuclear@1 381
nuclear@1 382 r0 = HALVE(x0[1] + x1[1]);
nuclear@1 383 r1 = HALVE(x0[0] - x1[0]);
nuclear@1 384
nuclear@1 385 w0[2] = r0 + r2;
nuclear@1 386 w1[0] = r0 - r2;
nuclear@1 387 w0[3] = r1 + r3;
nuclear@1 388 w1[1] = r3 - r1;
nuclear@1 389
nuclear@1 390 T += 4;
nuclear@1 391 bit += 4;
nuclear@1 392 w0 += 4;
nuclear@1 393
nuclear@1 394 }while(w0<w1);
nuclear@1 395 }
nuclear@1 396
nuclear@1 397 void mdct_backward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
nuclear@1 398 int n=init->n;
nuclear@1 399 int n2=n>>1;
nuclear@1 400 int n4=n>>2;
nuclear@1 401
nuclear@1 402 /* rotate */
nuclear@1 403
nuclear@1 404 DATA_TYPE *iX = in+n2-7;
nuclear@1 405 DATA_TYPE *oX = out+n2+n4;
nuclear@1 406 DATA_TYPE *T = init->trig+n4;
nuclear@1 407
nuclear@1 408 do{
nuclear@1 409 oX -= 4;
nuclear@1 410 oX[0] = MULT_NORM(-iX[2] * T[3] - iX[0] * T[2]);
nuclear@1 411 oX[1] = MULT_NORM (iX[0] * T[3] - iX[2] * T[2]);
nuclear@1 412 oX[2] = MULT_NORM(-iX[6] * T[1] - iX[4] * T[0]);
nuclear@1 413 oX[3] = MULT_NORM (iX[4] * T[1] - iX[6] * T[0]);
nuclear@1 414 iX -= 8;
nuclear@1 415 T += 4;
nuclear@1 416 }while(iX>=in);
nuclear@1 417
nuclear@1 418 iX = in+n2-8;
nuclear@1 419 oX = out+n2+n4;
nuclear@1 420 T = init->trig+n4;
nuclear@1 421
nuclear@1 422 do{
nuclear@1 423 T -= 4;
nuclear@1 424 oX[0] = MULT_NORM (iX[4] * T[3] + iX[6] * T[2]);
nuclear@1 425 oX[1] = MULT_NORM (iX[4] * T[2] - iX[6] * T[3]);
nuclear@1 426 oX[2] = MULT_NORM (iX[0] * T[1] + iX[2] * T[0]);
nuclear@1 427 oX[3] = MULT_NORM (iX[0] * T[0] - iX[2] * T[1]);
nuclear@1 428 iX -= 8;
nuclear@1 429 oX += 4;
nuclear@1 430 }while(iX>=in);
nuclear@1 431
nuclear@1 432 mdct_butterflies(init,out+n2,n2);
nuclear@1 433 mdct_bitreverse(init,out);
nuclear@1 434
nuclear@1 435 /* roatate + window */
nuclear@1 436
nuclear@1 437 {
nuclear@1 438 DATA_TYPE *oX1=out+n2+n4;
nuclear@1 439 DATA_TYPE *oX2=out+n2+n4;
nuclear@1 440 DATA_TYPE *iX =out;
nuclear@1 441 T =init->trig+n2;
nuclear@1 442
nuclear@1 443 do{
nuclear@1 444 oX1-=4;
nuclear@1 445
nuclear@1 446 oX1[3] = MULT_NORM (iX[0] * T[1] - iX[1] * T[0]);
nuclear@1 447 oX2[0] = -MULT_NORM (iX[0] * T[0] + iX[1] * T[1]);
nuclear@1 448
nuclear@1 449 oX1[2] = MULT_NORM (iX[2] * T[3] - iX[3] * T[2]);
nuclear@1 450 oX2[1] = -MULT_NORM (iX[2] * T[2] + iX[3] * T[3]);
nuclear@1 451
nuclear@1 452 oX1[1] = MULT_NORM (iX[4] * T[5] - iX[5] * T[4]);
nuclear@1 453 oX2[2] = -MULT_NORM (iX[4] * T[4] + iX[5] * T[5]);
nuclear@1 454
nuclear@1 455 oX1[0] = MULT_NORM (iX[6] * T[7] - iX[7] * T[6]);
nuclear@1 456 oX2[3] = -MULT_NORM (iX[6] * T[6] + iX[7] * T[7]);
nuclear@1 457
nuclear@1 458 oX2+=4;
nuclear@1 459 iX += 8;
nuclear@1 460 T += 8;
nuclear@1 461 }while(iX<oX1);
nuclear@1 462
nuclear@1 463 iX=out+n2+n4;
nuclear@1 464 oX1=out+n4;
nuclear@1 465 oX2=oX1;
nuclear@1 466
nuclear@1 467 do{
nuclear@1 468 oX1-=4;
nuclear@1 469 iX-=4;
nuclear@1 470
nuclear@1 471 oX2[0] = -(oX1[3] = iX[3]);
nuclear@1 472 oX2[1] = -(oX1[2] = iX[2]);
nuclear@1 473 oX2[2] = -(oX1[1] = iX[1]);
nuclear@1 474 oX2[3] = -(oX1[0] = iX[0]);
nuclear@1 475
nuclear@1 476 oX2+=4;
nuclear@1 477 }while(oX2<iX);
nuclear@1 478
nuclear@1 479 iX=out+n2+n4;
nuclear@1 480 oX1=out+n2+n4;
nuclear@1 481 oX2=out+n2;
nuclear@1 482 do{
nuclear@1 483 oX1-=4;
nuclear@1 484 oX1[0]= iX[3];
nuclear@1 485 oX1[1]= iX[2];
nuclear@1 486 oX1[2]= iX[1];
nuclear@1 487 oX1[3]= iX[0];
nuclear@1 488 iX+=4;
nuclear@1 489 }while(oX1>oX2);
nuclear@1 490 }
nuclear@1 491 }
nuclear@1 492
nuclear@1 493 void mdct_forward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
nuclear@1 494 int n=init->n;
nuclear@1 495 int n2=n>>1;
nuclear@1 496 int n4=n>>2;
nuclear@1 497 int n8=n>>3;
nuclear@1 498 DATA_TYPE *w=alloca(n*sizeof(*w)); /* forward needs working space */
nuclear@1 499 DATA_TYPE *w2=w+n2;
nuclear@1 500
nuclear@1 501 /* rotate */
nuclear@1 502
nuclear@1 503 /* window + rotate + step 1 */
nuclear@1 504
nuclear@1 505 REG_TYPE r0;
nuclear@1 506 REG_TYPE r1;
nuclear@1 507 DATA_TYPE *x0=in+n2+n4;
nuclear@1 508 DATA_TYPE *x1=x0+1;
nuclear@1 509 DATA_TYPE *T=init->trig+n2;
nuclear@1 510
nuclear@1 511 int i=0;
nuclear@1 512
nuclear@1 513 for(i=0;i<n8;i+=2){
nuclear@1 514 x0 -=4;
nuclear@1 515 T-=2;
nuclear@1 516 r0= x0[2] + x1[0];
nuclear@1 517 r1= x0[0] + x1[2];
nuclear@1 518 w2[i]= MULT_NORM(r1*T[1] + r0*T[0]);
nuclear@1 519 w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
nuclear@1 520 x1 +=4;
nuclear@1 521 }
nuclear@1 522
nuclear@1 523 x1=in+1;
nuclear@1 524
nuclear@1 525 for(;i<n2-n8;i+=2){
nuclear@1 526 T-=2;
nuclear@1 527 x0 -=4;
nuclear@1 528 r0= x0[2] - x1[0];
nuclear@1 529 r1= x0[0] - x1[2];
nuclear@1 530 w2[i]= MULT_NORM(r1*T[1] + r0*T[0]);
nuclear@1 531 w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
nuclear@1 532 x1 +=4;
nuclear@1 533 }
nuclear@1 534
nuclear@1 535 x0=in+n;
nuclear@1 536
nuclear@1 537 for(;i<n2;i+=2){
nuclear@1 538 T-=2;
nuclear@1 539 x0 -=4;
nuclear@1 540 r0= -x0[2] - x1[0];
nuclear@1 541 r1= -x0[0] - x1[2];
nuclear@1 542 w2[i]= MULT_NORM(r1*T[1] + r0*T[0]);
nuclear@1 543 w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
nuclear@1 544 x1 +=4;
nuclear@1 545 }
nuclear@1 546
nuclear@1 547
nuclear@1 548 mdct_butterflies(init,w+n2,n2);
nuclear@1 549 mdct_bitreverse(init,w);
nuclear@1 550
nuclear@1 551 /* roatate + window */
nuclear@1 552
nuclear@1 553 T=init->trig+n2;
nuclear@1 554 x0=out+n2;
nuclear@1 555
nuclear@1 556 for(i=0;i<n4;i++){
nuclear@1 557 x0--;
nuclear@1 558 out[i] =MULT_NORM((w[0]*T[0]+w[1]*T[1])*init->scale);
nuclear@1 559 x0[0] =MULT_NORM((w[0]*T[1]-w[1]*T[0])*init->scale);
nuclear@1 560 w+=2;
nuclear@1 561 T+=2;
nuclear@1 562 }
nuclear@1 563 }