[xiph-commits] r12948 - in trunk/ghost: . monty monty/lift
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Mon May 14 09:11:03 PDT 2007
Author: xiphmont
Date: 2007-05-14 09:11:02 -0700 (Mon, 14 May 2007)
New Revision: 12948
Added:
trunk/ghost/monty/
trunk/ghost/monty/lift/
trunk/ghost/monty/lift/liftsearch.c
trunk/ghost/monty/lift/smallft.c
trunk/ghost/monty/lift/smallft.h
Log:
No reason not to have some of the working lift searcher in SVN
Added: trunk/ghost/monty/lift/liftsearch.c
===================================================================
--- trunk/ghost/monty/lift/liftsearch.c (rev 0)
+++ trunk/ghost/monty/lift/liftsearch.c 2007-05-14 16:11:02 UTC (rev 12948)
@@ -0,0 +1,399 @@
+#include <math.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <stdio.h>
+#include <string.h>
+#include "smallft.h"
+
+int fixed_bits = 12;
+int num_coefficients = 14;
+int p_delay = 2;
+int q_delay = 3;
+double TB = .1;
+
+#define MAX_PZ 32
+#define FS 256
+drft_lookup fft;
+
+typedef struct local_minimum{
+ int *coefficients;
+ double cost;
+} lmin_t;
+
+lmin_t *minimum_list = NULL;
+int minimum_list_size = 0;
+
+static inline double todB(double x){
+ return ((x)==0?-400:log((x)*(x))*4.34294480f);
+}
+
+int log_minimum(int *c, double cost){
+ int i,j;
+
+ for(i=0;i<minimum_list_size;i++){
+ lmin_t *min = minimum_list+i;
+ for(j=0;j<num_coefficients;j++)
+ if(min->coefficients[j] != c[j])break;
+ if(j==num_coefficients) return i;
+ }
+
+ if(!minimum_list){
+ minimum_list = malloc(sizeof(*minimum_list));
+ }else{
+ minimum_list = realloc(minimum_list, (minimum_list_size+1)*sizeof(*minimum_list));
+ }
+
+ int *mc = malloc(num_coefficients * sizeof(*mc));
+ minimum_list[minimum_list_size].coefficients = mc;
+ for(j=0;j<num_coefficients;j++) mc[j] = c[j];
+ minimum_list_size++;
+
+ return minimum_list_size-1;
+}
+
+double fixed_to_double(int f){
+ return f*pow(2.,(double)-fixed_bits);
+}
+
+typedef struct {
+ int n;
+ double a[MAX_PZ]; // normalized filter; a[0] is actually a_1
+ double b[MAX_PZ+1];
+ double pre[MAX_PZ*2];
+} iir_t;
+
+double Az2(iir_t *f, double v){
+ double iir = v;
+ int i;
+ int n = f->n;
+
+ for(i=0;i<n;i++)
+ iir -= f->pre[i*2+1] * f->a[i];
+
+ double fir = f->b[0] * iir;
+ for(i=0;i<n;i++)
+ fir += f->pre[i*2+1] * f->b[i+1];
+
+ for(i = n*2-1; i>0; i--)
+ f->pre[i] = f->pre[i-1];
+ f->pre[0] = iir;
+
+ return fir;
+}
+
+void mag_e(double *d){
+ int i;
+
+ drft_forward(&fft,d);
+
+ d[0] = d[0]*d[0];
+ for(i=2;i<FS;i+=2)
+ d[i>>1] =d[i-1]*d[i-1] + d[i]*d[i];
+ d[FS/2] = d[FS-1]*d[FS-1];
+
+}
+
+void lift_run(int *c, double *lp, double *hp, int dump){
+ iir_t P;
+ iir_t Q;
+ int i,j;
+ int offset = (abs(p_delay) + abs(q_delay))*2+1;
+
+ memset(&P,0,sizeof(P));
+ memset(&Q,0,sizeof(Q));
+ lp[offset] = 1.;
+ hp[offset+1] = 1.;
+
+ P.n = Q.n = ((num_coefficients/2) -1)/2;
+
+ for(i=0,j=0; i<P.n; i++, j++)
+ P.a[i] = fixed_to_double(c[j]);
+ for(i=0; i<P.n+1; i++, j++)
+ P.b[i] = fixed_to_double(c[j]);
+
+ for(i=0; i<Q.n; i++, j++)
+ Q.a[i] = fixed_to_double(c[j]);
+ for(i=0; i<Q.n+1; i++, j++)
+ Q.b[i] = fixed_to_double(c[j]);
+
+ i=0; j=0;
+ if(p_delay>0) i += p_delay*2;
+ if(p_delay<0) j += -p_delay*2;
+
+ for(; i<FS && j<FS; i++, j++)
+ hp[j] -= Az2(&P, lp[i]);
+
+ i=0; j=0;
+ if(q_delay>0) i += q_delay*2;
+ if(q_delay<0) j += -q_delay*2;
+
+ for(; i<FS && j<FS; i++, j++)
+ lp[j] += .5*Az2(&Q, hp[i]);
+
+ // energy
+ mag_e(hp);
+ mag_e(lp);
+
+ if(dump){
+ fprintf(stdout,"Filter result: \n");
+
+ fprintf(stdout," P: (");
+ for(i=0;i<P.n+1;i++)
+ fprintf(stdout,"%s%f",(i==0?"":", "),P.b[i]);
+ fprintf(stdout,") / (1");
+ for(i=0;i<P.n;i++)
+ fprintf(stdout,", %f",P.a[i]);
+ fprintf(stdout,")\n");
+
+ fprintf(stdout," Q: (");
+ for(i=0;i<Q.n+1;i++)
+ fprintf(stdout,"%s%f",(i==0?"":", "),Q.b[i]);
+ fprintf(stdout,") / (1");
+ for(i=0;i<Q.n;i++)
+ fprintf(stdout,", %f",Q.a[i]);
+ fprintf(stdout,")\n");
+
+ fprintf(stdout," response: {\n");
+ for(i=0;i<=FS/2;i++)
+ fprintf(stdout,"%f %f\n",(double)i/FS, todB(lp[i])/2);
+ fprintf(stdout,"\n");
+ for(i=0;i<=FS/2;i++)
+ fprintf(stdout,"%f %f\n",(double)i/FS, todB(hp[i])/2);
+ fprintf(stdout,"}\n\n");
+ }
+
+}
+
+
+double lift_cost(int *c){
+ double *hp = calloc (FS,sizeof(*hp));
+ double *lp = calloc (FS,sizeof(*lp));
+ int i;
+
+ lift_run(c, lp, hp, 0);
+
+ // cost
+ double cost = 0;
+ int TBhi = (int)rint(FS/4*(1.+TB));
+ int TBlo = (int)rint(FS/4*(1.-TB));
+
+ // energy in stopband
+ for(i=TBhi; i<=FS/2; i++)
+ cost += lp[i];
+ for(i=0; i<=TBlo; i++)
+ cost += hp[i];
+
+ // deviation from 0dB in passband
+ //for(i=0; i<=TBlo; i++)
+ // cost += fabs(lp[i]-1.)/4096.;
+ //for(i=TBhi; i<=FS/2; i++)
+ // cost += fabs(hp[i]-1.)/4096.;
+
+ // overshoot from 0dB in TB
+ //for(i=TBlo; i<=TBhi; i++){
+ // double over = lp[i]-1.;
+ // if(over > 0.) cost+=over;
+ // over = hp[i]-1.;
+ // if(over > 0.) cost+=over;
+ //}
+
+ free(hp);
+ free(lp);
+
+ return cost;
+}
+
+void lift_dump(int *c){
+ double *hp = calloc (FS,sizeof(*hp));
+ double *lp = calloc (FS,sizeof(*lp));
+ int i;
+
+ lift_run(c, lp, hp, 1);
+
+ free(hp);
+ free(lp);
+}
+
+int walk_to_minimum(int *c){
+ int last_change = 0;
+ int cur_dim = 0;
+ double cost = lift_cost(c);
+
+ while(1){
+ c[cur_dim]++;
+ double up_cost = lift_cost(c);
+ c[cur_dim]-=2;
+ double down_cost = lift_cost(c);
+ c[cur_dim]++;
+
+ if(up_cost<cost && up_cost<down_cost){
+ c[cur_dim]++;
+ last_change = cur_dim;
+ cost = up_cost;
+
+ // keep walking until we hit the minimum in this dim
+ /*while(1){
+ c[cur_dim]++;
+ up_cost = lift_cost(c);
+ if(up_cost<cost){
+ cost = up_cost;
+ }else{
+ c[cur_dim]--;
+ break;
+ }
+ }*/
+
+ }else if (down_cost<cost){
+ c[cur_dim]--;
+ last_change = cur_dim;
+ cost = down_cost;
+
+ // keep walking until we hit the minimum in this dim
+ /*
+ while(1){
+ c[cur_dim]--;
+ down_cost = lift_cost(c);
+ if(down_cost<cost){
+ cost = down_cost;
+ }else{
+ c[cur_dim]++;
+ break;
+ }
+ }*/
+
+ }else{
+ if(cur_dim == last_change)break;
+ }
+
+ cur_dim++;
+ if(cur_dim >= num_coefficients)
+ cur_dim = 0;
+
+ fprintf(stderr,"\rwalking... current cost: %g ",cost);
+ }
+ fprintf(stderr,"\r ",cost);
+
+ log_minimum(c,cost);
+}
+
+
+int walk_to_minimum2(int *c){
+ double cost = lift_cost(c);
+ int i;
+
+ while(1){
+ double up_cost[num_coefficients];
+ double down_cost[num_coefficients];
+ double min_cost = cost;
+ int min_cost_i = -1;
+ int min_cost_d = -1;
+
+ for(i=0;i<num_coefficients;i++){
+ c[i]++;
+ up_cost[i] = lift_cost(c);
+ c[i]-=2;
+ down_cost[i] = lift_cost(c);
+ c[i]++;
+ }
+
+ for(i=0;i<num_coefficients;i++){
+ if(up_cost[i] < cost){
+ cost = up_cost[i];
+ min_cost_i = i;
+ min_cost_d = 1;
+ }
+ if(down_cost[i] < cost){
+ cost = down_cost[i];
+ min_cost_i = i;
+ min_cost_d = -1;
+ }
+ }
+
+ if(min_cost_i == -1) break;
+ c[min_cost_i] += min_cost_d;
+
+ fprintf(stderr,"\rwalking... current cost: %f ",cost);
+ }
+ fprintf(stderr,"\r ",cost);
+
+ log_minimum(c,cost);
+}
+
+
+void bisect_Ncube_recurse(int *lowrange, int *highrange, int min_num){
+
+
+
+
+
+}
+
+double n_choose_m(int n, int m){
+ double work[n+1];
+ int i,j;
+
+ work[0]=1;
+
+ for(i=0;i<n;i++){
+ work[i+1] = 1;
+ for(j=i;j>0;j--)
+ work[j] += work[j-1];
+ }
+
+ return work[m];
+}
+
+double def_a(int n, int sub, int delay){
+ double ret = n_choose_m(n,sub+1);
+ int i;
+
+ for(i=1;i<=sub+1;i++)
+ ret *= (n - delay - i +.5) / (delay + i + .5);
+
+ return ret;
+}
+
+int main(int argc, char *argv[]){
+ int *c;
+ double *fc;
+ int i,j;
+
+ drft_init(&fft,FS);
+
+ int order = 8;
+ num_coefficients = (order*2+1)*2;
+ TB = .1;
+ p_delay = order-1;
+ q_delay = order;
+
+ c = alloca(num_coefficients * sizeof(*c));
+ fc = alloca(num_coefficients * sizeof(*fc));
+
+ fixed_bits = 8;
+ for(i=0,j=0;i<order;i++,j++)
+ fc[j] = def_a(order, i, p_delay);
+ for(i=0;i<order;i++,j++)
+ fc[j] = def_a(order, order-i-1, p_delay);
+ fc[j++]=1.;
+ for(i=0;i<order;i++,j++)
+ fc[j] = def_a(order, i, p_delay);
+ for(i=0;i<order;i++,j++)
+ fc[j] = def_a(order, order-i-1, p_delay);
+ fc[j++]=1.;
+
+ for(i=0;i<num_coefficients;i++)
+ c[i] = rint(fc[i] * pow(2, fixed_bits));
+
+ walk_to_minimum(c);
+
+ for(j=9;j<=20;j++){
+ fixed_bits = j;
+ for(i=0;i<num_coefficients;i++)
+ c[i] <<=1;
+ walk_to_minimum(c);
+ }
+
+ lift_dump(c);
+
+ return 0;
+}
Added: trunk/ghost/monty/lift/smallft.c
===================================================================
--- trunk/ghost/monty/lift/smallft.c (rev 0)
+++ trunk/ghost/monty/lift/smallft.c 2007-05-14 16:11:02 UTC (rev 12948)
@@ -0,0 +1,335 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggGhost 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 OggGhost SOURCE CODE IS (C) COPYRIGHT 1994-2007 *
+ * the Xiph.Org FOundation http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: unnormalized powers-of-two forward 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, R_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 int drfti1(int n, double *wa, int *ifac){
+ static int ntryh[2] = { 4,2 };
+ static double tpi = 6.28318530717958648f;
+ double 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 < 2)
+ ntry=ntryh[j];
+ else // powers of two only in this one, sorry
+ return -1;
+
+ 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=(double)ld*argh;
+ fi=0.;
+ for (ii=2;ii<ido;ii+=2){
+ fi+=1.;
+ arg=fi*argld;
+ wa[i++]=cos(arg);
+ wa[i++]=sin(arg);
+ }
+ is+=ido;
+ }
+ l1=l2;
+ }
+ return 0;
+}
+
+static int fdrffti(int n, double *wsave, int *ifac){
+ if (n == 1) return 0;
+ return drfti1(n, wsave, ifac);
+}
+
+static void dradf2(int ido,int l1,double *cc,double *ch,double *wa1){
+ int i,k;
+ double 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,double *cc,double *ch,double *wa1,
+ double *wa2,double *wa3){
+ static double hsqt2 = .70710678118654752f;
+ int i,k,t0,t1,t2,t3,t4,t5,t6;
+ double 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 drftf1(int n,double *c,double *ch,double *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(na!=0)goto L103;
+
+ dradf2(ido,l1,c,ch,wa+iw-1);
+ goto L110;
+
+ L103:
+ dradf2(ido,l1,ch,c,wa+iw-1);
+
+ L110:
+ l2=l1;
+ }
+
+ if(na==1)return;
+
+ for(i=0;i<n;i++)c[i]=ch[i];
+}
+
+void drft_forward(drft_lookup *l,double *data){
+ double work[l->n];
+ if(l->n==1)return;
+ drftf1(l->n,data,work,l->trigcache,l->splitcache);
+}
+
+int drft_init(drft_lookup *l,int n){
+ l->n=n;
+ l->trigcache=calloc(2*n,sizeof(*l->trigcache));
+ l->splitcache=calloc(32,sizeof(*l->splitcache));
+ return 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/lift/smallft.h
===================================================================
--- trunk/ghost/monty/lift/smallft.h (rev 0)
+++ trunk/ghost/monty/lift/smallft.h 2007-05-14 16:11:02 UTC (rev 12948)
@@ -0,0 +1,44 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggGhost 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 OggGhost SOURCE CODE IS (C) COPYRIGHT 1994-2007 *
+ * the Xiph.Org FOundation http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: unnormalized powers-of-two forward 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, R_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
+ */
+
+#ifndef _V_SMFT_H_
+#define _V_SMFT_H_
+
+typedef struct {
+ int n;
+ double *trigcache;
+ int *splitcache;
+} drft_lookup;
+
+extern void drft_forward(drft_lookup *l,double *data);
+extern int drft_init(drft_lookup *l,int n);
+extern void drft_clear(drft_lookup *l);
+
+#endif
More information about the commits
mailing list