vrshoot

annotate libs/vorbis/mdct.c @ 0:b2f14e535253

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