X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=lib%2Fmdct.c;h=2a0ff8d01ba232634d259829eb2c0c25a5f2426d;hb=80e7915d0271dce18deddca3a468a7cb61b2e50f;hp=d24cd6165d9daac4f34fe7278ad401254fd51819;hpb=5e1fb9030d5cededd728d1e5798a3114c631b50e;p=platform%2Fupstream%2Flibvorbis.git diff --git a/lib/mdct.c b/lib/mdct.c index d24cd61..2a0ff8d 100644 --- a/lib/mdct.c +++ b/lib/mdct.c @@ -1,42 +1,46 @@ /******************************************************************** * * * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * - * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY * - * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH * - * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * * * - * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2000 * - * by Monty and the XIPHOPHORUS Company * - * http://www.xiph.org/ * + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 * + * by the Xiph.Org Foundation https://xiph.org/ * * * ******************************************************************** function: normalized modified discrete cosine transform - power of two length transform only [16 <= n ] - last mod: $Id: mdct.c,v 1.18 2000/11/06 00:07:01 xiphmont Exp $ + power of two length transform only [64 <= n ] - Algorithm adapted from _The use of multirate filter banks for coding - of high quality digital audio_, by T. Sporer, K. Brandenburg and - B. Edler, collection of the European Signal Processing Conference - (EUSIPCO), Amsterdam, June 1992, Vol.1, pp 211-214 + Original algorithm adapted long ago from _The use of multirate filter + banks for coding of high quality digital audio_, by T. Sporer, + K. Brandenburg and B. Edler, collection of the European Signal + Processing Conference (EUSIPCO), Amsterdam, June 1992, Vol.1, pp + 211-214 - Note that the below code won't make much sense without the paper; - The presented algorithm was already fairly polished, and the code - once followed it closely. The current code both corrects several - typos in the paper and goes beyond the presented optimizations - (steps 4 through 6 are, for example, entirely eliminated). + The below code implements an algorithm that no longer looks much like + that presented in the paper, but the basic structure remains if you + dig deep enough to see it. - This module DOES NOT INCLUDE code to generate the window function. - Everybody has their own weird favorite including me... I happen to - like the properties of y=sin(2PI*sin^2(x)), but others may vehemently - disagree. + This module DOES NOT INCLUDE code to generate/apply the window + function. Everybody has their own weird favorite including me... I + happen to like the properties of y=sin(.5PI*sin^2(x)), but others may + vehemently disagree. ********************************************************************/ +/* this can also be run as an integer transform by uncommenting a + define in mdct.h; the integerization is a first pass and although + it's likely stable for Vorbis, the dynamic range is constrained and + roundoff isn't done (so it's noisy). Consider it functional, but + only a starting point. There's no point on a machine with an FPU */ + #include #include #include #include +#include "vorbis/codec.h" #include "mdct.h" #include "os.h" #include "misc.h" @@ -45,32 +49,27 @@ some window function algebra. */ void mdct_init(mdct_lookup *lookup,int n){ - int *bitrev=_ogg_malloc(sizeof(int)*(n/4)); - float *trig=_ogg_malloc(sizeof(float)*(n+n/4)); - float *AE=trig; - float *AO=trig+1; - float *BE=AE+n/2; - float *BO=BE+1; - float *CE=BE+n/2; - float *CO=CE+1; - + int *bitrev=_ogg_malloc(sizeof(*bitrev)*(n/4)); + DATA_TYPE *T=_ogg_malloc(sizeof(*T)*(n+n/4)); + int i; - int log2n=lookup->log2n=rint(log(n)/log(2)); + int n2=n>>1; + int log2n=lookup->log2n=rint(log((float)n)/log(2.f)); lookup->n=n; - lookup->trig=trig; + lookup->trig=T; lookup->bitrev=bitrev; - /* trig lookups... */ +/* trig lookups... */ for(i=0;i>j;j++) - if((msb>>j)&i)acc|=1<>j)&i)acc|=1<scale=FLOAT_CONV(4.f/n); +} + +/* 8 point butterfly (in place, 4 register) */ +STIN void mdct_butterfly_8(DATA_TYPE *x){ + REG_TYPE r0 = x[6] + x[2]; + REG_TYPE r1 = x[6] - x[2]; + REG_TYPE r2 = x[4] + x[0]; + REG_TYPE r3 = x[4] - x[0]; + + x[6] = r0 + r2; + x[4] = r0 - r2; + + r0 = x[5] - x[1]; + r2 = x[7] - x[3]; + x[0] = r1 + r0; + x[2] = r1 - r0; + + r0 = x[5] + x[1]; + r1 = x[7] + x[3]; + x[3] = r2 + r3; + x[1] = r2 - r3; + x[7] = r1 + r0; + x[5] = r1 - r0; + +} + +/* 16 point butterfly (in place, 4 register) */ +STIN void mdct_butterfly_16(DATA_TYPE *x){ + REG_TYPE r0 = x[1] - x[9]; + REG_TYPE r1 = x[0] - x[8]; + + x[8] += x[0]; + x[9] += x[1]; + x[0] = MULT_NORM((r0 + r1) * cPI2_8); + x[1] = MULT_NORM((r0 - r1) * cPI2_8); + + r0 = x[3] - x[11]; + r1 = x[10] - x[2]; + x[10] += x[2]; + x[11] += x[3]; + x[2] = r0; + x[3] = r1; + + r0 = x[12] - x[4]; + r1 = x[13] - x[5]; + x[12] += x[4]; + x[13] += x[5]; + x[4] = MULT_NORM((r0 - r1) * cPI2_8); + x[5] = MULT_NORM((r0 + r1) * cPI2_8); + + r0 = x[14] - x[6]; + r1 = x[15] - x[7]; + x[14] += x[6]; + x[15] += x[7]; + x[6] = r0; + x[7] = r1; + + mdct_butterfly_8(x); + mdct_butterfly_8(x+8); +} + +/* 32 point butterfly (in place, 4 register) */ +STIN void mdct_butterfly_32(DATA_TYPE *x){ + REG_TYPE r0 = x[30] - x[14]; + REG_TYPE r1 = x[31] - x[15]; + + x[30] += x[14]; + x[31] += x[15]; + x[14] = r0; + x[15] = r1; + + r0 = x[28] - x[12]; + r1 = x[29] - x[13]; + x[28] += x[12]; + x[29] += x[13]; + x[12] = MULT_NORM( r0 * cPI1_8 - r1 * cPI3_8 ); + x[13] = MULT_NORM( r0 * cPI3_8 + r1 * cPI1_8 ); + + r0 = x[26] - x[10]; + r1 = x[27] - x[11]; + x[26] += x[10]; + x[27] += x[11]; + x[10] = MULT_NORM(( r0 - r1 ) * cPI2_8); + x[11] = MULT_NORM(( r0 + r1 ) * cPI2_8); + + r0 = x[24] - x[8]; + r1 = x[25] - x[9]; + x[24] += x[8]; + x[25] += x[9]; + x[8] = MULT_NORM( r0 * cPI3_8 - r1 * cPI1_8 ); + x[9] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 ); + + r0 = x[22] - x[6]; + r1 = x[7] - x[23]; + x[22] += x[6]; + x[23] += x[7]; + x[6] = r1; + x[7] = r0; + + r0 = x[4] - x[20]; + r1 = x[5] - x[21]; + x[20] += x[4]; + x[21] += x[5]; + x[4] = MULT_NORM( r1 * cPI1_8 + r0 * cPI3_8 ); + x[5] = MULT_NORM( r1 * cPI3_8 - r0 * cPI1_8 ); + + r0 = x[2] - x[18]; + r1 = x[3] - x[19]; + x[18] += x[2]; + x[19] += x[3]; + x[2] = MULT_NORM(( r1 + r0 ) * cPI2_8); + x[3] = MULT_NORM(( r1 - r0 ) * cPI2_8); + + r0 = x[0] - x[16]; + r1 = x[1] - x[17]; + x[16] += x[0]; + x[17] += x[1]; + x[0] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 ); + x[1] = MULT_NORM( r1 * cPI1_8 - r0 * cPI3_8 ); + + mdct_butterfly_16(x); + mdct_butterfly_16(x+16); + +} + +/* N point first stage butterfly (in place, 2 register) */ +STIN void mdct_butterfly_first(DATA_TYPE *T, + DATA_TYPE *x, + int points){ + + DATA_TYPE *x1 = x + points - 8; + DATA_TYPE *x2 = x + (points>>1) - 8; + REG_TYPE r0; + REG_TYPE r1; + + do{ + + r0 = x1[6] - x2[6]; + r1 = x1[7] - x2[7]; + x1[6] += x2[6]; + x1[7] += x2[7]; + x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]); + x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]); + + r0 = x1[4] - x2[4]; + r1 = x1[5] - x2[5]; + x1[4] += x2[4]; + x1[5] += x2[5]; + x2[4] = MULT_NORM(r1 * T[5] + r0 * T[4]); + x2[5] = MULT_NORM(r1 * T[4] - r0 * T[5]); + + r0 = x1[2] - x2[2]; + r1 = x1[3] - x2[3]; + x1[2] += x2[2]; + x1[3] += x2[3]; + x2[2] = MULT_NORM(r1 * T[9] + r0 * T[8]); + x2[3] = MULT_NORM(r1 * T[8] - r0 * T[9]); + + r0 = x1[0] - x2[0]; + r1 = x1[1] - x2[1]; + x1[0] += x2[0]; + x1[1] += x2[1]; + x2[0] = MULT_NORM(r1 * T[13] + r0 * T[12]); + x2[1] = MULT_NORM(r1 * T[12] - r0 * T[13]); + + x1-=8; + x2-=8; + T+=16; + + }while(x2>=x); +} + +/* N/stage point generic N stage butterfly (in place, 2 register) */ +STIN void mdct_butterfly_generic(DATA_TYPE *T, + DATA_TYPE *x, + int points, + int trigint){ + + DATA_TYPE *x1 = x + points - 8; + DATA_TYPE *x2 = x + (points>>1) - 8; + REG_TYPE r0; + REG_TYPE r1; + + do{ + + r0 = x1[6] - x2[6]; + r1 = x1[7] - x2[7]; + x1[6] += x2[6]; + x1[7] += x2[7]; + x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]); + x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]); + + T+=trigint; + + r0 = x1[4] - x2[4]; + r1 = x1[5] - x2[5]; + x1[4] += x2[4]; + x1[5] += x2[5]; + x2[4] = MULT_NORM(r1 * T[1] + r0 * T[0]); + x2[5] = MULT_NORM(r1 * T[0] - r0 * T[1]); + + T+=trigint; + + r0 = x1[2] - x2[2]; + r1 = x1[3] - x2[3]; + x1[2] += x2[2]; + x1[3] += x2[3]; + x2[2] = MULT_NORM(r1 * T[1] + r0 * T[0]); + x2[3] = MULT_NORM(r1 * T[0] - r0 * T[1]); + + T+=trigint; + + r0 = x1[0] - x2[0]; + r1 = x1[1] - x2[1]; + x1[0] += x2[0]; + x1[1] += x2[1]; + x2[0] = MULT_NORM(r1 * T[1] + r0 * T[0]); + x2[1] = MULT_NORM(r1 * T[0] - r0 * T[1]); + + T+=trigint; + x1-=8; + x2-=8; + + }while(x2>=x); +} + +STIN void mdct_butterflies(mdct_lookup *init, + DATA_TYPE *x, + int points){ + + DATA_TYPE *T=init->trig; + int stages=init->log2n-5; + int i,j; + + if(--stages>0){ + mdct_butterfly_first(T,x,points); + } + + for(i=1;--stages>0;i++){ + for(j=0;j<(1<>i)*j,points>>i,4<trig)free(l->trig); - if(l->bitrev)free(l->bitrev); - memset(l,0,sizeof(mdct_lookup)); + if(l->trig)_ogg_free(l->trig); + if(l->bitrev)_ogg_free(l->bitrev); + memset(l,0,sizeof(*l)); } } -static float *_mdct_kernel(float *x, float *w, - int n, int n2, int n4, int n8, - mdct_lookup *init){ - int i; - /* step 2 */ +STIN void mdct_bitreverse(mdct_lookup *init, + DATA_TYPE *x){ + int n = init->n; + int *bit = init->bitrev; + DATA_TYPE *w0 = x; + DATA_TYPE *w1 = x = w0+(n>>1); + DATA_TYPE *T = init->trig+n; - { - float *xA=x+n4; - float *xB=x; - float *w2=w+n4; - float *A=init->trig+n2; + do{ + DATA_TYPE *x0 = x+bit[0]; + DATA_TYPE *x1 = x+bit[1]; - float x0,x1; - i=0; - do{ - x0=*xA - *xB; - w2[i]= *xA++ + *xB++; - x1= *xA - *xB; - A-=4; - w[i++]= x0 * A[0] + x1 * A[1]; - w[i]= x1 * A[0] - x0 * A[1]; - w2[i++]= *xA++ + *xB++; - }while(ilog2n-3;i++){ - int k0=n>>(i+2); - int k1=1<<(i+3); - int wbase=n2-2; - float *A=init->trig; - float *temp; - - for(r=0;r<(k0>>2);r++){ - int w1=wbase; - int w2=w1-(k0>>1); - float AEv= A[0],wA; - float AOv= A[1],wB; - int unroll=i; - wbase-=2; - - k0++; - unroll--; - if(unroll>0){ - s=2<>=1; - do{ - wB =w[w1] -w[w2]; - x[w1] =w[w1] +w[w2]; - wA =w[++w1] -w[++w2]; - x[w1] =w[w1] +w[w2]; - x[w2] =wA*AEv - wB*AOv; - x[w2-1]=wB*AEv + wA*AOv; - w1-=k0; - w2-=k0; - wB =w[w1] -w[w2]; - x[w1] =w[w1] +w[w2]; - wA =w[++w1] -w[++w2]; - x[w1] =w[w1] +w[w2]; - x[w2] =wA*AEv - wB*AOv; - x[w2-1]=wB*AEv + wA*AOv; - w1-=k0; - w2-=k0; - wB =w[w1] -w[w2]; - x[w1] =w[w1] +w[w2]; - wA =w[++w1] -w[++w2]; - x[w1] =w[w1] +w[w2]; - x[w2] =wA*AEv - wB*AOv; - x[w2-1]=wB*AEv + wA*AOv; - w1-=k0; - w2-=k0; - wB =w[w1] -w[w2]; - x[w1] =w[w1] +w[w2]; - wA =w[++w1] -w[++w2]; - x[w1] =w[w1] +w[w2]; - x[w2] =wA*AEv - wB*AOv; - x[w2-1]=wB*AEv + wA*AOv; - w1-=k0; - w2-=k0; - }while(--s); - }else{ - s=2<trig+n; - int *bit=init->bitrev; - float *x1=x; - float *x2=x+n2-1; - i=n8-1; - do{ - int t1=*bit++; - int t2=*bit++; - - float wA=w[t1]-w[t2+1]; - float wB=w[t1-1]+w[t2]; - float wC=w[t1]+w[t2+1]; - float wD=w[t1-1]-w[t2]; - - float wACE=wA* *C; - float wBCE=wB* *C++; - float wACO=wA* *C; - float wBCO=wB* *C++; - - *x1++=( wC+wACO+wBCE)*.5; - *x2--=(-wD+wBCO-wACE)*.5; - *x1++=( wD+wBCO-wACE)*.5; - *x2--=( wC-wACO-wBCE)*.5; - }while(i--); - } - return(x); + w0[0] = r0 + r2; + w1[2] = r0 - r2; + w0[1] = r1 + r3; + w1[3] = r3 - r1; + + x0 = x+bit[2]; + x1 = x+bit[3]; + + r0 = x0[1] - x1[1]; + r1 = x0[0] + x1[0]; + r2 = MULT_NORM(r1 * T[2] + r0 * T[3]); + r3 = MULT_NORM(r1 * T[3] - r0 * T[2]); + + r0 = HALVE(x0[1] + x1[1]); + r1 = HALVE(x0[0] - x1[0]); + + w0[2] = r0 + r2; + w1[0] = r0 - r2; + w0[3] = r1 + r3; + w1[1] = r3 - r1; + + T += 4; + bit += 4; + w0 += 4; + + }while(w0n; - float *x=alloca(sizeof(float)*(n/2)); - float *w=alloca(sizeof(float)*(n/2)); - float *xx; int n2=n>>1; int n4=n>>2; - int n8=n>>3; - int i; - /* window + rotate + step 1 */ + /* rotate */ + + DATA_TYPE *iX = in+n2-7; + DATA_TYPE *oX = out+n2+n4; + DATA_TYPE *T = init->trig+n4; + + do{ + oX -= 4; + oX[0] = MULT_NORM(-iX[2] * T[3] - iX[0] * T[2]); + oX[1] = MULT_NORM (iX[0] * T[3] - iX[2] * T[2]); + oX[2] = MULT_NORM(-iX[6] * T[1] - iX[4] * T[0]); + oX[3] = MULT_NORM (iX[4] * T[1] - iX[6] * T[0]); + iX -= 8; + T += 4; + }while(iX>=in); + + iX = in+n2-8; + oX = out+n2+n4; + T = init->trig+n4; + + do{ + T -= 4; + oX[0] = MULT_NORM (iX[4] * T[3] + iX[6] * T[2]); + oX[1] = MULT_NORM (iX[4] * T[2] - iX[6] * T[3]); + oX[2] = MULT_NORM (iX[0] * T[1] + iX[2] * T[0]); + oX[3] = MULT_NORM (iX[0] * T[0] - iX[2] * T[1]); + iX -= 8; + oX += 4; + }while(iX>=in); + + mdct_butterflies(init,out+n2,n2); + mdct_bitreverse(init,out); + + /* roatate + window */ + { - float tempA,tempB; - int in1=n2+n4-4; - int in2=in1+5; - float *A=init->trig+n2; - - i=0; - - for(i=0;itrig+n2; - in2=1; + do{ + oX1-=4; - for(;itrig+n2; - float *out2=out+n2; - float scale=4./n; - for(i=0;ioX2); } } -void mdct_backward(mdct_lookup *init, float *in, float *out){ +void mdct_forward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){ int n=init->n; - float *x=alloca(sizeof(float)*(n/2)); - float *w=alloca(sizeof(float)*(n/2)); - float *xx; int n2=n>>1; int n4=n>>2; int n8=n>>3; - int i; + DATA_TYPE *w=alloca(n*sizeof(*w)); /* forward needs working space */ + DATA_TYPE *w2=w+n2; - /* rotate + step 1 */ - { - float *inO=in+1; - float *xO= x; - float *A=init->trig+n2; - - for(i=0;itrig+n2; + + int i=0; + + for(i=0;itrig+n2; - int o1=n4,o2=o1-1; - int o3=n4+n2,o4=o3-1; - - for(i=0;itrig+n2; + x0=out+n2; + + for(i=0;iscale); + x0[0] =MULT_NORM((w[0]*T[1]-w[1]*T[0])*init->scale); + w+=2; + T+=2; } }