Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | RSS feed

  1. /********************************************************************
  2.  *                                                                  *
  3.  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
  4.  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
  5.  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
  6.  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
  7.  *                                                                  *
  8.  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
  9.  * by the Xiph.Org Foundation https://xiph.org/                     *
  10.  *                                                                  *
  11.  ********************************************************************
  12.  
  13.  function: *unnormalized* fft transform
  14.  
  15.  ********************************************************************/
  16.  
  17. /* FFT implementation from OggSquish, minus cosine transforms,
  18.  * minus all but radix 2/4 case.  In Vorbis we only need this
  19.  * cut-down version.
  20.  *
  21.  * To do more than just power-of-two sized vectors, see the full
  22.  * version I wrote for NetLib.
  23.  *
  24.  * Note that the packing is a little strange; rather than the FFT r/i
  25.  * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
  26.  * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
  27.  * FORTRAN version
  28.  */
  29.  
  30. #include <stdlib.h>
  31. #include <string.h>
  32. #include <math.h>
  33. #include "smallft.h"
  34. #include "os.h"
  35. #include "misc.h"
  36.  
  37. static void drfti1(int n, float *wa, int *ifac){
  38.   static int ntryh[4] = { 4,2,3,5 };
  39.   static float tpi = 6.28318530717958648f;
  40.   float arg,argh,argld,fi;
  41.   int ntry=0,i,j=-1;
  42.   int k1, l1, l2, ib;
  43.   int ld, ii, ip, is, nq, nr;
  44.   int ido, ipm, nfm1;
  45.   int nl=n;
  46.   int nf=0;
  47.  
  48.  L101:
  49.   j++;
  50.   if (j < 4)
  51.     ntry=ntryh[j];
  52.   else
  53.     ntry+=2;
  54.  
  55.  L104:
  56.   nq=nl/ntry;
  57.   nr=nl-ntry*nq;
  58.   if (nr!=0) goto L101;
  59.  
  60.   nf++;
  61.   ifac[nf+1]=ntry;
  62.   nl=nq;
  63.   if(ntry!=2)goto L107;
  64.   if(nf==1)goto L107;
  65.  
  66.   for (i=1;i<nf;i++){
  67.     ib=nf-i+1;
  68.     ifac[ib+1]=ifac[ib];
  69.   }
  70.   ifac[2] = 2;
  71.  
  72.  L107:
  73.   if(nl!=1)goto L104;
  74.   ifac[0]=n;
  75.   ifac[1]=nf;
  76.   argh=tpi/n;
  77.   is=0;
  78.   nfm1=nf-1;
  79.   l1=1;
  80.  
  81.   if(nfm1==0)return;
  82.  
  83.   for (k1=0;k1<nfm1;k1++){
  84.     ip=ifac[k1+2];
  85.     ld=0;
  86.     l2=l1*ip;
  87.     ido=n/l2;
  88.     ipm=ip-1;
  89.  
  90.     for (j=0;j<ipm;j++){
  91.       ld+=l1;
  92.       i=is;
  93.       argld=(float)ld*argh;
  94.       fi=0.f;
  95.       for (ii=2;ii<ido;ii+=2){
  96.         fi+=1.f;
  97.         arg=fi*argld;
  98.         wa[i++]=cos(arg);
  99.         wa[i++]=sin(arg);
  100.       }
  101.       is+=ido;
  102.     }
  103.     l1=l2;
  104.   }
  105. }
  106.  
  107. static void fdrffti(int n, float *wsave, int *ifac){
  108.  
  109.   if (n == 1) return;
  110.   drfti1(n, wsave+n, ifac);
  111. }
  112.  
  113. static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){
  114.   int i,k;
  115.   float ti2,tr2;
  116.   int t0,t1,t2,t3,t4,t5,t6;
  117.  
  118.   t1=0;
  119.   t0=(t2=l1*ido);
  120.   t3=ido<<1;
  121.   for(k=0;k<l1;k++){
  122.     ch[t1<<1]=cc[t1]+cc[t2];
  123.     ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
  124.     t1+=ido;
  125.     t2+=ido;
  126.   }
  127.  
  128.   if(ido<2)return;
  129.   if(ido==2)goto L105;
  130.  
  131.   t1=0;
  132.   t2=t0;
  133.   for(k=0;k<l1;k++){
  134.     t3=t2;
  135.     t4=(t1<<1)+(ido<<1);
  136.     t5=t1;
  137.     t6=t1+t1;
  138.     for(i=2;i<ido;i+=2){
  139.       t3+=2;
  140.       t4-=2;
  141.       t5+=2;
  142.       t6+=2;
  143.       tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
  144.       ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
  145.       ch[t6]=cc[t5]+ti2;
  146.       ch[t4]=ti2-cc[t5];
  147.       ch[t6-1]=cc[t5-1]+tr2;
  148.       ch[t4-1]=cc[t5-1]-tr2;
  149.     }
  150.     t1+=ido;
  151.     t2+=ido;
  152.   }
  153.  
  154.   if(ido%2==1)return;
  155.  
  156.  L105:
  157.   t3=(t2=(t1=ido)-1);
  158.   t2+=t0;
  159.   for(k=0;k<l1;k++){
  160.     ch[t1]=-cc[t2];
  161.     ch[t1-1]=cc[t3];
  162.     t1+=ido<<1;
  163.     t2+=ido;
  164.     t3+=ido;
  165.   }
  166. }
  167.  
  168. static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,
  169.             float *wa2,float *wa3){
  170.   static float hsqt2 = .70710678118654752f;
  171.   int i,k,t0,t1,t2,t3,t4,t5,t6;
  172.   float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
  173.   t0=l1*ido;
  174.  
  175.   t1=t0;
  176.   t4=t1<<1;
  177.   t2=t1+(t1<<1);
  178.   t3=0;
  179.  
  180.   for(k=0;k<l1;k++){
  181.     tr1=cc[t1]+cc[t2];
  182.     tr2=cc[t3]+cc[t4];
  183.  
  184.     ch[t5=t3<<2]=tr1+tr2;
  185.     ch[(ido<<2)+t5-1]=tr2-tr1;
  186.     ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
  187.     ch[t5]=cc[t2]-cc[t1];
  188.  
  189.     t1+=ido;
  190.     t2+=ido;
  191.     t3+=ido;
  192.     t4+=ido;
  193.   }
  194.  
  195.   if(ido<2)return;
  196.   if(ido==2)goto L105;
  197.  
  198.  
  199.   t1=0;
  200.   for(k=0;k<l1;k++){
  201.     t2=t1;
  202.     t4=t1<<2;
  203.     t5=(t6=ido<<1)+t4;
  204.     for(i=2;i<ido;i+=2){
  205.       t3=(t2+=2);
  206.       t4+=2;
  207.       t5-=2;
  208.  
  209.       t3+=t0;
  210.       cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
  211.       ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
  212.       t3+=t0;
  213.       cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
  214.       ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
  215.       t3+=t0;
  216.       cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
  217.       ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
  218.  
  219.       tr1=cr2+cr4;
  220.       tr4=cr4-cr2;
  221.       ti1=ci2+ci4;
  222.       ti4=ci2-ci4;
  223.  
  224.       ti2=cc[t2]+ci3;
  225.       ti3=cc[t2]-ci3;
  226.       tr2=cc[t2-1]+cr3;
  227.       tr3=cc[t2-1]-cr3;
  228.  
  229.       ch[t4-1]=tr1+tr2;
  230.       ch[t4]=ti1+ti2;
  231.  
  232.       ch[t5-1]=tr3-ti4;
  233.       ch[t5]=tr4-ti3;
  234.  
  235.       ch[t4+t6-1]=ti4+tr3;
  236.       ch[t4+t6]=tr4+ti3;
  237.  
  238.       ch[t5+t6-1]=tr2-tr1;
  239.       ch[t5+t6]=ti1-ti2;
  240.     }
  241.     t1+=ido;
  242.   }
  243.   if(ido&1)return;
  244.  
  245.  L105:
  246.  
  247.   t2=(t1=t0+ido-1)+(t0<<1);
  248.   t3=ido<<2;
  249.   t4=ido;
  250.   t5=ido<<1;
  251.   t6=ido;
  252.  
  253.   for(k=0;k<l1;k++){
  254.     ti1=-hsqt2*(cc[t1]+cc[t2]);
  255.     tr1=hsqt2*(cc[t1]-cc[t2]);
  256.  
  257.     ch[t4-1]=tr1+cc[t6-1];
  258.     ch[t4+t5-1]=cc[t6-1]-tr1;
  259.  
  260.     ch[t4]=ti1-cc[t1+t0];
  261.     ch[t4+t5]=ti1+cc[t1+t0];
  262.  
  263.     t1+=ido;
  264.     t2+=ido;
  265.     t4+=t3;
  266.     t6+=ido;
  267.   }
  268. }
  269.  
  270. static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
  271.                           float *c2,float *ch,float *ch2,float *wa){
  272.  
  273.   static float tpi=6.283185307179586f;
  274.   int idij,ipph,i,j,k,l,ic,ik,is;
  275.   int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
  276.   float dc2,ai1,ai2,ar1,ar2,ds2;
  277.   int nbd;
  278.   float dcp,arg,dsp,ar1h,ar2h;
  279.   int idp2,ipp2;
  280.  
  281.   arg=tpi/(float)ip;
  282.   dcp=cos(arg);
  283.   dsp=sin(arg);
  284.   ipph=(ip+1)>>1;
  285.   ipp2=ip;
  286.   idp2=ido;
  287.   nbd=(ido-1)>>1;
  288.   t0=l1*ido;
  289.   t10=ip*ido;
  290.  
  291.   if(ido==1)goto L119;
  292.   for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
  293.  
  294.   t1=0;
  295.   for(j=1;j<ip;j++){
  296.     t1+=t0;
  297.     t2=t1;
  298.     for(k=0;k<l1;k++){
  299.       ch[t2]=c1[t2];
  300.       t2+=ido;
  301.     }
  302.   }
  303.  
  304.   is=-ido;
  305.   t1=0;
  306.   if(nbd>l1){
  307.     for(j=1;j<ip;j++){
  308.       t1+=t0;
  309.       is+=ido;
  310.       t2= -ido+t1;
  311.       for(k=0;k<l1;k++){
  312.         idij=is-1;
  313.         t2+=ido;
  314.         t3=t2;
  315.         for(i=2;i<ido;i+=2){
  316.           idij+=2;
  317.           t3+=2;
  318.           ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
  319.           ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
  320.         }
  321.       }
  322.     }
  323.   }else{
  324.  
  325.     for(j=1;j<ip;j++){
  326.       is+=ido;
  327.       idij=is-1;
  328.       t1+=t0;
  329.       t2=t1;
  330.       for(i=2;i<ido;i+=2){
  331.         idij+=2;
  332.         t2+=2;
  333.         t3=t2;
  334.         for(k=0;k<l1;k++){
  335.           ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
  336.           ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
  337.           t3+=ido;
  338.         }
  339.       }
  340.     }
  341.   }
  342.  
  343.   t1=0;
  344.   t2=ipp2*t0;
  345.   if(nbd<l1){
  346.     for(j=1;j<ipph;j++){
  347.       t1+=t0;
  348.       t2-=t0;
  349.       t3=t1;
  350.       t4=t2;
  351.       for(i=2;i<ido;i+=2){
  352.         t3+=2;
  353.         t4+=2;
  354.         t5=t3-ido;
  355.         t6=t4-ido;
  356.         for(k=0;k<l1;k++){
  357.           t5+=ido;
  358.           t6+=ido;
  359.           c1[t5-1]=ch[t5-1]+ch[t6-1];
  360.           c1[t6-1]=ch[t5]-ch[t6];
  361.           c1[t5]=ch[t5]+ch[t6];
  362.           c1[t6]=ch[t6-1]-ch[t5-1];
  363.         }
  364.       }
  365.     }
  366.   }else{
  367.     for(j=1;j<ipph;j++){
  368.       t1+=t0;
  369.       t2-=t0;
  370.       t3=t1;
  371.       t4=t2;
  372.       for(k=0;k<l1;k++){
  373.         t5=t3;
  374.         t6=t4;
  375.         for(i=2;i<ido;i+=2){
  376.           t5+=2;
  377.           t6+=2;
  378.           c1[t5-1]=ch[t5-1]+ch[t6-1];
  379.           c1[t6-1]=ch[t5]-ch[t6];
  380.           c1[t5]=ch[t5]+ch[t6];
  381.           c1[t6]=ch[t6-1]-ch[t5-1];
  382.         }
  383.         t3+=ido;
  384.         t4+=ido;
  385.       }
  386.     }
  387.   }
  388.  
  389. L119:
  390.   for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
  391.  
  392.   t1=0;
  393.   t2=ipp2*idl1;
  394.   for(j=1;j<ipph;j++){
  395.     t1+=t0;
  396.     t2-=t0;
  397.     t3=t1-ido;
  398.     t4=t2-ido;
  399.     for(k=0;k<l1;k++){
  400.       t3+=ido;
  401.       t4+=ido;
  402.       c1[t3]=ch[t3]+ch[t4];
  403.       c1[t4]=ch[t4]-ch[t3];
  404.     }
  405.   }
  406.  
  407.   ar1=1.f;
  408.   ai1=0.f;
  409.   t1=0;
  410.   t2=ipp2*idl1;
  411.   t3=(ip-1)*idl1;
  412.   for(l=1;l<ipph;l++){
  413.     t1+=idl1;
  414.     t2-=idl1;
  415.     ar1h=dcp*ar1-dsp*ai1;
  416.     ai1=dcp*ai1+dsp*ar1;
  417.     ar1=ar1h;
  418.     t4=t1;
  419.     t5=t2;
  420.     t6=t3;
  421.     t7=idl1;
  422.  
  423.     for(ik=0;ik<idl1;ik++){
  424.       ch2[t4++]=c2[ik]+ar1*c2[t7++];
  425.       ch2[t5++]=ai1*c2[t6++];
  426.     }
  427.  
  428.     dc2=ar1;
  429.     ds2=ai1;
  430.     ar2=ar1;
  431.     ai2=ai1;
  432.  
  433.     t4=idl1;
  434.     t5=(ipp2-1)*idl1;
  435.     for(j=2;j<ipph;j++){
  436.       t4+=idl1;
  437.       t5-=idl1;
  438.  
  439.       ar2h=dc2*ar2-ds2*ai2;
  440.       ai2=dc2*ai2+ds2*ar2;
  441.       ar2=ar2h;
  442.  
  443.       t6=t1;
  444.       t7=t2;
  445.       t8=t4;
  446.       t9=t5;
  447.       for(ik=0;ik<idl1;ik++){
  448.         ch2[t6++]+=ar2*c2[t8++];
  449.         ch2[t7++]+=ai2*c2[t9++];
  450.       }
  451.     }
  452.   }
  453.  
  454.   t1=0;
  455.   for(j=1;j<ipph;j++){
  456.     t1+=idl1;
  457.     t2=t1;
  458.     for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
  459.   }
  460.  
  461.   if(ido<l1)goto L132;
  462.  
  463.   t1=0;
  464.   t2=0;
  465.   for(k=0;k<l1;k++){
  466.     t3=t1;
  467.     t4=t2;
  468.     for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
  469.     t1+=ido;
  470.     t2+=t10;
  471.   }
  472.  
  473.   goto L135;
  474.  
  475.  L132:
  476.   for(i=0;i<ido;i++){
  477.     t1=i;
  478.     t2=i;
  479.     for(k=0;k<l1;k++){
  480.       cc[t2]=ch[t1];
  481.       t1+=ido;
  482.       t2+=t10;
  483.     }
  484.   }
  485.  
  486.  L135:
  487.   t1=0;
  488.   t2=ido<<1;
  489.   t3=0;
  490.   t4=ipp2*t0;
  491.   for(j=1;j<ipph;j++){
  492.  
  493.     t1+=t2;
  494.     t3+=t0;
  495.     t4-=t0;
  496.  
  497.     t5=t1;
  498.     t6=t3;
  499.     t7=t4;
  500.  
  501.     for(k=0;k<l1;k++){
  502.       cc[t5-1]=ch[t6];
  503.       cc[t5]=ch[t7];
  504.       t5+=t10;
  505.       t6+=ido;
  506.       t7+=ido;
  507.     }
  508.   }
  509.  
  510.   if(ido==1)return;
  511.   if(nbd<l1)goto L141;
  512.  
  513.   t1=-ido;
  514.   t3=0;
  515.   t4=0;
  516.   t5=ipp2*t0;
  517.   for(j=1;j<ipph;j++){
  518.     t1+=t2;
  519.     t3+=t2;
  520.     t4+=t0;
  521.     t5-=t0;
  522.     t6=t1;
  523.     t7=t3;
  524.     t8=t4;
  525.     t9=t5;
  526.     for(k=0;k<l1;k++){
  527.       for(i=2;i<ido;i+=2){
  528.         ic=idp2-i;
  529.         cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
  530.         cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
  531.         cc[i+t7]=ch[i+t8]+ch[i+t9];
  532.         cc[ic+t6]=ch[i+t9]-ch[i+t8];
  533.       }
  534.       t6+=t10;
  535.       t7+=t10;
  536.       t8+=ido;
  537.       t9+=ido;
  538.     }
  539.   }
  540.   return;
  541.  
  542.  L141:
  543.  
  544.   t1=-ido;
  545.   t3=0;
  546.   t4=0;
  547.   t5=ipp2*t0;
  548.   for(j=1;j<ipph;j++){
  549.     t1+=t2;
  550.     t3+=t2;
  551.     t4+=t0;
  552.     t5-=t0;
  553.     for(i=2;i<ido;i+=2){
  554.       t6=idp2+t1-i;
  555.       t7=i+t3;
  556.       t8=i+t4;
  557.       t9=i+t5;
  558.       for(k=0;k<l1;k++){
  559.         cc[t7-1]=ch[t8-1]+ch[t9-1];
  560.         cc[t6-1]=ch[t8-1]-ch[t9-1];
  561.         cc[t7]=ch[t8]+ch[t9];
  562.         cc[t6]=ch[t9]-ch[t8];
  563.         t6+=t10;
  564.         t7+=t10;
  565.         t8+=ido;
  566.         t9+=ido;
  567.       }
  568.     }
  569.   }
  570. }
  571.  
  572. static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){
  573.   int i,k1,l1,l2;
  574.   int na,kh,nf;
  575.   int ip,iw,ido,idl1,ix2,ix3;
  576.  
  577.   nf=ifac[1];
  578.   na=1;
  579.   l2=n;
  580.   iw=n;
  581.  
  582.   for(k1=0;k1<nf;k1++){
  583.     kh=nf-k1;
  584.     ip=ifac[kh+1];
  585.     l1=l2/ip;
  586.     ido=n/l2;
  587.     idl1=ido*l1;
  588.     iw-=(ip-1)*ido;
  589.     na=1-na;
  590.  
  591.     if(ip!=4)goto L102;
  592.  
  593.     ix2=iw+ido;
  594.     ix3=ix2+ido;
  595.     if(na!=0)
  596.       dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
  597.     else
  598.       dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
  599.     goto L110;
  600.  
  601.  L102:
  602.     if(ip!=2)goto L104;
  603.     if(na!=0)goto L103;
  604.  
  605.     dradf2(ido,l1,c,ch,wa+iw-1);
  606.     goto L110;
  607.  
  608.   L103:
  609.     dradf2(ido,l1,ch,c,wa+iw-1);
  610.     goto L110;
  611.  
  612.   L104:
  613.     if(ido==1)na=1-na;
  614.     if(na!=0)goto L109;
  615.  
  616.     dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
  617.     na=1;
  618.     goto L110;
  619.  
  620.   L109:
  621.     dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
  622.     na=0;
  623.  
  624.   L110:
  625.     l2=l1;
  626.   }
  627.  
  628.   if(na==1)return;
  629.  
  630.   for(i=0;i<n;i++)c[i]=ch[i];
  631. }
  632.  
  633. static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){
  634.   int i,k,t0,t1,t2,t3,t4,t5,t6;
  635.   float ti2,tr2;
  636.  
  637.   t0=l1*ido;
  638.  
  639.   t1=0;
  640.   t2=0;
  641.   t3=(ido<<1)-1;
  642.   for(k=0;k<l1;k++){
  643.     ch[t1]=cc[t2]+cc[t3+t2];
  644.     ch[t1+t0]=cc[t2]-cc[t3+t2];
  645.     t2=(t1+=ido)<<1;
  646.   }
  647.  
  648.   if(ido<2)return;
  649.   if(ido==2)goto L105;
  650.  
  651.   t1=0;
  652.   t2=0;
  653.   for(k=0;k<l1;k++){
  654.     t3=t1;
  655.     t5=(t4=t2)+(ido<<1);
  656.     t6=t0+t1;
  657.     for(i=2;i<ido;i+=2){
  658.       t3+=2;
  659.       t4+=2;
  660.       t5-=2;
  661.       t6+=2;
  662.       ch[t3-1]=cc[t4-1]+cc[t5-1];
  663.       tr2=cc[t4-1]-cc[t5-1];
  664.       ch[t3]=cc[t4]-cc[t5];
  665.       ti2=cc[t4]+cc[t5];
  666.       ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
  667.       ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
  668.     }
  669.     t2=(t1+=ido)<<1;
  670.   }
  671.  
  672.   if(ido%2==1)return;
  673.  
  674. L105:
  675.   t1=ido-1;
  676.   t2=ido-1;
  677.   for(k=0;k<l1;k++){
  678.     ch[t1]=cc[t2]+cc[t2];
  679.     ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
  680.     t1+=ido;
  681.     t2+=ido<<1;
  682.   }
  683. }
  684.  
  685. static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1,
  686.                           float *wa2){
  687.   static float taur = -.5f;
  688.   static float taui = .8660254037844386f;
  689.   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
  690.   float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
  691.   t0=l1*ido;
  692.  
  693.   t1=0;
  694.   t2=t0<<1;
  695.   t3=ido<<1;
  696.   t4=ido+(ido<<1);
  697.   t5=0;
  698.   for(k=0;k<l1;k++){
  699.     tr2=cc[t3-1]+cc[t3-1];
  700.     cr2=cc[t5]+(taur*tr2);
  701.     ch[t1]=cc[t5]+tr2;
  702.     ci3=taui*(cc[t3]+cc[t3]);
  703.     ch[t1+t0]=cr2-ci3;
  704.     ch[t1+t2]=cr2+ci3;
  705.     t1+=ido;
  706.     t3+=t4;
  707.     t5+=t4;
  708.   }
  709.  
  710.   if(ido==1)return;
  711.  
  712.   t1=0;
  713.   t3=ido<<1;
  714.   for(k=0;k<l1;k++){
  715.     t7=t1+(t1<<1);
  716.     t6=(t5=t7+t3);
  717.     t8=t1;
  718.     t10=(t9=t1+t0)+t0;
  719.  
  720.     for(i=2;i<ido;i+=2){
  721.       t5+=2;
  722.       t6-=2;
  723.       t7+=2;
  724.       t8+=2;
  725.       t9+=2;
  726.       t10+=2;
  727.       tr2=cc[t5-1]+cc[t6-1];
  728.       cr2=cc[t7-1]+(taur*tr2);
  729.       ch[t8-1]=cc[t7-1]+tr2;
  730.       ti2=cc[t5]-cc[t6];
  731.       ci2=cc[t7]+(taur*ti2);
  732.       ch[t8]=cc[t7]+ti2;
  733.       cr3=taui*(cc[t5-1]-cc[t6-1]);
  734.       ci3=taui*(cc[t5]+cc[t6]);
  735.       dr2=cr2-ci3;
  736.       dr3=cr2+ci3;
  737.       di2=ci2+cr3;
  738.       di3=ci2-cr3;
  739.       ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
  740.       ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
  741.       ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
  742.       ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
  743.     }
  744.     t1+=ido;
  745.   }
  746. }
  747.  
  748. static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,
  749.                           float *wa2,float *wa3){
  750.   static float sqrt2=1.414213562373095f;
  751.   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
  752.   float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
  753.   t0=l1*ido;
  754.  
  755.   t1=0;
  756.   t2=ido<<2;
  757.   t3=0;
  758.   t6=ido<<1;
  759.   for(k=0;k<l1;k++){
  760.     t4=t3+t6;
  761.     t5=t1;
  762.     tr3=cc[t4-1]+cc[t4-1];
  763.     tr4=cc[t4]+cc[t4];
  764.     tr1=cc[t3]-cc[(t4+=t6)-1];
  765.     tr2=cc[t3]+cc[t4-1];
  766.     ch[t5]=tr2+tr3;
  767.     ch[t5+=t0]=tr1-tr4;
  768.     ch[t5+=t0]=tr2-tr3;
  769.     ch[t5+=t0]=tr1+tr4;
  770.     t1+=ido;
  771.     t3+=t2;
  772.   }
  773.  
  774.   if(ido<2)return;
  775.   if(ido==2)goto L105;
  776.  
  777.   t1=0;
  778.   for(k=0;k<l1;k++){
  779.     t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
  780.     t7=t1;
  781.     for(i=2;i<ido;i+=2){
  782.       t2+=2;
  783.       t3+=2;
  784.       t4-=2;
  785.       t5-=2;
  786.       t7+=2;
  787.       ti1=cc[t2]+cc[t5];
  788.       ti2=cc[t2]-cc[t5];
  789.       ti3=cc[t3]-cc[t4];
  790.       tr4=cc[t3]+cc[t4];
  791.       tr1=cc[t2-1]-cc[t5-1];
  792.       tr2=cc[t2-1]+cc[t5-1];
  793.       ti4=cc[t3-1]-cc[t4-1];
  794.       tr3=cc[t3-1]+cc[t4-1];
  795.       ch[t7-1]=tr2+tr3;
  796.       cr3=tr2-tr3;
  797.       ch[t7]=ti2+ti3;
  798.       ci3=ti2-ti3;
  799.       cr2=tr1-tr4;
  800.       cr4=tr1+tr4;
  801.       ci2=ti1+ti4;
  802.       ci4=ti1-ti4;
  803.  
  804.       ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
  805.       ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
  806.       ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
  807.       ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
  808.       ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
  809.       ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
  810.     }
  811.     t1+=ido;
  812.   }
  813.  
  814.   if(ido%2 == 1)return;
  815.  
  816.  L105:
  817.  
  818.   t1=ido;
  819.   t2=ido<<2;
  820.   t3=ido-1;
  821.   t4=ido+(ido<<1);
  822.   for(k=0;k<l1;k++){
  823.     t5=t3;
  824.     ti1=cc[t1]+cc[t4];
  825.     ti2=cc[t4]-cc[t1];
  826.     tr1=cc[t1-1]-cc[t4-1];
  827.     tr2=cc[t1-1]+cc[t4-1];
  828.     ch[t5]=tr2+tr2;
  829.     ch[t5+=t0]=sqrt2*(tr1-ti1);
  830.     ch[t5+=t0]=ti2+ti2;
  831.     ch[t5+=t0]=-sqrt2*(tr1+ti1);
  832.  
  833.     t3+=ido;
  834.     t1+=t2;
  835.     t4+=t2;
  836.   }
  837. }
  838.  
  839. static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
  840.             float *c2,float *ch,float *ch2,float *wa){
  841.   static float tpi=6.283185307179586f;
  842.   int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
  843.       t11,t12;
  844.   float dc2,ai1,ai2,ar1,ar2,ds2;
  845.   int nbd;
  846.   float dcp,arg,dsp,ar1h,ar2h;
  847.   int ipp2;
  848.  
  849.   t10=ip*ido;
  850.   t0=l1*ido;
  851.   arg=tpi/(float)ip;
  852.   dcp=cos(arg);
  853.   dsp=sin(arg);
  854.   nbd=(ido-1)>>1;
  855.   ipp2=ip;
  856.   ipph=(ip+1)>>1;
  857.   if(ido<l1)goto L103;
  858.  
  859.   t1=0;
  860.   t2=0;
  861.   for(k=0;k<l1;k++){
  862.     t3=t1;
  863.     t4=t2;
  864.     for(i=0;i<ido;i++){
  865.       ch[t3]=cc[t4];
  866.       t3++;
  867.       t4++;
  868.     }
  869.     t1+=ido;
  870.     t2+=t10;
  871.   }
  872.   goto L106;
  873.  
  874.  L103:
  875.   t1=0;
  876.   for(i=0;i<ido;i++){
  877.     t2=t1;
  878.     t3=t1;
  879.     for(k=0;k<l1;k++){
  880.       ch[t2]=cc[t3];
  881.       t2+=ido;
  882.       t3+=t10;
  883.     }
  884.     t1++;
  885.   }
  886.  
  887.  L106:
  888.   t1=0;
  889.   t2=ipp2*t0;
  890.   t7=(t5=ido<<1);
  891.   for(j=1;j<ipph;j++){
  892.     t1+=t0;
  893.     t2-=t0;
  894.     t3=t1;
  895.     t4=t2;
  896.     t6=t5;
  897.     for(k=0;k<l1;k++){
  898.       ch[t3]=cc[t6-1]+cc[t6-1];
  899.       ch[t4]=cc[t6]+cc[t6];
  900.       t3+=ido;
  901.       t4+=ido;
  902.       t6+=t10;
  903.     }
  904.     t5+=t7;
  905.   }
  906.  
  907.   if (ido == 1)goto L116;
  908.   if(nbd<l1)goto L112;
  909.  
  910.   t1=0;
  911.   t2=ipp2*t0;
  912.   t7=0;
  913.   for(j=1;j<ipph;j++){
  914.     t1+=t0;
  915.     t2-=t0;
  916.     t3=t1;
  917.     t4=t2;
  918.  
  919.     t7+=(ido<<1);
  920.     t8=t7;
  921.     for(k=0;k<l1;k++){
  922.       t5=t3;
  923.       t6=t4;
  924.       t9=t8;
  925.       t11=t8;
  926.       for(i=2;i<ido;i+=2){
  927.         t5+=2;
  928.         t6+=2;
  929.         t9+=2;
  930.         t11-=2;
  931.         ch[t5-1]=cc[t9-1]+cc[t11-1];
  932.         ch[t6-1]=cc[t9-1]-cc[t11-1];
  933.         ch[t5]=cc[t9]-cc[t11];
  934.         ch[t6]=cc[t9]+cc[t11];
  935.       }
  936.       t3+=ido;
  937.       t4+=ido;
  938.       t8+=t10;
  939.     }
  940.   }
  941.   goto L116;
  942.  
  943.  L112:
  944.   t1=0;
  945.   t2=ipp2*t0;
  946.   t7=0;
  947.   for(j=1;j<ipph;j++){
  948.     t1+=t0;
  949.     t2-=t0;
  950.     t3=t1;
  951.     t4=t2;
  952.     t7+=(ido<<1);
  953.     t8=t7;
  954.     t9=t7;
  955.     for(i=2;i<ido;i+=2){
  956.       t3+=2;
  957.       t4+=2;
  958.       t8+=2;
  959.       t9-=2;
  960.       t5=t3;
  961.       t6=t4;
  962.       t11=t8;
  963.       t12=t9;
  964.       for(k=0;k<l1;k++){
  965.         ch[t5-1]=cc[t11-1]+cc[t12-1];
  966.         ch[t6-1]=cc[t11-1]-cc[t12-1];
  967.         ch[t5]=cc[t11]-cc[t12];
  968.         ch[t6]=cc[t11]+cc[t12];
  969.         t5+=ido;
  970.         t6+=ido;
  971.         t11+=t10;
  972.         t12+=t10;
  973.       }
  974.     }
  975.   }
  976.  
  977. L116:
  978.   ar1=1.f;
  979.   ai1=0.f;
  980.   t1=0;
  981.   t9=(t2=ipp2*idl1);
  982.   t3=(ip-1)*idl1;
  983.   for(l=1;l<ipph;l++){
  984.     t1+=idl1;
  985.     t2-=idl1;
  986.  
  987.     ar1h=dcp*ar1-dsp*ai1;
  988.     ai1=dcp*ai1+dsp*ar1;
  989.     ar1=ar1h;
  990.     t4=t1;
  991.     t5=t2;
  992.     t6=0;
  993.     t7=idl1;
  994.     t8=t3;
  995.     for(ik=0;ik<idl1;ik++){
  996.       c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
  997.       c2[t5++]=ai1*ch2[t8++];
  998.     }
  999.     dc2=ar1;
  1000.     ds2=ai1;
  1001.     ar2=ar1;
  1002.     ai2=ai1;
  1003.  
  1004.     t6=idl1;
  1005.     t7=t9-idl1;
  1006.     for(j=2;j<ipph;j++){
  1007.       t6+=idl1;
  1008.       t7-=idl1;
  1009.       ar2h=dc2*ar2-ds2*ai2;
  1010.       ai2=dc2*ai2+ds2*ar2;
  1011.       ar2=ar2h;
  1012.       t4=t1;
  1013.       t5=t2;
  1014.       t11=t6;
  1015.       t12=t7;
  1016.       for(ik=0;ik<idl1;ik++){
  1017.         c2[t4++]+=ar2*ch2[t11++];
  1018.         c2[t5++]+=ai2*ch2[t12++];
  1019.       }
  1020.     }
  1021.   }
  1022.  
  1023.   t1=0;
  1024.   for(j=1;j<ipph;j++){
  1025.     t1+=idl1;
  1026.     t2=t1;
  1027.     for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
  1028.   }
  1029.  
  1030.   t1=0;
  1031.   t2=ipp2*t0;
  1032.   for(j=1;j<ipph;j++){
  1033.     t1+=t0;
  1034.     t2-=t0;
  1035.     t3=t1;
  1036.     t4=t2;
  1037.     for(k=0;k<l1;k++){
  1038.       ch[t3]=c1[t3]-c1[t4];
  1039.       ch[t4]=c1[t3]+c1[t4];
  1040.       t3+=ido;
  1041.       t4+=ido;
  1042.     }
  1043.   }
  1044.  
  1045.   if(ido==1)goto L132;
  1046.   if(nbd<l1)goto L128;
  1047.  
  1048.   t1=0;
  1049.   t2=ipp2*t0;
  1050.   for(j=1;j<ipph;j++){
  1051.     t1+=t0;
  1052.     t2-=t0;
  1053.     t3=t1;
  1054.     t4=t2;
  1055.     for(k=0;k<l1;k++){
  1056.       t5=t3;
  1057.       t6=t4;
  1058.       for(i=2;i<ido;i+=2){
  1059.         t5+=2;
  1060.         t6+=2;
  1061.         ch[t5-1]=c1[t5-1]-c1[t6];
  1062.         ch[t6-1]=c1[t5-1]+c1[t6];
  1063.         ch[t5]=c1[t5]+c1[t6-1];
  1064.         ch[t6]=c1[t5]-c1[t6-1];
  1065.       }
  1066.       t3+=ido;
  1067.       t4+=ido;
  1068.     }
  1069.   }
  1070.   goto L132;
  1071.  
  1072.  L128:
  1073.   t1=0;
  1074.   t2=ipp2*t0;
  1075.   for(j=1;j<ipph;j++){
  1076.     t1+=t0;
  1077.     t2-=t0;
  1078.     t3=t1;
  1079.     t4=t2;
  1080.     for(i=2;i<ido;i+=2){
  1081.       t3+=2;
  1082.       t4+=2;
  1083.       t5=t3;
  1084.       t6=t4;
  1085.       for(k=0;k<l1;k++){
  1086.         ch[t5-1]=c1[t5-1]-c1[t6];
  1087.         ch[t6-1]=c1[t5-1]+c1[t6];
  1088.         ch[t5]=c1[t5]+c1[t6-1];
  1089.         ch[t6]=c1[t5]-c1[t6-1];
  1090.         t5+=ido;
  1091.         t6+=ido;
  1092.       }
  1093.     }
  1094.   }
  1095.  
  1096. L132:
  1097.   if(ido==1)return;
  1098.  
  1099.   for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
  1100.  
  1101.   t1=0;
  1102.   for(j=1;j<ip;j++){
  1103.     t2=(t1+=t0);
  1104.     for(k=0;k<l1;k++){
  1105.       c1[t2]=ch[t2];
  1106.       t2+=ido;
  1107.     }
  1108.   }
  1109.  
  1110.   if(nbd>l1)goto L139;
  1111.  
  1112.   is= -ido-1;
  1113.   t1=0;
  1114.   for(j=1;j<ip;j++){
  1115.     is+=ido;
  1116.     t1+=t0;
  1117.     idij=is;
  1118.     t2=t1;
  1119.     for(i=2;i<ido;i+=2){
  1120.       t2+=2;
  1121.       idij+=2;
  1122.       t3=t2;
  1123.       for(k=0;k<l1;k++){
  1124.         c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
  1125.         c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
  1126.         t3+=ido;
  1127.       }
  1128.     }
  1129.   }
  1130.   return;
  1131.  
  1132.  L139:
  1133.   is= -ido-1;
  1134.   t1=0;
  1135.   for(j=1;j<ip;j++){
  1136.     is+=ido;
  1137.     t1+=t0;
  1138.     t2=t1;
  1139.     for(k=0;k<l1;k++){
  1140.       idij=is;
  1141.       t3=t2;
  1142.       for(i=2;i<ido;i+=2){
  1143.         idij+=2;
  1144.         t3+=2;
  1145.         c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
  1146.         c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
  1147.       }
  1148.       t2+=ido;
  1149.     }
  1150.   }
  1151. }
  1152.  
  1153. static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){
  1154.   int i,k1,l1,l2;
  1155.   int na;
  1156.   int nf,ip,iw,ix2,ix3,ido,idl1;
  1157.  
  1158.   nf=ifac[1];
  1159.   na=0;
  1160.   l1=1;
  1161.   iw=1;
  1162.  
  1163.   for(k1=0;k1<nf;k1++){
  1164.     ip=ifac[k1 + 2];
  1165.     l2=ip*l1;
  1166.     ido=n/l2;
  1167.     idl1=ido*l1;
  1168.     if(ip!=4)goto L103;
  1169.     ix2=iw+ido;
  1170.     ix3=ix2+ido;
  1171.  
  1172.     if(na!=0)
  1173.       dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
  1174.     else
  1175.       dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
  1176.     na=1-na;
  1177.     goto L115;
  1178.  
  1179.   L103:
  1180.     if(ip!=2)goto L106;
  1181.  
  1182.     if(na!=0)
  1183.       dradb2(ido,l1,ch,c,wa+iw-1);
  1184.     else
  1185.       dradb2(ido,l1,c,ch,wa+iw-1);
  1186.     na=1-na;
  1187.     goto L115;
  1188.  
  1189.   L106:
  1190.     if(ip!=3)goto L109;
  1191.  
  1192.     ix2=iw+ido;
  1193.     if(na!=0)
  1194.       dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
  1195.     else
  1196.       dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
  1197.     na=1-na;
  1198.     goto L115;
  1199.  
  1200.   L109:
  1201. /*    The radix five case can be translated later..... */
  1202. /*    if(ip!=5)goto L112;
  1203.  
  1204.     ix2=iw+ido;
  1205.     ix3=ix2+ido;
  1206.     ix4=ix3+ido;
  1207.     if(na!=0)
  1208.       dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
  1209.     else
  1210.       dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
  1211.     na=1-na;
  1212.     goto L115;
  1213.  
  1214.   L112:*/
  1215.     if(na!=0)
  1216.       dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
  1217.     else
  1218.       dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
  1219.     if(ido==1)na=1-na;
  1220.  
  1221.   L115:
  1222.     l1=l2;
  1223.     iw+=(ip-1)*ido;
  1224.   }
  1225.  
  1226.   if(na==0)return;
  1227.  
  1228.   for(i=0;i<n;i++)c[i]=ch[i];
  1229. }
  1230.  
  1231. void drft_forward(drft_lookup *l,float *data){
  1232.   if(l->n==1)return;
  1233.   drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
  1234. }
  1235.  
  1236. void drft_backward(drft_lookup *l,float *data){
  1237.   if (l->n==1)return;
  1238.   drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
  1239. }
  1240.  
  1241. void drft_init(drft_lookup *l,int n){
  1242.   l->n=n;
  1243.   l->trigcache=_ogg_calloc(3*n,sizeof(*l->trigcache));
  1244.   l->splitcache=_ogg_calloc(32,sizeof(*l->splitcache));
  1245.   fdrffti(n, l->trigcache, l->splitcache);
  1246. }
  1247.  
  1248. void drft_clear(drft_lookup *l){
  1249.   if(l){
  1250.     if(l->trigcache)_ogg_free(l->trigcache);
  1251.     if(l->splitcache)_ogg_free(l->splitcache);
  1252.     memset(l,0,sizeof(*l));
  1253.   }
  1254. }
  1255.