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