media/libvorbis/lib/vorbis_mdct.c

Tue, 06 Jan 2015 21:39:09 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Tue, 06 Jan 2015 21:39:09 +0100
branch
TOR_BUG_9701
changeset 8
97036ab72558
permissions
-rw-r--r--

Conditionally force memory storage according to privacy.thirdparty.isolate;
This solves Tor bug #9701, complying with disk avoidance documented in
https://www.torproject.org/projects/torbrowser/design/#disk-avoidance.

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

mercurial