1 /******************************************************************
2 * CopyPolicy: GNU Public License 2 applies
3 * Copyright (C) 1998 Monty xiphmont@mit.edu, monty@xiph.org
5 * FFT implementation from OggSquish, minus cosine transforms,
6 * minus all but radix 2/4 case. In Vorbis we only need this
7 * cut-down version (used by the encoder, not decoder).
9 * To do more than just power-of-two sized vectors, see the full
10 * version I wrote for NetLib.
12 ******************************************************************/
19 static void drfti1(int n, double *wa, int *ifac){
20 static int ntryh[4] = { 4,2,3,5 };
21 static double tpi = 6.28318530717958647692528676655900577;
22 double arg,argh,argld,fi;
25 int ld, ii, ip, is, nq, nr;
65 for (k1=0;k1<nfm1;k1++){
75 argld=(double)ld*argh;
77 for (ii=2;ii<ido;ii+=2){
89 static void fdrffti(int n, double *wsave, int *ifac){
92 drfti1(n, wsave+n, ifac);
95 static void dradf2(int ido,int l1,double *cc,double *ch,double *wa1){
98 int t0,t1,t2,t3,t4,t5,t6;
104 ch[t1<<1]=cc[t1]+cc[t2];
105 ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
125 tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
126 ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
129 ch[t6-1]=cc[t5-1]+tr2;
130 ch[t4-1]=cc[t5-1]-tr2;
150 static void dradf4(int ido,int l1,double *cc,double *ch,double *wa1,
151 double *wa2,double *wa3){
152 static double hsqt2 = .70710678118654752440084436210485;
153 int i,k,t0,t1,t2,t3,t4,t5,t6;
154 double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
166 ch[t5=t3<<2]=tr1+tr2;
167 ch[(ido<<2)+t5-1]=tr2-tr1;
168 ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
169 ch[t5]=cc[t2]-cc[t1];
192 cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
193 ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
195 cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
196 ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
198 cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
199 ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
229 t2=(t1=t0+ido-1)+(t0<<1);
236 ti1=-hsqt2*(cc[t1]+cc[t2]);
237 tr1=hsqt2*(cc[t1]-cc[t2]);
239 ch[t4-1]=tr1+cc[t6-1];
240 ch[t4+t5-1]=cc[t6-1]-tr1;
242 ch[t4]=ti1-cc[t1+t0];
243 ch[t4+t5]=ti1+cc[t1+t0];
252 static void drftf1(int n,double *c,double *ch,double *wa,int *ifac){
255 int ip,iw,ido,idl1,ix2,ix3;
262 for(k1=0;k1<nf;k1++){
276 dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
278 dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
285 dradf2(ido,l1,c,ch,wa+iw-1);
289 dradf2(ido,l1,ch,c,wa+iw-1);
293 return; /* We're restricted to powers of two. just fail */
301 for(i=0;i<n;i++)c[i]=ch[i];
304 static void dradb2(int ido,int l1,double *cc,double *ch,double *wa1){
305 int i,k,t0,t1,t2,t3,t4,t5,t6;
314 ch[t1]=cc[t2]+cc[t3+t2];
315 ch[t1+t0]=cc[t2]-cc[t3+t2];
333 ch[t3-1]=cc[t4-1]+cc[t5-1];
334 tr2=cc[t4-1]-cc[t5-1];
335 ch[t3]=cc[t4]-cc[t5];
337 ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
338 ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
349 ch[t1]=cc[t2]+cc[t2];
350 ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
356 static void dradb4(int ido,int l1,double *cc,double *ch,double *wa1,
357 double *wa2,double *wa3){
358 static double sqrt2=1.4142135623730950488016887242097;
359 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
360 double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
370 tr3=cc[t4-1]+cc[t4-1];
372 tr1=cc[t3]-cc[(t4+=t6)-1];
387 t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
399 tr1=cc[t2-1]-cc[t5-1];
400 tr2=cc[t2-1]+cc[t5-1];
401 ti4=cc[t3-1]-cc[t4-1];
402 tr3=cc[t3-1]+cc[t4-1];
412 ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
413 ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
414 ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
415 ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
416 ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
417 ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
422 if(ido%2 == 1)return;
434 tr1=cc[t1-1]-cc[t4-1];
435 tr2=cc[t1-1]+cc[t4-1];
437 ch[t5+=t0]=sqrt2*(tr1-ti1);
439 ch[t5+=t0]=-sqrt2*(tr1+ti1);
447 static void drftb1(int n, double *c, double *ch, double *wa, int *ifac){
450 int nf,ip,iw,ix2,ix3,ido,idl1;
457 for(k1=0;k1<nf;k1++){
467 dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
469 dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
477 dradb2(ido,l1,ch,c,wa+iw-1);
479 dradb2(ido,l1,c,ch,wa+iw-1);
484 return; /* silently fail. we only do powers of two in this version */
493 for(i=0;i<n;i++)c[i]=ch[i];
496 void drft_forward(drft_lookup *l,double *data){
498 drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
501 void drft_backward(drft_lookup *l,double *data){
504 drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
505 for(i=0;i<l->n;i++)data[i]/=l->n;
508 void drft_init(drft_lookup *l,int n){
510 l->trigcache=calloc(3*n,sizeof(double));
511 l->splitcache=calloc(32,sizeof(int));
512 fdrffti(n, l->trigcache, l->splitcache);
515 void drft_clear(drft_lookup *l){
517 if(l->trigcache)free(l->trigcache);
518 if(l->splitcache)free(l->splitcache);
519 memset(l,0,sizeof(drft_lookup));