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 dradfg(int ido,int ip,int l1,int idl1,double *cc,double *c1,
253 double *c2,double *ch,double *ch2,double *wa){
255 static double tpi=6.28318530717958647692528676655900577;
256 int idij,ipph,i,j,k,l,ic,ik,is;
257 int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
258 double dc2,ai1,ai2,ar1,ar2,ds2;
260 double dcp,arg,dsp,ar1h,ar2h;
274 for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
300 ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
301 ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
317 ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
318 ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
341 c1[t5-1]=ch[t5-1]+ch[t6-1];
342 c1[t6-1]=ch[t5]-ch[t6];
343 c1[t5]=ch[t5]+ch[t6];
344 c1[t6]=ch[t6-1]-ch[t5-1];
360 c1[t5-1]=ch[t5-1]+ch[t6-1];
361 c1[t6-1]=ch[t5]-ch[t6];
362 c1[t5]=ch[t5]+ch[t6];
363 c1[t6]=ch[t6-1]-ch[t5-1];
372 for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
384 c1[t3]=ch[t3]+ch[t4];
385 c1[t4]=ch[t4]-ch[t3];
397 ar1h=dcp*ar1-dsp*ai1;
405 for(ik=0;ik<idl1;ik++){
406 ch2[t4++]=c2[ik]+ar1*c2[t7++];
407 ch2[t5++]=ai1*c2[t6++];
421 ar2h=dc2*ar2-ds2*ai2;
429 for(ik=0;ik<idl1;ik++){
430 ch2[t6++]+=ar2*c2[t8++];
431 ch2[t7++]+=ai2*c2[t9++];
440 for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
450 for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
511 cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
512 cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
513 cc[i+t7]=ch[i+t8]+ch[i+t9];
514 cc[ic+t6]=ch[i+t9]-ch[i+t8];
541 cc[t7-1]=ch[t8-1]+ch[t9-1];
542 cc[t6-1]=ch[t8-1]-ch[t9-1];
543 cc[t7]=ch[t8]+ch[t9];
544 cc[t6]=ch[t9]-ch[t8];
554 static void drftf1(int n,double *c,double *ch,double *wa,int *ifac){
557 int ip,iw,ido,idl1,ix2,ix3;
564 for(k1=0;k1<nf;k1++){
578 dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
580 dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
587 dradf2(ido,l1,c,ch,wa+iw-1);
591 dradf2(ido,l1,ch,c,wa+iw-1);
598 dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
603 dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
612 for(i=0;i<n;i++)c[i]=ch[i];
615 static void dradb2(int ido,int l1,double *cc,double *ch,double *wa1){
616 int i,k,t0,t1,t2,t3,t4,t5,t6;
625 ch[t1]=cc[t2]+cc[t3+t2];
626 ch[t1+t0]=cc[t2]-cc[t3+t2];
644 ch[t3-1]=cc[t4-1]+cc[t5-1];
645 tr2=cc[t4-1]-cc[t5-1];
646 ch[t3]=cc[t4]-cc[t5];
648 ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
649 ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
660 ch[t1]=cc[t2]+cc[t2];
661 ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
667 static void dradb3(int ido,int l1,double *cc,double *ch,double *wa1,
669 static double taur = -.5;
670 static double taui = .86602540378443864676372317075293618;
671 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
672 double ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
681 tr2=cc[t3-1]+cc[t3-1];
682 cr2=cc[t5]+(taur*tr2);
684 ci3=taui*(cc[t3]+cc[t3]);
709 tr2=cc[t5-1]+cc[t6-1];
710 cr2=cc[t7-1]+(taur*tr2);
711 ch[t8-1]=cc[t7-1]+tr2;
713 ci2=cc[t7]+(taur*ti2);
715 cr3=taui*(cc[t5-1]-cc[t6-1]);
716 ci3=taui*(cc[t5]+cc[t6]);
721 ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
722 ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
723 ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
724 ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
730 static void dradb4(int ido,int l1,double *cc,double *ch,double *wa1,
731 double *wa2,double *wa3){
732 static double sqrt2=1.4142135623730950488016887242097;
733 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
734 double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
744 tr3=cc[t4-1]+cc[t4-1];
746 tr1=cc[t3]-cc[(t4+=t6)-1];
761 t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
773 tr1=cc[t2-1]-cc[t5-1];
774 tr2=cc[t2-1]+cc[t5-1];
775 ti4=cc[t3-1]-cc[t4-1];
776 tr3=cc[t3-1]+cc[t4-1];
786 ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
787 ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
788 ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
789 ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
790 ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
791 ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
796 if(ido%2 == 1)return;
808 tr1=cc[t1-1]-cc[t4-1];
809 tr2=cc[t1-1]+cc[t4-1];
811 ch[t5+=t0]=sqrt2*(tr1-ti1);
813 ch[t5+=t0]=-sqrt2*(tr1+ti1);
821 static void dradbg(int ido,int ip,int l1,int idl1,double *cc,double *c1,
822 double *c2,double *ch,double *ch2,double *wa){
823 static double tpi=6.28318530717958647692528676655900577;
824 int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
826 double dc2,ai1,ai2,ar1,ar2,ds2;
828 double dcp,arg,dsp,ar1h,ar2h;
880 ch[t3]=cc[t6-1]+cc[t6-1];
881 ch[t4]=cc[t6]+cc[t6];
889 if (ido == 1)goto L116;
913 ch[t5-1]=cc[t9-1]+cc[t11-1];
914 ch[t6-1]=cc[t9-1]-cc[t11-1];
915 ch[t5]=cc[t9]-cc[t11];
916 ch[t6]=cc[t9]+cc[t11];
947 ch[t5-1]=cc[t11-1]+cc[t12-1];
948 ch[t6-1]=cc[t11-1]-cc[t12-1];
949 ch[t5]=cc[t11]-cc[t12];
950 ch[t6]=cc[t11]+cc[t12];
969 ar1h=dcp*ar1-dsp*ai1;
977 for(ik=0;ik<idl1;ik++){
978 c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
979 c2[t5++]=ai1*ch2[t8++];
991 ar2h=dc2*ar2-ds2*ai2;
998 for(ik=0;ik<idl1;ik++){
999 c2[t4++]+=ar2*ch2[t11++];
1000 c2[t5++]+=ai2*ch2[t12++];
1006 for(j=1;j<ipph;j++){
1009 for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
1014 for(j=1;j<ipph;j++){
1020 ch[t3]=c1[t3]-c1[t4];
1021 ch[t4]=c1[t3]+c1[t4];
1027 if(ido==1)goto L132;
1028 if(nbd<l1)goto L128;
1032 for(j=1;j<ipph;j++){
1040 for(i=2;i<ido;i+=2){
1043 ch[t5-1]=c1[t5-1]-c1[t6];
1044 ch[t6-1]=c1[t5-1]+c1[t6];
1045 ch[t5]=c1[t5]+c1[t6-1];
1046 ch[t6]=c1[t5]-c1[t6-1];
1057 for(j=1;j<ipph;j++){
1062 for(i=2;i<ido;i+=2){
1068 ch[t5-1]=c1[t5-1]-c1[t6];
1069 ch[t6-1]=c1[t5-1]+c1[t6];
1070 ch[t5]=c1[t5]+c1[t6-1];
1071 ch[t6]=c1[t5]-c1[t6-1];
1081 for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
1092 if(nbd>l1)goto L139;
1101 for(i=2;i<ido;i+=2){
1106 c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
1107 c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
1124 for(i=2;i<ido;i+=2){
1127 c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
1128 c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
1135 static void drftb1(int n, double *c, double *ch, double *wa, int *ifac){
1138 int nf,ip,iw,ix2,ix3,ido,idl1;
1145 for(k1=0;k1<nf;k1++){
1155 dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
1157 dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
1165 dradb2(ido,l1,ch,c,wa+iw-1);
1167 dradb2(ido,l1,c,ch,wa+iw-1);
1176 dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
1178 dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
1183 /* The radix five case can be translated later..... */
1184 /* if(ip!=5)goto L112;
1190 dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
1192 dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
1198 dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
1200 dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
1210 for(i=0;i<n;i++)c[i]=ch[i];
1213 void drft_forward(drft_lookup *l,double *data){
1215 drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
1218 void drft_backward(drft_lookup *l,double *data){
1222 drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
1223 for(i=0;i<l->n;i++)data[i]*=n1;
1226 void drft_init(drft_lookup *l,int n){
1228 l->trigcache=calloc(3*n,sizeof(double));
1229 l->splitcache=calloc(32,sizeof(int));
1230 fdrffti(n, l->trigcache, l->splitcache);
1233 void drft_clear(drft_lookup *l){
1235 if(l->trigcache)free(l->trigcache);
1236 if(l->splitcache)free(l->splitcache);
1237 memset(l,0,sizeof(drft_lookup));