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 ******************************************************************/
18 static void drfti1(int n, double *wa, int *ifac){
19 static int ntryh[4] = { 4,2,3,5 };
20 static double tpi = 6.28318530717958647692528676655900577;
21 double arg,argh,argld,fi;
24 int ld, ii, ip, is, nq, nr;
64 for (k1=0;k1<nfm1;k1++){
74 argld=(double)ld*argh;
76 for (ii=2;ii<ido;ii+=2){
88 static void fdrffti(int n, double *wsave, int *ifac){
91 drfti1(n, wsave+n, ifac);
94 static void dradf2(int ido,int l1,double *cc,double *ch,double *wa1){
97 int t0,t1,t2,t3,t4,t5,t6;
103 ch[t1<<1]=cc[t1]+cc[t2];
104 ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
124 tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
125 ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
128 ch[t6-1]=cc[t5-1]+tr2;
129 ch[t4-1]=cc[t5-1]-tr2;
149 static void dradf4(int ido,int l1,double *cc,double *ch,double *wa1,
150 double *wa2,double *wa3){
151 static double hsqt2 = .70710678118654752440084436210485;
152 int i,k,t0,t1,t2,t3,t4,t5,t6;
153 double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
165 ch[t5=t3<<2]=tr1+tr2;
166 ch[(ido<<2)+t5-1]=tr2-tr1;
167 ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
168 ch[t5]=cc[t2]-cc[t1];
191 cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
192 ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
194 cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
195 ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
197 cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
198 ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
228 t2=(t1=t0+ido-1)+(t0<<1);
235 ti1=-hsqt2*(cc[t1]+cc[t2]);
236 tr1=hsqt2*(cc[t1]-cc[t2]);
238 ch[t4-1]=tr1+cc[t6-1];
239 ch[t4+t5-1]=cc[t6-1]-tr1;
241 ch[t4]=ti1-cc[t1+t0];
242 ch[t4+t5]=ti1+cc[t1+t0];
251 static void drftf1(int n,double *c,double *ch,double *wa,int *ifac){
254 int ip,iw,ido,idl1,ix2,ix3;
261 for(k1=0;k1<nf;k1++){
275 dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
277 dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
284 dradf2(ido,l1,c,ch,wa+iw-1);
288 dradf2(ido,l1,ch,c,wa+iw-1);
292 return; /* We're restricted to powers of two. just fail */
300 for(i=0;i<n;i++)c[i]=ch[i];
303 static void fdrfftf(int n,double *r,double *wsave,int *ifac){
305 drftf1(n,r,wsave,wsave+n,ifac);
308 static void dradb2(int ido,int l1,double *cc,double *ch,double *wa1){
309 int i,k,t0,t1,t2,t3,t4,t5,t6;
318 ch[t1]=cc[t2]+cc[t3+t2];
319 ch[t1+t0]=cc[t2]-cc[t3+t2];
337 ch[t3-1]=cc[t4-1]+cc[t5-1];
338 tr2=cc[t4-1]-cc[t5-1];
339 ch[t3]=cc[t4]-cc[t5];
341 ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
342 ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
353 ch[t1]=cc[t2]+cc[t2];
354 ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
360 static void dradb4(int ido,int l1,double *cc,double *ch,double *wa1,
361 double *wa2,double *wa3){
362 static double sqrt2=1.4142135623730950488016887242097;
363 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
364 double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
374 tr3=cc[t4-1]+cc[t4-1];
376 tr1=cc[t3]-cc[(t4+=t6)-1];
391 t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
403 tr1=cc[t2-1]-cc[t5-1];
404 tr2=cc[t2-1]+cc[t5-1];
405 ti4=cc[t3-1]-cc[t4-1];
406 tr3=cc[t3-1]+cc[t4-1];
416 ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
417 ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
418 ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
419 ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
420 ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
421 ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
426 if(ido%2 == 1)return;
438 tr1=cc[t1-1]-cc[t4-1];
439 tr2=cc[t1-1]+cc[t4-1];
441 ch[t5+=t0]=sqrt2*(tr1-ti1);
443 ch[t5+=t0]=-sqrt2*(tr1+ti1);
451 static void drftb1(int n, double *c, double *ch, double *wa, int *ifac){
454 int nf,ip,iw,ix2,ix3,ido,idl1;
461 for(k1=0;k1<nf;k1++){
471 dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
473 dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
481 dradb2(ido,l1,ch,c,wa+iw-1);
483 dradb2(ido,l1,c,ch,wa+iw-1);
488 return; /* silently fail. we only do powers of two in this version */
497 for(i=0;i<n;i++)c[i]=ch[i];
500 static void fdrfftb(int n, double *r, double *wsave, int *ifac){
502 drftb1(n, r, wsave, wsave+n, ifac);
505 void fft_forward(int n, double *buf,double *trigcache,int *splitcache){
508 if(!trigcache || !splitcache){
509 trigcache=calloc(3*n,sizeof(double));
510 splitcache=calloc(32,sizeof(int));
511 fdrffti(n, trigcache, splitcache);
515 fdrfftf(n, buf, trigcache, splitcache);
523 void fft_backward(int n, double *buf, double *trigcache,int *splitcache){
527 if(!trigcache || !splitcache){
528 trigcache=calloc(3*n,sizeof(double));
529 splitcache=calloc(32,sizeof(int));
530 fdrffti(n, trigcache, splitcache);
534 fdrfftb(n, buf, trigcache, splitcache);
536 for(i=0;i<n;i++)buf[i]/=n;
544 void fft_i(int n, double **trigcache, int **splitcache){
545 *trigcache=calloc(3*n,sizeof(double));
546 *splitcache=calloc(32,sizeof(int));
547 fdrffti(n, *trigcache, *splitcache);