[xiph-commits] r17897 - in trunk/ghost/monty: . window
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Thu Mar 17 23:42:25 PDT 2011
Author: xiphmont
Date: 2011-03-17 23:42:25 -0700 (Thu, 17 Mar 2011)
New Revision: 17897
Added:
trunk/ghost/monty/window/
trunk/ghost/monty/window/mdct.c
trunk/ghost/monty/window/mdct.h
trunk/ghost/monty/window/noise.c
trunk/ghost/monty/window/smallft.c
trunk/ghost/monty/window/smallft.h
trunk/ghost/monty/window/window.c
Log:
Add some window comparison test code that was sitting in local
Added: trunk/ghost/monty/window/mdct.c
===================================================================
--- trunk/ghost/monty/window/mdct.c (rev 0)
+++ trunk/ghost/monty/window/mdct.c 2011-03-18 06:42:25 UTC (rev 17897)
@@ -0,0 +1,564 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
+ * 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-2002 *
+ * by the XIPHOPHORUS Company http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: normalized modified discrete cosine transform
+ power of two length transform only [64 <= n ]
+ last mod: $Id: mdct.c 7187 2004-07-20 07:24:27Z xiphmont $
+
+ 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
+
+ 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/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 <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "vorbis/codec.h"
+#include "mdct.h"
+#include "os.h"
+#include "misc.h"
+
+/* build lookups for trig functions; also pre-figure scaling and
+ some window function algebra. */
+
+void mdct_init(mdct_lookup *lookup,int n){
+ int *bitrev=_ogg_malloc(sizeof(*bitrev)*(n/4));
+ DATA_TYPE *T=_ogg_malloc(sizeof(*T)*(n+n/4));
+
+ int i;
+ int n2=n>>1;
+ int log2n=lookup->log2n=rint(log((float)n)/log(2.f));
+ lookup->n=n;
+ lookup->trig=T;
+ lookup->bitrev=bitrev;
+
+/* trig lookups... */
+
+ for(i=0;i<n/4;i++){
+ T[i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i)));
+ T[i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i)));
+ T[n2+i*2]=FLOAT_CONV(cos((M_PI/(2*n))*(2*i+1)));
+ T[n2+i*2+1]=FLOAT_CONV(sin((M_PI/(2*n))*(2*i+1)));
+ }
+ for(i=0;i<n/8;i++){
+ T[n+i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i+2))*.5);
+ T[n+i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i+2))*.5);
+ }
+
+ /* bitreverse lookup... */
+
+ {
+ int mask=(1<<(log2n-1))-1,i,j;
+ int msb=1<<(log2n-2);
+ for(i=0;i<n/8;i++){
+ int acc=0;
+ for(j=0;msb>>j;j++)
+ if((msb>>j)&i)acc|=1<<j;
+ bitrev[i*2]=((~acc)&mask)-1;
+ bitrev[i*2+1]=acc;
+
+ }
+ }
+ lookup->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++)
+ mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i);
+ }
+
+ for(j=0;j<points;j+=32)
+ mdct_butterfly_32(x+j);
+
+}
+
+void mdct_clear(mdct_lookup *l){
+ if(l){
+ if(l->trig)_ogg_free(l->trig);
+ if(l->bitrev)_ogg_free(l->bitrev);
+ memset(l,0,sizeof(*l));
+ }
+}
+
+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;
+
+ do{
+ DATA_TYPE *x0 = x+bit[0];
+ DATA_TYPE *x1 = x+bit[1];
+
+ REG_TYPE r0 = x0[1] - x1[1];
+ REG_TYPE r1 = x0[0] + x1[0];
+ REG_TYPE r2 = MULT_NORM(r1 * T[0] + r0 * T[1]);
+ REG_TYPE r3 = MULT_NORM(r1 * T[1] - r0 * T[0]);
+
+ w1 -= 4;
+
+ r0 = HALVE(x0[1] + x1[1]);
+ r1 = HALVE(x0[0] - x1[0]);
+
+ 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(w0<w1);
+}
+
+void mdct_backward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
+ int n=init->n;
+ int n2=n>>1;
+ int n4=n>>2;
+
+ /* 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 */
+
+ {
+ DATA_TYPE *oX1=out+n2+n4;
+ DATA_TYPE *oX2=out+n2+n4;
+ DATA_TYPE *iX =out;
+ T =init->trig+n2;
+
+ do{
+ oX1-=4;
+
+ oX1[3] = MULT_NORM (iX[0] * T[1] - iX[1] * T[0]);
+ oX2[0] = -MULT_NORM (iX[0] * T[0] + iX[1] * T[1]);
+
+ oX1[2] = MULT_NORM (iX[2] * T[3] - iX[3] * T[2]);
+ oX2[1] = -MULT_NORM (iX[2] * T[2] + iX[3] * T[3]);
+
+ oX1[1] = MULT_NORM (iX[4] * T[5] - iX[5] * T[4]);
+ oX2[2] = -MULT_NORM (iX[4] * T[4] + iX[5] * T[5]);
+
+ oX1[0] = MULT_NORM (iX[6] * T[7] - iX[7] * T[6]);
+ oX2[3] = -MULT_NORM (iX[6] * T[6] + iX[7] * T[7]);
+
+ oX2+=4;
+ iX += 8;
+ T += 8;
+ }while(iX<oX1);
+
+ iX=out+n2+n4;
+ oX1=out+n4;
+ oX2=oX1;
+
+ do{
+ oX1-=4;
+ iX-=4;
+
+ oX2[0] = -(oX1[3] = iX[3]);
+ oX2[1] = -(oX1[2] = iX[2]);
+ oX2[2] = -(oX1[1] = iX[1]);
+ oX2[3] = -(oX1[0] = iX[0]);
+
+ oX2+=4;
+ }while(oX2<iX);
+
+ iX=out+n2+n4;
+ oX1=out+n2+n4;
+ oX2=out+n2;
+ do{
+ oX1-=4;
+ oX1[0]= iX[3];
+ oX1[1]= iX[2];
+ oX1[2]= iX[1];
+ oX1[3]= iX[0];
+ iX+=4;
+ }while(oX1>oX2);
+ }
+}
+
+void mdct_forward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
+ int n=init->n;
+ int n2=n>>1;
+ int n4=n>>2;
+ int n8=n>>3;
+ DATA_TYPE *w=alloca(n*sizeof(*w)); /* forward needs working space */
+ DATA_TYPE *w2=w+n2;
+
+ /* rotate */
+
+ /* window + rotate + step 1 */
+
+ REG_TYPE r0;
+ REG_TYPE r1;
+ DATA_TYPE *x0=in+n2+n4;
+ DATA_TYPE *x1=x0+1;
+ DATA_TYPE *T=init->trig+n2;
+
+ int i=0;
+
+ for(i=0;i<n8;i+=2){
+ x0 -=4;
+ T-=2;
+ r0= x0[2] + x1[0];
+ r1= x0[0] + x1[2];
+ w2[i]= MULT_NORM(r1*T[1] + r0*T[0]);
+ w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
+ x1 +=4;
+ }
+
+ x1=in+1;
+
+ for(;i<n2-n8;i+=2){
+ T-=2;
+ x0 -=4;
+ r0= x0[2] - x1[0];
+ r1= x0[0] - x1[2];
+ w2[i]= MULT_NORM(r1*T[1] + r0*T[0]);
+ w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
+ x1 +=4;
+ }
+
+ x0=in+n;
+
+ for(;i<n2;i+=2){
+ T-=2;
+ x0 -=4;
+ r0= -x0[2] - x1[0];
+ r1= -x0[0] - x1[2];
+ w2[i]= MULT_NORM(r1*T[1] + r0*T[0]);
+ w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
+ x1 +=4;
+ }
+
+
+ mdct_butterflies(init,w+n2,n2);
+ mdct_bitreverse(init,w);
+
+ /* roatate + window */
+
+ T=init->trig+n2;
+ x0=out+n2;
+
+ for(i=0;i<n4;i++){
+ x0--;
+ out[i] =MULT_NORM((w[0]*T[0]+w[1]*T[1])*init->scale);
+ x0[0] =MULT_NORM((w[0]*T[1]-w[1]*T[0])*init->scale);
+ w+=2;
+ T+=2;
+ }
+}
+
Added: trunk/ghost/monty/window/mdct.h
===================================================================
--- trunk/ghost/monty/window/mdct.h (rev 0)
+++ trunk/ghost/monty/window/mdct.h 2011-03-18 06:42:25 UTC (rev 17897)
@@ -0,0 +1,83 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
+ * 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-2002 *
+ * by the XIPHOPHORUS Company http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: modified discrete cosine transform prototypes
+ last mod: $Id: mdct.h 7187 2004-07-20 07:24:27Z xiphmont $
+
+ ********************************************************************/
+
+#ifndef _OGG_mdct_H_
+#define _OGG_mdct_H_
+
+#include "vorbis/codec.h"
+
+
+
+
+
+/*#define MDCT_INTEGERIZED <- be warned there could be some hurt left here*/
+#ifdef MDCT_INTEGERIZED
+
+#define DATA_TYPE int
+#define REG_TYPE register int
+#define TRIGBITS 14
+#define cPI3_8 6270
+#define cPI2_8 11585
+#define cPI1_8 15137
+
+#define FLOAT_CONV(x) ((int)((x)*(1<<TRIGBITS)+.5))
+#define MULT_NORM(x) ((x)>>TRIGBITS)
+#define HALVE(x) ((x)>>1)
+
+#else
+
+#define DATA_TYPE float
+#define REG_TYPE float
+#define cPI3_8 .38268343236508977175F
+#define cPI2_8 .70710678118654752441F
+#define cPI1_8 .92387953251128675613F
+
+#define FLOAT_CONV(x) (x)
+#define MULT_NORM(x) (x)
+#define HALVE(x) ((x)*.5f)
+
+#endif
+
+
+typedef struct {
+ int n;
+ int log2n;
+
+ DATA_TYPE *trig;
+ int *bitrev;
+
+ DATA_TYPE scale;
+} mdct_lookup;
+
+extern void mdct_init(mdct_lookup *lookup,int n);
+extern void mdct_clear(mdct_lookup *l);
+extern void mdct_forward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out);
+extern void mdct_backward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out);
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
Added: trunk/ghost/monty/window/noise.c
===================================================================
--- trunk/ghost/monty/window/noise.c (rev 0)
+++ trunk/ghost/monty/window/noise.c 2011-03-18 06:42:25 UTC (rev 17897)
@@ -0,0 +1,13 @@
+#include <stdio.h>
+
+main(){
+ int i,j;
+ for(i=0;i<44100;i++){
+ int rand1 = (rand()%41)-20;
+ int rand2 = rand()%256;
+ putchar(rand2);
+ putchar(rand1);
+ putchar(rand2);
+ putchar(rand1);
+ }
+}
Added: trunk/ghost/monty/window/smallft.c
===================================================================
--- trunk/ghost/monty/window/smallft.c (rev 0)
+++ trunk/ghost/monty/window/smallft.c 2011-03-18 06:42:25 UTC (rev 17897)
@@ -0,0 +1,1253 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
+ * 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-2002 *
+ * by the XIPHOPHORUS Company http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: *unnormalized* fft transform
+ last mod: $Id: smallft.c 7573 2004-08-16 01:26:52Z conrad $
+
+ ********************************************************************/
+
+/* FFT implementation from OggSquish, minus cosine transforms,
+ * minus all but radix 2/4 case. In Vorbis we only need this
+ * cut-down version.
+ *
+ * To do more than just power-of-two sized vectors, see the full
+ * version I wrote for NetLib.
+ *
+ * Note that the packing is a little strange; rather than the FFT r/i
+ * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
+ * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
+ * FORTRAN version
+ */
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "smallft.h"
+
+static void drfti1(int n, float *wa, int *ifac){
+ static int ntryh[4] = { 4,2,3,5 };
+ static float tpi = 6.28318530717958648f;
+ float arg,argh,argld,fi;
+ int ntry=0,i,j=-1;
+ int k1, l1, l2, ib;
+ int ld, ii, ip, is, nq, nr;
+ int ido, ipm, nfm1;
+ int nl=n;
+ int nf=0;
+
+ L101:
+ j++;
+ if (j < 4)
+ ntry=ntryh[j];
+ else
+ ntry+=2;
+
+ L104:
+ nq=nl/ntry;
+ nr=nl-ntry*nq;
+ if (nr!=0) goto L101;
+
+ nf++;
+ ifac[nf+1]=ntry;
+ nl=nq;
+ if(ntry!=2)goto L107;
+ if(nf==1)goto L107;
+
+ for (i=1;i<nf;i++){
+ ib=nf-i+1;
+ ifac[ib+1]=ifac[ib];
+ }
+ ifac[2] = 2;
+
+ L107:
+ if(nl!=1)goto L104;
+ ifac[0]=n;
+ ifac[1]=nf;
+ argh=tpi/n;
+ is=0;
+ nfm1=nf-1;
+ l1=1;
+
+ if(nfm1==0)return;
+
+ for (k1=0;k1<nfm1;k1++){
+ ip=ifac[k1+2];
+ ld=0;
+ l2=l1*ip;
+ ido=n/l2;
+ ipm=ip-1;
+
+ for (j=0;j<ipm;j++){
+ ld+=l1;
+ i=is;
+ argld=(float)ld*argh;
+ fi=0.f;
+ for (ii=2;ii<ido;ii+=2){
+ fi+=1.f;
+ arg=fi*argld;
+ wa[i++]=cos(arg);
+ wa[i++]=sin(arg);
+ }
+ is+=ido;
+ }
+ l1=l2;
+ }
+}
+
+static void fdrffti(int n, float *wsave, int *ifac){
+
+ if (n == 1) return;
+ drfti1(n, wsave+n, ifac);
+}
+
+static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){
+ int i,k;
+ float ti2,tr2;
+ int t0,t1,t2,t3,t4,t5,t6;
+
+ t1=0;
+ t0=(t2=l1*ido);
+ t3=ido<<1;
+ for(k=0;k<l1;k++){
+ ch[t1<<1]=cc[t1]+cc[t2];
+ ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
+ t1+=ido;
+ t2+=ido;
+ }
+
+ if(ido<2)return;
+ if(ido==2)goto L105;
+
+ t1=0;
+ t2=t0;
+ for(k=0;k<l1;k++){
+ t3=t2;
+ t4=(t1<<1)+(ido<<1);
+ t5=t1;
+ t6=t1+t1;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4-=2;
+ t5+=2;
+ t6+=2;
+ tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
+ ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
+ ch[t6]=cc[t5]+ti2;
+ ch[t4]=ti2-cc[t5];
+ ch[t6-1]=cc[t5-1]+tr2;
+ ch[t4-1]=cc[t5-1]-tr2;
+ }
+ t1+=ido;
+ t2+=ido;
+ }
+
+ if(ido%2==1)return;
+
+ L105:
+ t3=(t2=(t1=ido)-1);
+ t2+=t0;
+ for(k=0;k<l1;k++){
+ ch[t1]=-cc[t2];
+ ch[t1-1]=cc[t3];
+ t1+=ido<<1;
+ t2+=ido;
+ t3+=ido;
+ }
+}
+
+static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,
+ float *wa2,float *wa3){
+ static float hsqt2 = .70710678118654752f;
+ int i,k,t0,t1,t2,t3,t4,t5,t6;
+ float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
+ t0=l1*ido;
+
+ t1=t0;
+ t4=t1<<1;
+ t2=t1+(t1<<1);
+ t3=0;
+
+ for(k=0;k<l1;k++){
+ tr1=cc[t1]+cc[t2];
+ tr2=cc[t3]+cc[t4];
+
+ ch[t5=t3<<2]=tr1+tr2;
+ ch[(ido<<2)+t5-1]=tr2-tr1;
+ ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
+ ch[t5]=cc[t2]-cc[t1];
+
+ t1+=ido;
+ t2+=ido;
+ t3+=ido;
+ t4+=ido;
+ }
+
+ if(ido<2)return;
+ if(ido==2)goto L105;
+
+
+ t1=0;
+ for(k=0;k<l1;k++){
+ t2=t1;
+ t4=t1<<2;
+ t5=(t6=ido<<1)+t4;
+ for(i=2;i<ido;i+=2){
+ t3=(t2+=2);
+ t4+=2;
+ t5-=2;
+
+ t3+=t0;
+ cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
+ ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
+ t3+=t0;
+ cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
+ ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
+ t3+=t0;
+ cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
+ ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
+
+ tr1=cr2+cr4;
+ tr4=cr4-cr2;
+ ti1=ci2+ci4;
+ ti4=ci2-ci4;
+
+ ti2=cc[t2]+ci3;
+ ti3=cc[t2]-ci3;
+ tr2=cc[t2-1]+cr3;
+ tr3=cc[t2-1]-cr3;
+
+ ch[t4-1]=tr1+tr2;
+ ch[t4]=ti1+ti2;
+
+ ch[t5-1]=tr3-ti4;
+ ch[t5]=tr4-ti3;
+
+ ch[t4+t6-1]=ti4+tr3;
+ ch[t4+t6]=tr4+ti3;
+
+ ch[t5+t6-1]=tr2-tr1;
+ ch[t5+t6]=ti1-ti2;
+ }
+ t1+=ido;
+ }
+ if(ido&1)return;
+
+ L105:
+
+ t2=(t1=t0+ido-1)+(t0<<1);
+ t3=ido<<2;
+ t4=ido;
+ t5=ido<<1;
+ t6=ido;
+
+ for(k=0;k<l1;k++){
+ ti1=-hsqt2*(cc[t1]+cc[t2]);
+ tr1=hsqt2*(cc[t1]-cc[t2]);
+
+ ch[t4-1]=tr1+cc[t6-1];
+ ch[t4+t5-1]=cc[t6-1]-tr1;
+
+ ch[t4]=ti1-cc[t1+t0];
+ ch[t4+t5]=ti1+cc[t1+t0];
+
+ t1+=ido;
+ t2+=ido;
+ t4+=t3;
+ t6+=ido;
+ }
+}
+
+static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
+ float *c2,float *ch,float *ch2,float *wa){
+
+ static float tpi=6.283185307179586f;
+ int idij,ipph,i,j,k,l,ic,ik,is;
+ int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
+ float dc2,ai1,ai2,ar1,ar2,ds2;
+ int nbd;
+ float dcp,arg,dsp,ar1h,ar2h;
+ int idp2,ipp2;
+
+ arg=tpi/(float)ip;
+ dcp=cos(arg);
+ dsp=sin(arg);
+ ipph=(ip+1)>>1;
+ ipp2=ip;
+ idp2=ido;
+ nbd=(ido-1)>>1;
+ t0=l1*ido;
+ t10=ip*ido;
+
+ if(ido==1)goto L119;
+ for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
+
+ t1=0;
+ for(j=1;j<ip;j++){
+ t1+=t0;
+ t2=t1;
+ for(k=0;k<l1;k++){
+ ch[t2]=c1[t2];
+ t2+=ido;
+ }
+ }
+
+ is=-ido;
+ t1=0;
+ if(nbd>l1){
+ for(j=1;j<ip;j++){
+ t1+=t0;
+ is+=ido;
+ t2= -ido+t1;
+ for(k=0;k<l1;k++){
+ idij=is-1;
+ t2+=ido;
+ t3=t2;
+ for(i=2;i<ido;i+=2){
+ idij+=2;
+ t3+=2;
+ ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
+ ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
+ }
+ }
+ }
+ }else{
+
+ for(j=1;j<ip;j++){
+ is+=ido;
+ idij=is-1;
+ t1+=t0;
+ t2=t1;
+ for(i=2;i<ido;i+=2){
+ idij+=2;
+ t2+=2;
+ t3=t2;
+ for(k=0;k<l1;k++){
+ ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
+ ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
+ t3+=ido;
+ }
+ }
+ }
+ }
+
+ t1=0;
+ t2=ipp2*t0;
+ if(nbd<l1){
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t5=t3-ido;
+ t6=t4-ido;
+ for(k=0;k<l1;k++){
+ t5+=ido;
+ t6+=ido;
+ c1[t5-1]=ch[t5-1]+ch[t6-1];
+ c1[t6-1]=ch[t5]-ch[t6];
+ c1[t5]=ch[t5]+ch[t6];
+ c1[t6]=ch[t6-1]-ch[t5-1];
+ }
+ }
+ }
+ }else{
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(k=0;k<l1;k++){
+ t5=t3;
+ t6=t4;
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6+=2;
+ c1[t5-1]=ch[t5-1]+ch[t6-1];
+ c1[t6-1]=ch[t5]-ch[t6];
+ c1[t5]=ch[t5]+ch[t6];
+ c1[t6]=ch[t6-1]-ch[t5-1];
+ }
+ t3+=ido;
+ t4+=ido;
+ }
+ }
+ }
+
+L119:
+ for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
+
+ t1=0;
+ t2=ipp2*idl1;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1-ido;
+ t4=t2-ido;
+ for(k=0;k<l1;k++){
+ t3+=ido;
+ t4+=ido;
+ c1[t3]=ch[t3]+ch[t4];
+ c1[t4]=ch[t4]-ch[t3];
+ }
+ }
+
+ ar1=1.f;
+ ai1=0.f;
+ t1=0;
+ t2=ipp2*idl1;
+ t3=(ip-1)*idl1;
+ for(l=1;l<ipph;l++){
+ t1+=idl1;
+ t2-=idl1;
+ ar1h=dcp*ar1-dsp*ai1;
+ ai1=dcp*ai1+dsp*ar1;
+ ar1=ar1h;
+ t4=t1;
+ t5=t2;
+ t6=t3;
+ t7=idl1;
+
+ for(ik=0;ik<idl1;ik++){
+ ch2[t4++]=c2[ik]+ar1*c2[t7++];
+ ch2[t5++]=ai1*c2[t6++];
+ }
+
+ dc2=ar1;
+ ds2=ai1;
+ ar2=ar1;
+ ai2=ai1;
+
+ t4=idl1;
+ t5=(ipp2-1)*idl1;
+ for(j=2;j<ipph;j++){
+ t4+=idl1;
+ t5-=idl1;
+
+ ar2h=dc2*ar2-ds2*ai2;
+ ai2=dc2*ai2+ds2*ar2;
+ ar2=ar2h;
+
+ t6=t1;
+ t7=t2;
+ t8=t4;
+ t9=t5;
+ for(ik=0;ik<idl1;ik++){
+ ch2[t6++]+=ar2*c2[t8++];
+ ch2[t7++]+=ai2*c2[t9++];
+ }
+ }
+ }
+
+ t1=0;
+ for(j=1;j<ipph;j++){
+ t1+=idl1;
+ t2=t1;
+ for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
+ }
+
+ if(ido<l1)goto L132;
+
+ t1=0;
+ t2=0;
+ for(k=0;k<l1;k++){
+ t3=t1;
+ t4=t2;
+ for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
+ t1+=ido;
+ t2+=t10;
+ }
+
+ goto L135;
+
+ L132:
+ for(i=0;i<ido;i++){
+ t1=i;
+ t2=i;
+ for(k=0;k<l1;k++){
+ cc[t2]=ch[t1];
+ t1+=ido;
+ t2+=t10;
+ }
+ }
+
+ L135:
+ t1=0;
+ t2=ido<<1;
+ t3=0;
+ t4=ipp2*t0;
+ for(j=1;j<ipph;j++){
+
+ t1+=t2;
+ t3+=t0;
+ t4-=t0;
+
+ t5=t1;
+ t6=t3;
+ t7=t4;
+
+ for(k=0;k<l1;k++){
+ cc[t5-1]=ch[t6];
+ cc[t5]=ch[t7];
+ t5+=t10;
+ t6+=ido;
+ t7+=ido;
+ }
+ }
+
+ if(ido==1)return;
+ if(nbd<l1)goto L141;
+
+ t1=-ido;
+ t3=0;
+ t4=0;
+ t5=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t2;
+ t3+=t2;
+ t4+=t0;
+ t5-=t0;
+ t6=t1;
+ t7=t3;
+ t8=t4;
+ t9=t5;
+ for(k=0;k<l1;k++){
+ for(i=2;i<ido;i+=2){
+ ic=idp2-i;
+ cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
+ cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
+ cc[i+t7]=ch[i+t8]+ch[i+t9];
+ cc[ic+t6]=ch[i+t9]-ch[i+t8];
+ }
+ t6+=t10;
+ t7+=t10;
+ t8+=ido;
+ t9+=ido;
+ }
+ }
+ return;
+
+ L141:
+
+ t1=-ido;
+ t3=0;
+ t4=0;
+ t5=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t2;
+ t3+=t2;
+ t4+=t0;
+ t5-=t0;
+ for(i=2;i<ido;i+=2){
+ t6=idp2+t1-i;
+ t7=i+t3;
+ t8=i+t4;
+ t9=i+t5;
+ for(k=0;k<l1;k++){
+ cc[t7-1]=ch[t8-1]+ch[t9-1];
+ cc[t6-1]=ch[t8-1]-ch[t9-1];
+ cc[t7]=ch[t8]+ch[t9];
+ cc[t6]=ch[t9]-ch[t8];
+ t6+=t10;
+ t7+=t10;
+ t8+=ido;
+ t9+=ido;
+ }
+ }
+ }
+}
+
+static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){
+ int i,k1,l1,l2;
+ int na,kh,nf;
+ int ip,iw,ido,idl1,ix2,ix3;
+
+ nf=ifac[1];
+ na=1;
+ l2=n;
+ iw=n;
+
+ for(k1=0;k1<nf;k1++){
+ kh=nf-k1;
+ ip=ifac[kh+1];
+ l1=l2/ip;
+ ido=n/l2;
+ idl1=ido*l1;
+ iw-=(ip-1)*ido;
+ na=1-na;
+
+ if(ip!=4)goto L102;
+
+ ix2=iw+ido;
+ ix3=ix2+ido;
+ if(na!=0)
+ dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
+ else
+ dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
+ goto L110;
+
+ L102:
+ if(ip!=2)goto L104;
+ if(na!=0)goto L103;
+
+ dradf2(ido,l1,c,ch,wa+iw-1);
+ goto L110;
+
+ L103:
+ dradf2(ido,l1,ch,c,wa+iw-1);
+ goto L110;
+
+ L104:
+ if(ido==1)na=1-na;
+ if(na!=0)goto L109;
+
+ dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
+ na=1;
+ goto L110;
+
+ L109:
+ dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
+ na=0;
+
+ L110:
+ l2=l1;
+ }
+
+ if(na==1)return;
+
+ for(i=0;i<n;i++)c[i]=ch[i];
+}
+
+static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){
+ int i,k,t0,t1,t2,t3,t4,t5,t6;
+ float ti2,tr2;
+
+ t0=l1*ido;
+
+ t1=0;
+ t2=0;
+ t3=(ido<<1)-1;
+ for(k=0;k<l1;k++){
+ ch[t1]=cc[t2]+cc[t3+t2];
+ ch[t1+t0]=cc[t2]-cc[t3+t2];
+ t2=(t1+=ido)<<1;
+ }
+
+ if(ido<2)return;
+ if(ido==2)goto L105;
+
+ t1=0;
+ t2=0;
+ for(k=0;k<l1;k++){
+ t3=t1;
+ t5=(t4=t2)+(ido<<1);
+ t6=t0+t1;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t5-=2;
+ t6+=2;
+ ch[t3-1]=cc[t4-1]+cc[t5-1];
+ tr2=cc[t4-1]-cc[t5-1];
+ ch[t3]=cc[t4]-cc[t5];
+ ti2=cc[t4]+cc[t5];
+ ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
+ ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
+ }
+ t2=(t1+=ido)<<1;
+ }
+
+ if(ido%2==1)return;
+
+L105:
+ t1=ido-1;
+ t2=ido-1;
+ for(k=0;k<l1;k++){
+ ch[t1]=cc[t2]+cc[t2];
+ ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
+ t1+=ido;
+ t2+=ido<<1;
+ }
+}
+
+static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1,
+ float *wa2){
+ static float taur = -.5f;
+ static float taui = .8660254037844386f;
+ int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
+ float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
+ t0=l1*ido;
+
+ t1=0;
+ t2=t0<<1;
+ t3=ido<<1;
+ t4=ido+(ido<<1);
+ t5=0;
+ for(k=0;k<l1;k++){
+ tr2=cc[t3-1]+cc[t3-1];
+ cr2=cc[t5]+(taur*tr2);
+ ch[t1]=cc[t5]+tr2;
+ ci3=taui*(cc[t3]+cc[t3]);
+ ch[t1+t0]=cr2-ci3;
+ ch[t1+t2]=cr2+ci3;
+ t1+=ido;
+ t3+=t4;
+ t5+=t4;
+ }
+
+ if(ido==1)return;
+
+ t1=0;
+ t3=ido<<1;
+ for(k=0;k<l1;k++){
+ t7=t1+(t1<<1);
+ t6=(t5=t7+t3);
+ t8=t1;
+ t10=(t9=t1+t0)+t0;
+
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6-=2;
+ t7+=2;
+ t8+=2;
+ t9+=2;
+ t10+=2;
+ tr2=cc[t5-1]+cc[t6-1];
+ cr2=cc[t7-1]+(taur*tr2);
+ ch[t8-1]=cc[t7-1]+tr2;
+ ti2=cc[t5]-cc[t6];
+ ci2=cc[t7]+(taur*ti2);
+ ch[t8]=cc[t7]+ti2;
+ cr3=taui*(cc[t5-1]-cc[t6-1]);
+ ci3=taui*(cc[t5]+cc[t6]);
+ dr2=cr2-ci3;
+ dr3=cr2+ci3;
+ di2=ci2+cr3;
+ di3=ci2-cr3;
+ ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
+ ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
+ ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
+ ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
+ }
+ t1+=ido;
+ }
+}
+
+static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,
+ float *wa2,float *wa3){
+ static float sqrt2=1.414213562373095f;
+ int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
+ float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
+ t0=l1*ido;
+
+ t1=0;
+ t2=ido<<2;
+ t3=0;
+ t6=ido<<1;
+ for(k=0;k<l1;k++){
+ t4=t3+t6;
+ t5=t1;
+ tr3=cc[t4-1]+cc[t4-1];
+ tr4=cc[t4]+cc[t4];
+ tr1=cc[t3]-cc[(t4+=t6)-1];
+ tr2=cc[t3]+cc[t4-1];
+ ch[t5]=tr2+tr3;
+ ch[t5+=t0]=tr1-tr4;
+ ch[t5+=t0]=tr2-tr3;
+ ch[t5+=t0]=tr1+tr4;
+ t1+=ido;
+ t3+=t2;
+ }
+
+ if(ido<2)return;
+ if(ido==2)goto L105;
+
+ t1=0;
+ for(k=0;k<l1;k++){
+ t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
+ t7=t1;
+ for(i=2;i<ido;i+=2){
+ t2+=2;
+ t3+=2;
+ t4-=2;
+ t5-=2;
+ t7+=2;
+ ti1=cc[t2]+cc[t5];
+ ti2=cc[t2]-cc[t5];
+ ti3=cc[t3]-cc[t4];
+ tr4=cc[t3]+cc[t4];
+ tr1=cc[t2-1]-cc[t5-1];
+ tr2=cc[t2-1]+cc[t5-1];
+ ti4=cc[t3-1]-cc[t4-1];
+ tr3=cc[t3-1]+cc[t4-1];
+ ch[t7-1]=tr2+tr3;
+ cr3=tr2-tr3;
+ ch[t7]=ti2+ti3;
+ ci3=ti2-ti3;
+ cr2=tr1-tr4;
+ cr4=tr1+tr4;
+ ci2=ti1+ti4;
+ ci4=ti1-ti4;
+
+ ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
+ ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
+ ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
+ ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
+ ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
+ ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
+ }
+ t1+=ido;
+ }
+
+ if(ido%2 == 1)return;
+
+ L105:
+
+ t1=ido;
+ t2=ido<<2;
+ t3=ido-1;
+ t4=ido+(ido<<1);
+ for(k=0;k<l1;k++){
+ t5=t3;
+ ti1=cc[t1]+cc[t4];
+ ti2=cc[t4]-cc[t1];
+ tr1=cc[t1-1]-cc[t4-1];
+ tr2=cc[t1-1]+cc[t4-1];
+ ch[t5]=tr2+tr2;
+ ch[t5+=t0]=sqrt2*(tr1-ti1);
+ ch[t5+=t0]=ti2+ti2;
+ ch[t5+=t0]=-sqrt2*(tr1+ti1);
+
+ t3+=ido;
+ t1+=t2;
+ t4+=t2;
+ }
+}
+
+static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
+ float *c2,float *ch,float *ch2,float *wa){
+ static float tpi=6.283185307179586f;
+ int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
+ t11,t12;
+ float dc2,ai1,ai2,ar1,ar2,ds2;
+ int nbd;
+ float dcp,arg,dsp,ar1h,ar2h;
+ int ipp2;
+
+ t10=ip*ido;
+ t0=l1*ido;
+ arg=tpi/(float)ip;
+ dcp=cos(arg);
+ dsp=sin(arg);
+ nbd=(ido-1)>>1;
+ ipp2=ip;
+ ipph=(ip+1)>>1;
+ if(ido<l1)goto L103;
+
+ t1=0;
+ t2=0;
+ for(k=0;k<l1;k++){
+ t3=t1;
+ t4=t2;
+ for(i=0;i<ido;i++){
+ ch[t3]=cc[t4];
+ t3++;
+ t4++;
+ }
+ t1+=ido;
+ t2+=t10;
+ }
+ goto L106;
+
+ L103:
+ t1=0;
+ for(i=0;i<ido;i++){
+ t2=t1;
+ t3=t1;
+ for(k=0;k<l1;k++){
+ ch[t2]=cc[t3];
+ t2+=ido;
+ t3+=t10;
+ }
+ t1++;
+ }
+
+ L106:
+ t1=0;
+ t2=ipp2*t0;
+ t7=(t5=ido<<1);
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ t6=t5;
+ for(k=0;k<l1;k++){
+ ch[t3]=cc[t6-1]+cc[t6-1];
+ ch[t4]=cc[t6]+cc[t6];
+ t3+=ido;
+ t4+=ido;
+ t6+=t10;
+ }
+ t5+=t7;
+ }
+
+ if (ido == 1)goto L116;
+ if(nbd<l1)goto L112;
+
+ t1=0;
+ t2=ipp2*t0;
+ t7=0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+
+ t7+=(ido<<1);
+ t8=t7;
+ for(k=0;k<l1;k++){
+ t5=t3;
+ t6=t4;
+ t9=t8;
+ t11=t8;
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6+=2;
+ t9+=2;
+ t11-=2;
+ ch[t5-1]=cc[t9-1]+cc[t11-1];
+ ch[t6-1]=cc[t9-1]-cc[t11-1];
+ ch[t5]=cc[t9]-cc[t11];
+ ch[t6]=cc[t9]+cc[t11];
+ }
+ t3+=ido;
+ t4+=ido;
+ t8+=t10;
+ }
+ }
+ goto L116;
+
+ L112:
+ t1=0;
+ t2=ipp2*t0;
+ t7=0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ t7+=(ido<<1);
+ t8=t7;
+ t9=t7;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t8+=2;
+ t9-=2;
+ t5=t3;
+ t6=t4;
+ t11=t8;
+ t12=t9;
+ for(k=0;k<l1;k++){
+ ch[t5-1]=cc[t11-1]+cc[t12-1];
+ ch[t6-1]=cc[t11-1]-cc[t12-1];
+ ch[t5]=cc[t11]-cc[t12];
+ ch[t6]=cc[t11]+cc[t12];
+ t5+=ido;
+ t6+=ido;
+ t11+=t10;
+ t12+=t10;
+ }
+ }
+ }
+
+L116:
+ ar1=1.f;
+ ai1=0.f;
+ t1=0;
+ t9=(t2=ipp2*idl1);
+ t3=(ip-1)*idl1;
+ for(l=1;l<ipph;l++){
+ t1+=idl1;
+ t2-=idl1;
+
+ ar1h=dcp*ar1-dsp*ai1;
+ ai1=dcp*ai1+dsp*ar1;
+ ar1=ar1h;
+ t4=t1;
+ t5=t2;
+ t6=0;
+ t7=idl1;
+ t8=t3;
+ for(ik=0;ik<idl1;ik++){
+ c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
+ c2[t5++]=ai1*ch2[t8++];
+ }
+ dc2=ar1;
+ ds2=ai1;
+ ar2=ar1;
+ ai2=ai1;
+
+ t6=idl1;
+ t7=t9-idl1;
+ for(j=2;j<ipph;j++){
+ t6+=idl1;
+ t7-=idl1;
+ ar2h=dc2*ar2-ds2*ai2;
+ ai2=dc2*ai2+ds2*ar2;
+ ar2=ar2h;
+ t4=t1;
+ t5=t2;
+ t11=t6;
+ t12=t7;
+ for(ik=0;ik<idl1;ik++){
+ c2[t4++]+=ar2*ch2[t11++];
+ c2[t5++]+=ai2*ch2[t12++];
+ }
+ }
+ }
+
+ t1=0;
+ for(j=1;j<ipph;j++){
+ t1+=idl1;
+ t2=t1;
+ for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
+ }
+
+ t1=0;
+ t2=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(k=0;k<l1;k++){
+ ch[t3]=c1[t3]-c1[t4];
+ ch[t4]=c1[t3]+c1[t4];
+ t3+=ido;
+ t4+=ido;
+ }
+ }
+
+ if(ido==1)goto L132;
+ if(nbd<l1)goto L128;
+
+ t1=0;
+ t2=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(k=0;k<l1;k++){
+ t5=t3;
+ t6=t4;
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6+=2;
+ ch[t5-1]=c1[t5-1]-c1[t6];
+ ch[t6-1]=c1[t5-1]+c1[t6];
+ ch[t5]=c1[t5]+c1[t6-1];
+ ch[t6]=c1[t5]-c1[t6-1];
+ }
+ t3+=ido;
+ t4+=ido;
+ }
+ }
+ goto L132;
+
+ L128:
+ t1=0;
+ t2=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t5=t3;
+ t6=t4;
+ for(k=0;k<l1;k++){
+ ch[t5-1]=c1[t5-1]-c1[t6];
+ ch[t6-1]=c1[t5-1]+c1[t6];
+ ch[t5]=c1[t5]+c1[t6-1];
+ ch[t6]=c1[t5]-c1[t6-1];
+ t5+=ido;
+ t6+=ido;
+ }
+ }
+ }
+
+L132:
+ if(ido==1)return;
+
+ for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
+
+ t1=0;
+ for(j=1;j<ip;j++){
+ t2=(t1+=t0);
+ for(k=0;k<l1;k++){
+ c1[t2]=ch[t2];
+ t2+=ido;
+ }
+ }
+
+ if(nbd>l1)goto L139;
+
+ is= -ido-1;
+ t1=0;
+ for(j=1;j<ip;j++){
+ is+=ido;
+ t1+=t0;
+ idij=is;
+ t2=t1;
+ for(i=2;i<ido;i+=2){
+ t2+=2;
+ idij+=2;
+ t3=t2;
+ for(k=0;k<l1;k++){
+ c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
+ c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
+ t3+=ido;
+ }
+ }
+ }
+ return;
+
+ L139:
+ is= -ido-1;
+ t1=0;
+ for(j=1;j<ip;j++){
+ is+=ido;
+ t1+=t0;
+ t2=t1;
+ for(k=0;k<l1;k++){
+ idij=is;
+ t3=t2;
+ for(i=2;i<ido;i+=2){
+ idij+=2;
+ t3+=2;
+ c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
+ c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
+ }
+ t2+=ido;
+ }
+ }
+}
+
+static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){
+ int i,k1,l1,l2;
+ int na;
+ int nf,ip,iw,ix2,ix3,ido,idl1;
+
+ nf=ifac[1];
+ na=0;
+ l1=1;
+ iw=1;
+
+ for(k1=0;k1<nf;k1++){
+ ip=ifac[k1 + 2];
+ l2=ip*l1;
+ ido=n/l2;
+ idl1=ido*l1;
+ if(ip!=4)goto L103;
+ ix2=iw+ido;
+ ix3=ix2+ido;
+
+ if(na!=0)
+ dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
+ else
+ dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
+ na=1-na;
+ goto L115;
+
+ L103:
+ if(ip!=2)goto L106;
+
+ if(na!=0)
+ dradb2(ido,l1,ch,c,wa+iw-1);
+ else
+ dradb2(ido,l1,c,ch,wa+iw-1);
+ na=1-na;
+ goto L115;
+
+ L106:
+ if(ip!=3)goto L109;
+
+ ix2=iw+ido;
+ if(na!=0)
+ dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
+ else
+ dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
+ na=1-na;
+ goto L115;
+
+ L109:
+/* The radix five case can be translated later..... */
+/* if(ip!=5)goto L112;
+
+ ix2=iw+ido;
+ ix3=ix2+ido;
+ ix4=ix3+ido;
+ if(na!=0)
+ dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
+ else
+ dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
+ na=1-na;
+ goto L115;
+
+ L112:*/
+ if(na!=0)
+ dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
+ else
+ dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
+ if(ido==1)na=1-na;
+
+ L115:
+ l1=l2;
+ iw+=(ip-1)*ido;
+ }
+
+ if(na==0)return;
+
+ for(i=0;i<n;i++)c[i]=ch[i];
+}
+
+void drft_forward(drft_lookup *l,float *data){
+ if(l->n==1)return;
+ drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
+}
+
+void drft_backward(drft_lookup *l,float *data){
+ if (l->n==1)return;
+ drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
+}
+
+void drft_init(drft_lookup *l,int n){
+ l->n=n;
+ l->trigcache=calloc(3*n,sizeof(*l->trigcache));
+ l->splitcache=calloc(32,sizeof(*l->splitcache));
+ fdrffti(n, l->trigcache, l->splitcache);
+}
+
+void drft_clear(drft_lookup *l){
+ if(l){
+ if(l->trigcache)free(l->trigcache);
+ if(l->splitcache)free(l->splitcache);
+ memset(l,0,sizeof(*l));
+ }
+}
Added: trunk/ghost/monty/window/smallft.h
===================================================================
--- trunk/ghost/monty/window/smallft.h (rev 0)
+++ trunk/ghost/monty/window/smallft.h 2011-03-18 06:42:25 UTC (rev 17897)
@@ -0,0 +1,32 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
+ * 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-2002 *
+ * by the XIPHOPHORUS Company http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: fft transform
+ last mod: $Id: smallft.h 7187 2004-07-20 07:24:27Z xiphmont $
+
+ ********************************************************************/
+
+#ifndef _V_SMFT_H_
+#define _V_SMFT_H_
+
+typedef struct {
+ int n;
+ float *trigcache;
+ int *splitcache;
+} drft_lookup;
+
+extern void drft_forward(drft_lookup *l,float *data);
+extern void drft_backward(drft_lookup *l,float *data);
+extern void drft_init(drft_lookup *l,int n);
+extern void drft_clear(drft_lookup *l);
+
+#endif
Added: trunk/ghost/monty/window/window.c
===================================================================
--- trunk/ghost/monty/window/window.c (rev 0)
+++ trunk/ghost/monty/window/window.c 2011-03-18 06:42:25 UTC (rev 17897)
@@ -0,0 +1,218 @@
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include <sys/types.h>
+#include "smallft.h"
+
+#define todB(x) ((x)==0?-400.f:log((x)*(x))*4.34294480f)
+#define todB_nn(x) ((x)==0.f?-400.f:log((x))*8.6858896f)
+
+static void _analysis(char *base,int seq, float *data, int n,int dB,
+ off_t offset){
+
+ FILE *f;
+ char buf[80];
+
+ sprintf(buf,"%s_%d.m",base,seq);
+
+ f=fopen(buf,"w");
+ if(f){
+ int i;
+ for(i=0;i<n;i++)
+ if(dB)
+ fprintf(f,"%d %f\n",(int)(i+offset),todB(data[i]));
+ else
+ fprintf(f,"%d %f\n",(int)(i+offset),(data[i]));
+
+ }
+ fclose(f);
+}
+
+void compute_bh4(float *d, int n){
+ int i;
+
+ float a0 = .35875f;
+ float a1 = .48829f;
+ float a2 = .14128f;
+ float a3 = .01168f;
+
+ for(i=0;i<n;i++){
+ float phi =2.*M_PI/n*(i+.5);
+ d[i] = sin((i+.5)/n * M_PI)*sin((i+.5)/n * M_PI);
+ }
+}
+
+void compute_vorbis(float *d, int n){
+ int i;
+ for(i=0;i<n;i++)
+ d[i] = sin(0.5 * M_PI * sin((i+.5)/n * M_PI)*sin((i+.5)/n * M_PI));
+}
+
+float beta(int n, float alpha){
+ return cosh (acosh(pow(10,alpha))/(n-1));
+}
+
+double T(double n, double x){
+ if(fabs(x)<=1){
+ return cos(n*acos(x));
+ }else{
+ return cosh(n*acosh(x));
+ }
+}
+
+void compute_dolphcheb(float *d, int n){
+ int i,k;
+ float a = 6.;
+ int M=n/2;
+ int N=M*2;
+ double b = beta(N,a);
+
+ for(i=0;i<n;i++){
+ double sum=0;
+ for(k=0;k<M;k++)
+ sum += (k&1?-1:1)*T(N,b*cos(M_PI*k/N)) * cos (2*i*k*M_PI/N);
+
+ sum /= T(N,b);
+
+ sum-=.5;
+ d[i]=sum;
+ }
+}
+
+
+void compute_hanning(float *x, int n){
+ float scale = 2*M_PI/n;
+ int i;
+ for(i=0;i<n;i++){
+ float i5 = i+.5;
+ x[i] = .5-.5*cos(scale*i5);
+ }
+}
+
+void compute_inverse(float *d, float *inv, int n){
+ int i;
+
+ float temp[n];
+ compute_vorbis(temp,n);
+
+ for(i=0;i<n/2;i++)
+ inv[i] = d[i]*temp[i] + d[n/2-i-1]*temp[n/2-i-1];
+
+ for(;i<n;i++)
+ inv[i] = d[i]*temp[i] + d[n-(i-n/2)-1]*temp[n-(i-n/2)-1];
+
+
+ for(i=0;i<n;i++)
+ inv[i] = temp[i] /inv[i];
+
+
+}
+
+void multiply(float *acc, float *x, int n){
+ int i;
+ for(i=0;i<n;i++)
+ acc[i]*=x[i];
+}
+
+void lap(float *d, float *dest, int n){
+ int i;
+
+ for(i=0;i<n/2;i++)
+ dest[i] = d[i] + d[n/2-i-1];
+
+ for(;i<n;i++)
+ dest[i] = d[i] + d[n-(i-n/2)-1];
+
+
+}
+
+void fftmag(float *in, float *mag, int n){
+ int i;
+ drft_lookup fft;
+ drft_init(&fft,n);
+ memcpy(mag,in,n*sizeof(*mag));
+ drft_forward(&fft,mag);
+ for(i=1;i+1<n;i+=2)
+ mag[(i+1)>>1] = hypot(mag[i],mag[i+1]);
+}
+
+#define width 256
+main(){
+ int i;
+ float win[width];
+ float inverse[width];
+ float measure[width];
+
+ compute_bh4(win,width);
+ compute_inverse(win,inverse,width);
+
+ _analysis("bh4",0,win,width,0,0);
+ fftmag(win,measure,width);
+ _analysis("bh4_H",0,measure,width/2,1,0);
+
+ _analysis("bh4i",0,inverse,width,0,0);
+ fftmag(inverse,measure,width);
+ _analysis("bh4i_H",0,measure,width/2,1,0);
+
+ multiply(win,inverse,width);
+ _analysis("bh4sq",0,win,width,0,0);
+ fftmag(win,measure,width);
+ _analysis("bh4sq_H",0,measure,width/2,1,0);
+
+ lap(win,inverse,width);
+ _analysis("bh4lap",0,inverse,width,0,0);
+ fftmag(inverse,measure,width);
+ _analysis("bh4lap_H",0,measure,width/2,1,0);
+
+
+ compute_vorbis(win,width);
+ compute_inverse(win,inverse,width);
+
+ _analysis("v",0,win,width,0,0);
+ fftmag(win,measure,width);
+ _analysis("v_H",0,measure,width/2,1,0);
+
+ _analysis("vi",0,inverse,width,0,0);
+ fftmag(inverse,measure,width);
+ _analysis("vi_H",0,measure,width/2,1,0);
+
+ multiply(win,inverse,width);
+ _analysis("vsq",0,win,width,0,0);
+ fftmag(win,measure,width);
+ _analysis("vsq_H",0,measure,width/2,1,0);
+
+
+ compute_dolphcheb(win,width);
+ compute_inverse(win,inverse,width);
+
+ _analysis("d",0,win,width,0,0);
+ fftmag(win,measure,width);
+ _analysis("d_H",0,measure,width/2,1,0);
+
+ _analysis("di",0,inverse,width,0,0);
+ fftmag(inverse,measure,width);
+ _analysis("di_H",0,measure,width/2,1,0);
+
+ multiply(win,inverse,width);
+ _analysis("dsq",0,win,width,0,0);
+ fftmag(win,measure,width);
+ _analysis("dsq_H",0,measure,width/2,1,0);
+
+
+ compute_hanning(win,width);
+ compute_inverse(win,inverse,width);
+
+ _analysis("h",0,win,width,0,0);
+ fftmag(win,measure,width);
+ _analysis("h_H",0,measure,width/2,1,0);
+
+ _analysis("hi",0,inverse,width,0,0);
+ fftmag(inverse,measure,width);
+ _analysis("hi_H",0,measure,width/2,1,0);
+
+ multiply(win,inverse,width);
+ _analysis("hsq",0,win,width,0,0);
+ fftmag(win,measure,width);
+ _analysis("hsq_H",0,measure,width/2,1,0);
+
+}
More information about the commits
mailing list