12a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)/******************************************************************** 22a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) * * 32a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * 42a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * 52a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * 62a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * 72a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) * * 85d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * 9c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) * by the XIPHOPHORUS Company http://www.xiph.org/ * 10f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) * * 112a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ******************************************************************** 12f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) 132a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) function: *unnormalized* fft transform 142a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) last mod: $Id: smallft.c,v 1.19 2003/10/08 05:12:37 jm Exp $ 152a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 162a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ********************************************************************/ 172a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 182a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)/* FFT implementation from OggSquish, minus cosine transforms, 19a3f6a49ab37290eeeb8db0f41ec0f1cb74a68be7Torne (Richard Coles) * minus all but radix 2/4 case. In Vorbis we only need this 20c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) * cut-down version. 212a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) * 2268043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) * To do more than just power-of-two sized vectors, see the full 23a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) * version I wrote for NetLib. 244e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) * 252a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) * Note that the packing is a little strange; rather than the FFT r/i 26a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, 272a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the 28868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) * FORTRAN version 292a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) */ 302a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 31a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles)#ifdef HAVE_CONFIG_H 322a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)#include "config.h" 332a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)#endif 342a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 352a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)#include <math.h> 362a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)#include "smallft.h" 372a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)#include "arch.h" 382a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)#include "os_support.h" 392a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 402a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)static void drfti1(int n, float *wa, int *ifac){ 412a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) static int ntryh[4] = { 4,2,3,5 }; 42b2df76ea8fec9e32f6f3718986dba0d95315b29cTorne (Richard Coles) static float tpi = 6.28318530717958648f; 432a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) float arg,argh,argld,fi; 44f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) int ntry=0,i,j=-1; 45a3f6a49ab37290eeeb8db0f41ec0f1cb74a68be7Torne (Richard Coles) int k1, l1, l2, ib; 462a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) int ld, ii, ip, is, nq, nr; 472a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) int ido, ipm, nfm1; 482a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) int nl=n; 492a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) int nf=0; 502a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 512a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) L101: 522a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) j++; 537d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) if (j < 4) 542a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ntry=ntryh[j]; 552a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) else 567d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) ntry+=2; 572a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 582a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) L104: 592a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) nq=nl/ntry; 602a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) nr=nl-ntry*nq; 612a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) if (nr!=0) goto L101; 624e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) 634e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) nf++; 642a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ifac[nf+1]=ntry; 655d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) nl=nq; 665d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) if(ntry!=2)goto L107; 672a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) if(nf==1)goto L107; 682a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 692a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for (i=1;i<nf;i++){ 702a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ib=nf-i+1; 717d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) ifac[ib+1]=ifac[ib]; 727d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) } 737d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) ifac[2] = 2; 747d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) 757d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) L107: 762a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) if(nl!=1)goto L104; 772a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ifac[0]=n; 782a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ifac[1]=nf; 792a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) argh=tpi/n; 802a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) is=0; 812a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) nfm1=nf-1; 822a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) l1=1; 832a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 845d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) if(nfm1==0)return; 852a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 862a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for (k1=0;k1<nfm1;k1++){ 872a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ip=ifac[k1+2]; 88c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ld=0; 89c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) l2=l1*ip; 90c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ido=n/l2; 91c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ipm=ip-1; 92c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 932a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for (j=0;j<ipm;j++){ 942a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ld+=l1; 95c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) i=is; 9668043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) argld=(float)ld*argh; 9768043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) fi=0.f; 981e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) for (ii=2;ii<ido;ii+=2){ 99c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch fi+=1.f; 100c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) arg=fi*argld; 101c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) wa[i++]=cos(arg); 1022a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) wa[i++]=sin(arg); 1032a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) } 1042a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) is+=ido; 1052a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) } 106c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) l1=l2; 107c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 1082a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)} 1092a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 1102a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)static void fdrffti(int n, float *wsave, int *ifac){ 111c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 112c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) if (n == 1) return; 1132a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) drfti1(n, wsave+n, ifac); 1142a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)} 1152a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 116c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles)static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){ 11768043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) int i,k; 118c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch float ti2,tr2; 119c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) int t0,t1,t2,t3,t4,t5,t6; 1202a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 121c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t1=0; 1222a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t0=(t2=l1*ido); 1232a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t3=ido<<1; 1242a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for(k=0;k<l1;k++){ 1252a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t1<<1]=cc[t1]+cc[t2]; 1262a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[(t1<<1)+t3-1]=cc[t1]-cc[t2]; 1272a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t1+=ido; 12868043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) t2+=ido; 129c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 130c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 131c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch if(ido<2)return; 1322a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) if(ido==2)goto L105; 1332a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 1342a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t1=0; 1352a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t2=t0; 1362a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for(k=0;k<l1;k++){ 1372a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t3=t2; 1382a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t4=(t1<<1)+(ido<<1); 1392a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t5=t1; 1405d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) t6=t1+t1; 1412a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for(i=2;i<ido;i+=2){ 1422a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t3+=2; 1432a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t4-=2; 1442a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t5+=2; 1452a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t6+=2; 1462a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3]; 1472a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1]; 1484e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ch[t6]=cc[t5]+ti2; 1494e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ch[t4]=ti2-cc[t5]; 1504e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ch[t6-1]=cc[t5-1]+tr2; 1512a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t4-1]=cc[t5-1]-tr2; 1522a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) } 1532a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t1+=ido; 1542a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t2+=ido; 1552a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) } 1562a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 1572a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) if(ido%2==1)return; 1582a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 1592a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) L105: 1602a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t3=(t2=(t1=ido)-1); 1612a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t2+=t0; 162c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(k=0;k<l1;k++){ 1632a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t1]=-cc[t2]; 1642a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t1-1]=cc[t3]; 1652a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t1+=ido<<1; 1662a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t2+=ido; 1672a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t3+=ido; 1682a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) } 1692a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles)} 1702a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 171c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles)static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1, 1722a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) float *wa2,float *wa3){ 1732a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) static float hsqt2 = .70710678118654752f; 1742a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) int i,k,t0,t1,t2,t3,t4,t5,t6; 1752a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4; 1762a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t0=l1*ido; 1772a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 1782a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t1=t0; 1797d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t4=t1<<1; 180c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t2=t1+(t1<<1); 181f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) t3=0; 182a3f6a49ab37290eeeb8db0f41ec0f1cb74a68be7Torne (Richard Coles) 1832a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for(k=0;k<l1;k++){ 1842a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) tr1=cc[t1]+cc[t2]; 1852a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) tr2=cc[t3]+cc[t4]; 1862a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 187868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) ch[t5=t3<<2]=tr1+tr2; 188868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) ch[(ido<<2)+t5-1]=tr2-tr1; 189868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4]; 190868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) ch[t5]=cc[t2]-cc[t1]; 191868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) 192868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t1+=ido; 193868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t2+=ido; 194868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t3+=ido; 195868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t4+=ido; 196868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) } 197868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) 198868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) if(ido<2)return; 199868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) if(ido==2)goto L105; 200868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) 201868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) 202868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t1=0; 203868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) for(k=0;k<l1;k++){ 204868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t2=t1; 205868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t4=t1<<2; 206868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t5=(t6=ido<<1)+t4; 207868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) for(i=2;i<ido;i+=2){ 208868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t3=(t2+=2); 209868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t4+=2; 210868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t5-=2; 211868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) 212868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t3+=t0; 213868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3]; 214868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1]; 2154e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t3+=t0; 2164e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3]; 2174e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1]; 2184e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t3+=t0; 2195d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3]; 2204e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1]; 2214e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) 2224e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) tr1=cr2+cr4; 2234e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) tr4=cr4-cr2; 2244e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ti1=ci2+ci4; 2255d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) ti4=ci2-ci4; 2264e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) 2274e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ti2=cc[t2]+ci3; 2284e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ti3=cc[t2]-ci3; 2294e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) tr2=cc[t2-1]+cr3; 2304e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) tr3=cc[t2-1]-cr3; 2312a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 2322a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t4-1]=tr1+tr2; 2332a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t4]=ti1+ti2; 234eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch 235eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch ch[t5-1]=tr3-ti4; 236eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch ch[t5]=tr4-ti3; 237c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 2384e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ch[t4+t6-1]=ti4+tr3; 239c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ch[t4+t6]=tr4+ti3; 2402a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 2412a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t5+t6-1]=tr2-tr1; 2422a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t5+t6]=ti1-ti2; 243eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch } 244eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch t1+=ido; 245eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch } 2462a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) if(ido&1)return; 2472a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 248c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) L105: 2494e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) 250c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2=(t1=t0+ido-1)+(t0<<1); 2512a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t3=ido<<2; 2522a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t4=ido; 2537dbb3d5cf0c15f500944d211057644d6a2f37371Ben Murdoch t5=ido<<1; 2547dbb3d5cf0c15f500944d211057644d6a2f37371Ben Murdoch t6=ido; 2557dbb3d5cf0c15f500944d211057644d6a2f37371Ben Murdoch 2567dbb3d5cf0c15f500944d211057644d6a2f37371Ben Murdoch for(k=0;k<l1;k++){ 2572a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ti1=-hsqt2*(cc[t1]+cc[t2]); 2582a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) tr1=hsqt2*(cc[t1]-cc[t2]); 2592a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 260c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ch[t4-1]=tr1+cc[t6-1]; 2614e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ch[t4+t5-1]=cc[t6-1]-tr1; 262c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 2632a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t4]=ti1-cc[t1+t0]; 2642a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) ch[t4+t5]=ti1+cc[t1+t0]; 2652a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 266eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch t1+=ido; 267eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch t2+=ido; 268eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch t4+=t3; 2692a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t6+=ido; 2702a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) } 271c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles)} 2724e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) 273c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles)static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1, 274c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) float *c2,float *ch,float *ch2,float *wa){ 275c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 2761320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci static float tpi=6.283185307179586f; 2771320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci int idij,ipph,i,j,k,l,ic,ik,is; 2781320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10; 2791320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci float dc2,ai1,ai2,ar1,ar2,ds2; 2801320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci int nbd; 2811320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci float dcp,arg,dsp,ar1h,ar2h; 2821320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci int idp2,ipp2; 2831320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci 2841320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci arg=tpi/(float)ip; 2851320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci dcp=cos(arg); 2861320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci dsp=sin(arg); 2871320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci ipph=(ip+1)>>1; 2881320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci ipp2=ip; 2891320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci idp2=ido; 2901320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci nbd=(ido-1)>>1; 2911320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci t0=l1*ido; 2921320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci t10=ip*ido; 2931320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci 2941320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci if(ido==1)goto L119; 2951320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik]; 2961320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci 2971320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci t1=0; 2981320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci for(j=1;j<ip;j++){ 2991320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci t1+=t0; 3001320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci t2=t1; 3011320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci for(k=0;k<l1;k++){ 3021320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci ch[t2]=c1[t2]; 3031320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci t2+=ido; 3041320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci } 3051320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci } 3061320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci 3071320f92c476a1ad9d19dba2a48c72b75566198e9Primiano Tucci is=-ido; 3087d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t1=0; 3097d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) if(nbd>l1){ 3107d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) for(j=1;j<ip;j++){ 311c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t1+=t0; 312c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) is+=ido; 313f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) t2= -ido+t1; 314c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(k=0;k<l1;k++){ 3154e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) idij=is-1; 316c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2+=ido; 317c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t3=t2; 318c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(i=2;i<ido;i+=2){ 319c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) idij+=2; 3207d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t3+=2; 3217d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3]; 3227d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1]; 3237d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) } 324b2df76ea8fec9e32f6f3718986dba0d95315b29cTorne (Richard Coles) } 325f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) } 326c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) }else{ 3274e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) 328c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(j=1;j<ip;j++){ 329c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) is+=ido; 330c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) idij=is-1; 331c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t1+=t0; 332c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2=t1; 333c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(i=2;i<ido;i+=2){ 334c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) idij+=2; 335c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2+=2; 33668043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) t3=t2; 337c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch for(k=0;k<l1;k++){ 338c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3]; 339c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1]; 340c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t3+=ido; 341c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 342c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 343c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 344c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 34568043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) 346c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t1=0; 347c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2=ipp2*t0; 348c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) if(nbd<l1){ 3497d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) for(j=1;j<ipph;j++){ 3507d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t1+=t0; 3517d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t2-=t0; 3527d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t3=t1; 353b2df76ea8fec9e32f6f3718986dba0d95315b29cTorne (Richard Coles) t4=t2; 354f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) for(i=2;i<ido;i+=2){ 355c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t3+=2; 3564e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t4+=2; 357c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t5=t3-ido; 358c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t6=t4-ido; 359c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(k=0;k<l1;k++){ 360c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t5+=ido; 361c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t6+=ido; 362c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) c1[t5-1]=ch[t5-1]+ch[t6-1]; 363c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) c1[t6-1]=ch[t5]-ch[t6]; 364c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) c1[t5]=ch[t5]+ch[t6]; 36568043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) c1[t6]=ch[t6-1]-ch[t5-1]; 366c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch } 367c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 368c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 3697d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) }else{ 3707d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) for(j=1;j<ipph;j++){ 3717d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t1+=t0; 3727d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t2-=t0; 373b2df76ea8fec9e32f6f3718986dba0d95315b29cTorne (Richard Coles) t3=t1; 374f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) t4=t2; 375c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(k=0;k<l1;k++){ 3764e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t5=t3; 377c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t6=t4; 378c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(i=2;i<ido;i+=2){ 379c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t5+=2; 380c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t6+=2; 381c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) c1[t5-1]=ch[t5-1]+ch[t6-1]; 382c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) c1[t6-1]=ch[t5]-ch[t6]; 383c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) c1[t5]=ch[t5]+ch[t6]; 384c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) c1[t6]=ch[t6-1]-ch[t5-1]; 38568043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) } 386c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t3+=ido; 387c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t4+=ido; 388c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 389c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 390c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 391c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 392c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles)L119: 393c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik]; 394c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 395c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t1=0; 396c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2=ipp2*idl1; 397c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(j=1;j<ipph;j++){ 39868043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) t1+=t0; 399c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t2-=t0; 400c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t3=t1-ido; 401c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t4=t2-ido; 4027d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) for(k=0;k<l1;k++){ 4037d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t3+=ido; 4047d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) t4+=ido; 4057d4cd473f85ac64c3747c96c277f9e506a0d2246Torne (Richard Coles) c1[t3]=ch[t3]+ch[t4]; 406b2df76ea8fec9e32f6f3718986dba0d95315b29cTorne (Richard Coles) c1[t4]=ch[t4]-ch[t3]; 407f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) } 408c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 4094e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) 410c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ar1=1.f; 411c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ai1=0.f; 412c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t1=0; 413c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2=ipp2*idl1; 414c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t3=(ip-1)*idl1; 415c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(l=1;l<ipph;l++){ 416c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t1+=idl1; 417c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2-=idl1; 418c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ar1h=dcp*ar1-dsp*ai1; 419c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ai1=dcp*ai1+dsp*ar1; 420c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ar1=ar1h; 421c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t4=t1; 422c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t5=t2; 423c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t6=t3; 42468043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) t7=idl1; 425c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch 426c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(ik=0;ik<idl1;ik++){ 427c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) ch2[t4++]=c2[ik]+ar1*c2[t7++]; 428a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) ch2[t5++]=ai1*c2[t6++]; 429a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) } 430a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) 431a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) dc2=ar1; 432a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) ds2=ai1; 433a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) ar2=ar1; 434a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) ai2=ai1; 435a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) 436a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t4=idl1; 437a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t5=(ipp2-1)*idl1; 438a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) for(j=2;j<ipph;j++){ 439a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t4+=idl1; 440a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t5-=idl1; 441a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) 442a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) ar2h=dc2*ar2-ds2*ai2; 443a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) ai2=dc2*ai2+ds2*ar2; 444a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) ar2=ar2h; 445a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) 446a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t6=t1; 447a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t7=t2; 448a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t8=t4; 449a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t9=t5; 450a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) for(ik=0;ik<idl1;ik++){ 451a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) ch2[t6++]+=ar2*c2[t8++]; 452a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) ch2[t7++]+=ai2*c2[t9++]; 453a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) } 454a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) } 455a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) } 456a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) 457a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t1=0; 458c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch for(j=1;j<ipph;j++){ 459a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t1+=idl1; 460a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t2=t1; 461a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++]; 462a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) } 463a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) 464c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) if(ido<l1)goto L132; 465c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 4664e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t1=0; 46768043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) t2=0; 468c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch for(k=0;k<l1;k++){ 469c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t3=t1; 470c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t4=t2; 471c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(i=0;i<ido;i++)cc[t4++]=ch[t3++]; 47268043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) t1+=ido; 473c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t2+=t10; 474c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 475c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 476c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) goto L135; 477c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) 478c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) L132: 47968043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) for(i=0;i<ido;i++){ 480c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t1=i; 481c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2=i; 482c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) for(k=0;k<l1;k++){ 483c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) cc[t2]=ch[t1]; 484c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t1+=ido; 485c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) t2+=t10; 486c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 487c2e0dbddbe15c98d52c4786dac06cb8952a8ae6dTorne (Richard Coles) } 48868043e1e95eeb07d5cae7aca370b26518b0867d6Torne (Richard Coles) 489c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch L135: 4902a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t1=0; 4912a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t2=ido<<1; 4922a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t3=0; 4934e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t4=ipp2*t0; 4942a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for(j=1;j<ipph;j++){ 4952a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 4962a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t1+=t2; 497c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t3+=t0; 4982a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t4-=t0; 4992a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 5002a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t5=t1; 5012a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t6=t3; 502c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t7=t4; 5032a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 5042a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for(k=0;k<l1;k++){ 5052a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) cc[t5-1]=ch[t6]; 506c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch cc[t5]=ch[t7]; 5072a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t5+=t10; 5082a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t6+=ido; 509868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) t7+=ido; 510eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch } 511eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch } 512eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch 513868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) if(ido==1)return; 5142a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) if(nbd<l1)goto L141; 5152a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 5162a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t1=-ido; 5177dbb3d5cf0c15f500944d211057644d6a2f37371Ben Murdoch t3=0; 5187dbb3d5cf0c15f500944d211057644d6a2f37371Ben Murdoch t4=0; 5197dbb3d5cf0c15f500944d211057644d6a2f37371Ben Murdoch t5=ipp2*t0; 5202a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) for(j=1;j<ipph;j++){ 5212a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t1+=t2; 5222a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t3+=t2; 5232a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t4+=t0; 5242a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) t5-=t0; 5254e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t6=t1; 5264e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t7=t3; 5274e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t8=t4; 5284e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t9=t5; 5294e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) for(k=0;k<l1;k++){ 5304e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) for(i=2;i<ido;i+=2){ 5314e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) ic=idp2-i; 5324e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1]; 5334e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1]; 534a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) cc[i+t7]=ch[i+t8]+ch[i+t9]; 5354e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) cc[ic+t6]=ch[i+t9]-ch[i+t8]; 5364e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) } 537a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) t6+=t10; 5384e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t7+=t10; 5394e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t8+=ido; 5404e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t9+=ido; 5414e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) } 5424e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) } 5434e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) return; 5444e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) 5454e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) L141: 546c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch 5474e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t1=-ido; 5484e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t3=0; 5494e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t4=0; 5504e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t5=ipp2*t0; 5514e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) for(j=1;j<ipph;j++){ 5524e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t1+=t2; 5534e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t3+=t2; 5544e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t4+=t0; 5554e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t5-=t0; 5564e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) for(i=2;i<ido;i+=2){ 5574e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t6=idp2+t1-i; 558c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t7=i+t3; 5594e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t8=i+t4; 5604e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t9=i+t5; 5614e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) for(k=0;k<l1;k++){ 5624e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) cc[t7-1]=ch[t8-1]+ch[t9-1]; 5634e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) cc[t6-1]=ch[t8-1]-ch[t9-1]; 5644e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) cc[t7]=ch[t8]+ch[t9]; 5654e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) cc[t6]=ch[t9]-ch[t8]; 5664e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t6+=t10; 5674e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) t7+=t10; 568c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t8+=ido; 569c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch t9+=ido; 5704e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) } 5714e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) } 572c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch } 573c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch} 5744e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) 5754e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles)static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){ 5764e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) int i,k1,l1,l2; 5774e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) int na,kh,nf; 578c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch int ip,iw,ido,idl1,ix2,ix3; 579c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch 5804e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) nf=ifac[1]; 5814e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) na=1; 5824e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) l2=n; 583c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch iw=n; 584c5cede9ae108bb15f6b7a8aea21c7e1fefa2834cBen Murdoch 5854e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) for(k1=0;k1<nf;k1++){ 5864e180b6a0b4720a9b8e9e959a882386f690f08ffTorne (Richard Coles) kh=nf-k1; 587868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) ip=ifac[kh+1]; 588868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) l1=l2/ip; 5895d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) ido=n/l2; 5905d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) idl1=ido*l1; 5915d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) iw-=(ip-1)*ido; 5925d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) na=1-na; 5935d1f7b1de12d16ceb2c938c56701a3e8bfa558f7Torne (Richard Coles) 594868fa2fe829687343ffae624259930155e16dbd8Torne (Richard Coles) if(ip!=4)goto L102; 5952a99a7e74a7f215066514fe81d2bfa6639d9edddTorne (Richard Coles) 596 ix2=iw+ido; 597 ix3=ix2+ido; 598 if(na!=0) 599 dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1); 600 else 601 dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1); 602 goto L110; 603 604 L102: 605 if(ip!=2)goto L104; 606 if(na!=0)goto L103; 607 608 dradf2(ido,l1,c,ch,wa+iw-1); 609 goto L110; 610 611 L103: 612 dradf2(ido,l1,ch,c,wa+iw-1); 613 goto L110; 614 615 L104: 616 if(ido==1)na=1-na; 617 if(na!=0)goto L109; 618 619 dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1); 620 na=1; 621 goto L110; 622 623 L109: 624 dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1); 625 na=0; 626 627 L110: 628 l2=l1; 629 } 630 631 if(na==1)return; 632 633 for(i=0;i<n;i++)c[i]=ch[i]; 634} 635 636static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){ 637 int i,k,t0,t1,t2,t3,t4,t5,t6; 638 float ti2,tr2; 639 640 t0=l1*ido; 641 642 t1=0; 643 t2=0; 644 t3=(ido<<1)-1; 645 for(k=0;k<l1;k++){ 646 ch[t1]=cc[t2]+cc[t3+t2]; 647 ch[t1+t0]=cc[t2]-cc[t3+t2]; 648 t2=(t1+=ido)<<1; 649 } 650 651 if(ido<2)return; 652 if(ido==2)goto L105; 653 654 t1=0; 655 t2=0; 656 for(k=0;k<l1;k++){ 657 t3=t1; 658 t5=(t4=t2)+(ido<<1); 659 t6=t0+t1; 660 for(i=2;i<ido;i+=2){ 661 t3+=2; 662 t4+=2; 663 t5-=2; 664 t6+=2; 665 ch[t3-1]=cc[t4-1]+cc[t5-1]; 666 tr2=cc[t4-1]-cc[t5-1]; 667 ch[t3]=cc[t4]-cc[t5]; 668 ti2=cc[t4]+cc[t5]; 669 ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2; 670 ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2; 671 } 672 t2=(t1+=ido)<<1; 673 } 674 675 if(ido%2==1)return; 676 677L105: 678 t1=ido-1; 679 t2=ido-1; 680 for(k=0;k<l1;k++){ 681 ch[t1]=cc[t2]+cc[t2]; 682 ch[t1+t0]=-(cc[t2+1]+cc[t2+1]); 683 t1+=ido; 684 t2+=ido<<1; 685 } 686} 687 688static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1, 689 float *wa2){ 690 static float taur = -.5f; 691 static float taui = .8660254037844386f; 692 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10; 693 float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2; 694 t0=l1*ido; 695 696 t1=0; 697 t2=t0<<1; 698 t3=ido<<1; 699 t4=ido+(ido<<1); 700 t5=0; 701 for(k=0;k<l1;k++){ 702 tr2=cc[t3-1]+cc[t3-1]; 703 cr2=cc[t5]+(taur*tr2); 704 ch[t1]=cc[t5]+tr2; 705 ci3=taui*(cc[t3]+cc[t3]); 706 ch[t1+t0]=cr2-ci3; 707 ch[t1+t2]=cr2+ci3; 708 t1+=ido; 709 t3+=t4; 710 t5+=t4; 711 } 712 713 if(ido==1)return; 714 715 t1=0; 716 t3=ido<<1; 717 for(k=0;k<l1;k++){ 718 t7=t1+(t1<<1); 719 t6=(t5=t7+t3); 720 t8=t1; 721 t10=(t9=t1+t0)+t0; 722 723 for(i=2;i<ido;i+=2){ 724 t5+=2; 725 t6-=2; 726 t7+=2; 727 t8+=2; 728 t9+=2; 729 t10+=2; 730 tr2=cc[t5-1]+cc[t6-1]; 731 cr2=cc[t7-1]+(taur*tr2); 732 ch[t8-1]=cc[t7-1]+tr2; 733 ti2=cc[t5]-cc[t6]; 734 ci2=cc[t7]+(taur*ti2); 735 ch[t8]=cc[t7]+ti2; 736 cr3=taui*(cc[t5-1]-cc[t6-1]); 737 ci3=taui*(cc[t5]+cc[t6]); 738 dr2=cr2-ci3; 739 dr3=cr2+ci3; 740 di2=ci2+cr3; 741 di3=ci2-cr3; 742 ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2; 743 ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2; 744 ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3; 745 ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3; 746 } 747 t1+=ido; 748 } 749} 750 751static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1, 752 float *wa2,float *wa3){ 753 static float sqrt2=1.414213562373095f; 754 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8; 755 float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4; 756 t0=l1*ido; 757 758 t1=0; 759 t2=ido<<2; 760 t3=0; 761 t6=ido<<1; 762 for(k=0;k<l1;k++){ 763 t4=t3+t6; 764 t5=t1; 765 tr3=cc[t4-1]+cc[t4-1]; 766 tr4=cc[t4]+cc[t4]; 767 tr1=cc[t3]-cc[(t4+=t6)-1]; 768 tr2=cc[t3]+cc[t4-1]; 769 ch[t5]=tr2+tr3; 770 ch[t5+=t0]=tr1-tr4; 771 ch[t5+=t0]=tr2-tr3; 772 ch[t5+=t0]=tr1+tr4; 773 t1+=ido; 774 t3+=t2; 775 } 776 777 if(ido<2)return; 778 if(ido==2)goto L105; 779 780 t1=0; 781 for(k=0;k<l1;k++){ 782 t5=(t4=(t3=(t2=t1<<2)+t6))+t6; 783 t7=t1; 784 for(i=2;i<ido;i+=2){ 785 t2+=2; 786 t3+=2; 787 t4-=2; 788 t5-=2; 789 t7+=2; 790 ti1=cc[t2]+cc[t5]; 791 ti2=cc[t2]-cc[t5]; 792 ti3=cc[t3]-cc[t4]; 793 tr4=cc[t3]+cc[t4]; 794 tr1=cc[t2-1]-cc[t5-1]; 795 tr2=cc[t2-1]+cc[t5-1]; 796 ti4=cc[t3-1]-cc[t4-1]; 797 tr3=cc[t3-1]+cc[t4-1]; 798 ch[t7-1]=tr2+tr3; 799 cr3=tr2-tr3; 800 ch[t7]=ti2+ti3; 801 ci3=ti2-ti3; 802 cr2=tr1-tr4; 803 cr4=tr1+tr4; 804 ci2=ti1+ti4; 805 ci4=ti1-ti4; 806 807 ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2; 808 ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2; 809 ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3; 810 ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3; 811 ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4; 812 ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4; 813 } 814 t1+=ido; 815 } 816 817 if(ido%2 == 1)return; 818 819 L105: 820 821 t1=ido; 822 t2=ido<<2; 823 t3=ido-1; 824 t4=ido+(ido<<1); 825 for(k=0;k<l1;k++){ 826 t5=t3; 827 ti1=cc[t1]+cc[t4]; 828 ti2=cc[t4]-cc[t1]; 829 tr1=cc[t1-1]-cc[t4-1]; 830 tr2=cc[t1-1]+cc[t4-1]; 831 ch[t5]=tr2+tr2; 832 ch[t5+=t0]=sqrt2*(tr1-ti1); 833 ch[t5+=t0]=ti2+ti2; 834 ch[t5+=t0]=-sqrt2*(tr1+ti1); 835 836 t3+=ido; 837 t1+=t2; 838 t4+=t2; 839 } 840} 841 842static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1, 843 float *c2,float *ch,float *ch2,float *wa){ 844 static float tpi=6.283185307179586f; 845 int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10, 846 t11,t12; 847 float dc2,ai1,ai2,ar1,ar2,ds2; 848 int nbd; 849 float dcp,arg,dsp,ar1h,ar2h; 850 int ipp2; 851 852 t10=ip*ido; 853 t0=l1*ido; 854 arg=tpi/(float)ip; 855 dcp=cos(arg); 856 dsp=sin(arg); 857 nbd=(ido-1)>>1; 858 ipp2=ip; 859 ipph=(ip+1)>>1; 860 if(ido<l1)goto L103; 861 862 t1=0; 863 t2=0; 864 for(k=0;k<l1;k++){ 865 t3=t1; 866 t4=t2; 867 for(i=0;i<ido;i++){ 868 ch[t3]=cc[t4]; 869 t3++; 870 t4++; 871 } 872 t1+=ido; 873 t2+=t10; 874 } 875 goto L106; 876 877 L103: 878 t1=0; 879 for(i=0;i<ido;i++){ 880 t2=t1; 881 t3=t1; 882 for(k=0;k<l1;k++){ 883 ch[t2]=cc[t3]; 884 t2+=ido; 885 t3+=t10; 886 } 887 t1++; 888 } 889 890 L106: 891 t1=0; 892 t2=ipp2*t0; 893 t7=(t5=ido<<1); 894 for(j=1;j<ipph;j++){ 895 t1+=t0; 896 t2-=t0; 897 t3=t1; 898 t4=t2; 899 t6=t5; 900 for(k=0;k<l1;k++){ 901 ch[t3]=cc[t6-1]+cc[t6-1]; 902 ch[t4]=cc[t6]+cc[t6]; 903 t3+=ido; 904 t4+=ido; 905 t6+=t10; 906 } 907 t5+=t7; 908 } 909 910 if (ido == 1)goto L116; 911 if(nbd<l1)goto L112; 912 913 t1=0; 914 t2=ipp2*t0; 915 t7=0; 916 for(j=1;j<ipph;j++){ 917 t1+=t0; 918 t2-=t0; 919 t3=t1; 920 t4=t2; 921 922 t7+=(ido<<1); 923 t8=t7; 924 for(k=0;k<l1;k++){ 925 t5=t3; 926 t6=t4; 927 t9=t8; 928 t11=t8; 929 for(i=2;i<ido;i+=2){ 930 t5+=2; 931 t6+=2; 932 t9+=2; 933 t11-=2; 934 ch[t5-1]=cc[t9-1]+cc[t11-1]; 935 ch[t6-1]=cc[t9-1]-cc[t11-1]; 936 ch[t5]=cc[t9]-cc[t11]; 937 ch[t6]=cc[t9]+cc[t11]; 938 } 939 t3+=ido; 940 t4+=ido; 941 t8+=t10; 942 } 943 } 944 goto L116; 945 946 L112: 947 t1=0; 948 t2=ipp2*t0; 949 t7=0; 950 for(j=1;j<ipph;j++){ 951 t1+=t0; 952 t2-=t0; 953 t3=t1; 954 t4=t2; 955 t7+=(ido<<1); 956 t8=t7; 957 t9=t7; 958 for(i=2;i<ido;i+=2){ 959 t3+=2; 960 t4+=2; 961 t8+=2; 962 t9-=2; 963 t5=t3; 964 t6=t4; 965 t11=t8; 966 t12=t9; 967 for(k=0;k<l1;k++){ 968 ch[t5-1]=cc[t11-1]+cc[t12-1]; 969 ch[t6-1]=cc[t11-1]-cc[t12-1]; 970 ch[t5]=cc[t11]-cc[t12]; 971 ch[t6]=cc[t11]+cc[t12]; 972 t5+=ido; 973 t6+=ido; 974 t11+=t10; 975 t12+=t10; 976 } 977 } 978 } 979 980L116: 981 ar1=1.f; 982 ai1=0.f; 983 t1=0; 984 t9=(t2=ipp2*idl1); 985 t3=(ip-1)*idl1; 986 for(l=1;l<ipph;l++){ 987 t1+=idl1; 988 t2-=idl1; 989 990 ar1h=dcp*ar1-dsp*ai1; 991 ai1=dcp*ai1+dsp*ar1; 992 ar1=ar1h; 993 t4=t1; 994 t5=t2; 995 t6=0; 996 t7=idl1; 997 t8=t3; 998 for(ik=0;ik<idl1;ik++){ 999 c2[t4++]=ch2[t6++]+ar1*ch2[t7++]; 1000 c2[t5++]=ai1*ch2[t8++]; 1001 } 1002 dc2=ar1; 1003 ds2=ai1; 1004 ar2=ar1; 1005 ai2=ai1; 1006 1007 t6=idl1; 1008 t7=t9-idl1; 1009 for(j=2;j<ipph;j++){ 1010 t6+=idl1; 1011 t7-=idl1; 1012 ar2h=dc2*ar2-ds2*ai2; 1013 ai2=dc2*ai2+ds2*ar2; 1014 ar2=ar2h; 1015 t4=t1; 1016 t5=t2; 1017 t11=t6; 1018 t12=t7; 1019 for(ik=0;ik<idl1;ik++){ 1020 c2[t4++]+=ar2*ch2[t11++]; 1021 c2[t5++]+=ai2*ch2[t12++]; 1022 } 1023 } 1024 } 1025 1026 t1=0; 1027 for(j=1;j<ipph;j++){ 1028 t1+=idl1; 1029 t2=t1; 1030 for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++]; 1031 } 1032 1033 t1=0; 1034 t2=ipp2*t0; 1035 for(j=1;j<ipph;j++){ 1036 t1+=t0; 1037 t2-=t0; 1038 t3=t1; 1039 t4=t2; 1040 for(k=0;k<l1;k++){ 1041 ch[t3]=c1[t3]-c1[t4]; 1042 ch[t4]=c1[t3]+c1[t4]; 1043 t3+=ido; 1044 t4+=ido; 1045 } 1046 } 1047 1048 if(ido==1)goto L132; 1049 if(nbd<l1)goto L128; 1050 1051 t1=0; 1052 t2=ipp2*t0; 1053 for(j=1;j<ipph;j++){ 1054 t1+=t0; 1055 t2-=t0; 1056 t3=t1; 1057 t4=t2; 1058 for(k=0;k<l1;k++){ 1059 t5=t3; 1060 t6=t4; 1061 for(i=2;i<ido;i+=2){ 1062 t5+=2; 1063 t6+=2; 1064 ch[t5-1]=c1[t5-1]-c1[t6]; 1065 ch[t6-1]=c1[t5-1]+c1[t6]; 1066 ch[t5]=c1[t5]+c1[t6-1]; 1067 ch[t6]=c1[t5]-c1[t6-1]; 1068 } 1069 t3+=ido; 1070 t4+=ido; 1071 } 1072 } 1073 goto L132; 1074 1075 L128: 1076 t1=0; 1077 t2=ipp2*t0; 1078 for(j=1;j<ipph;j++){ 1079 t1+=t0; 1080 t2-=t0; 1081 t3=t1; 1082 t4=t2; 1083 for(i=2;i<ido;i+=2){ 1084 t3+=2; 1085 t4+=2; 1086 t5=t3; 1087 t6=t4; 1088 for(k=0;k<l1;k++){ 1089 ch[t5-1]=c1[t5-1]-c1[t6]; 1090 ch[t6-1]=c1[t5-1]+c1[t6]; 1091 ch[t5]=c1[t5]+c1[t6-1]; 1092 ch[t6]=c1[t5]-c1[t6-1]; 1093 t5+=ido; 1094 t6+=ido; 1095 } 1096 } 1097 } 1098 1099L132: 1100 if(ido==1)return; 1101 1102 for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik]; 1103 1104 t1=0; 1105 for(j=1;j<ip;j++){ 1106 t2=(t1+=t0); 1107 for(k=0;k<l1;k++){ 1108 c1[t2]=ch[t2]; 1109 t2+=ido; 1110 } 1111 } 1112 1113 if(nbd>l1)goto L139; 1114 1115 is= -ido-1; 1116 t1=0; 1117 for(j=1;j<ip;j++){ 1118 is+=ido; 1119 t1+=t0; 1120 idij=is; 1121 t2=t1; 1122 for(i=2;i<ido;i+=2){ 1123 t2+=2; 1124 idij+=2; 1125 t3=t2; 1126 for(k=0;k<l1;k++){ 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]; 1129 t3+=ido; 1130 } 1131 } 1132 } 1133 return; 1134 1135 L139: 1136 is= -ido-1; 1137 t1=0; 1138 for(j=1;j<ip;j++){ 1139 is+=ido; 1140 t1+=t0; 1141 t2=t1; 1142 for(k=0;k<l1;k++){ 1143 idij=is; 1144 t3=t2; 1145 for(i=2;i<ido;i+=2){ 1146 idij+=2; 1147 t3+=2; 1148 c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3]; 1149 c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1]; 1150 } 1151 t2+=ido; 1152 } 1153 } 1154} 1155 1156static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){ 1157 int i,k1,l1,l2; 1158 int na; 1159 int nf,ip,iw,ix2,ix3,ido,idl1; 1160 1161 nf=ifac[1]; 1162 na=0; 1163 l1=1; 1164 iw=1; 1165 1166 for(k1=0;k1<nf;k1++){ 1167 ip=ifac[k1 + 2]; 1168 l2=ip*l1; 1169 ido=n/l2; 1170 idl1=ido*l1; 1171 if(ip!=4)goto L103; 1172 ix2=iw+ido; 1173 ix3=ix2+ido; 1174 1175 if(na!=0) 1176 dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1); 1177 else 1178 dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1); 1179 na=1-na; 1180 goto L115; 1181 1182 L103: 1183 if(ip!=2)goto L106; 1184 1185 if(na!=0) 1186 dradb2(ido,l1,ch,c,wa+iw-1); 1187 else 1188 dradb2(ido,l1,c,ch,wa+iw-1); 1189 na=1-na; 1190 goto L115; 1191 1192 L106: 1193 if(ip!=3)goto L109; 1194 1195 ix2=iw+ido; 1196 if(na!=0) 1197 dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1); 1198 else 1199 dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1); 1200 na=1-na; 1201 goto L115; 1202 1203 L109: 1204/* The radix five case can be translated later..... */ 1205/* if(ip!=5)goto L112; 1206 1207 ix2=iw+ido; 1208 ix3=ix2+ido; 1209 ix4=ix3+ido; 1210 if(na!=0) 1211 dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); 1212 else 1213 dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); 1214 na=1-na; 1215 goto L115; 1216 1217 L112:*/ 1218 if(na!=0) 1219 dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1); 1220 else 1221 dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1); 1222 if(ido==1)na=1-na; 1223 1224 L115: 1225 l1=l2; 1226 iw+=(ip-1)*ido; 1227 } 1228 1229 if(na==0)return; 1230 1231 for(i=0;i<n;i++)c[i]=ch[i]; 1232} 1233 1234void spx_drft_forward(struct drft_lookup *l,float *data){ 1235 if(l->n==1)return; 1236 drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache); 1237} 1238 1239void spx_drft_backward(struct drft_lookup *l,float *data){ 1240 if (l->n==1)return; 1241 drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache); 1242} 1243 1244void spx_drft_init(struct drft_lookup *l,int n) 1245{ 1246 l->n=n; 1247 l->trigcache=(float*)speex_alloc(3*n*sizeof(*l->trigcache)); 1248 l->splitcache=(int*)speex_alloc(32*sizeof(*l->splitcache)); 1249 fdrffti(n, l->trigcache, l->splitcache); 1250} 1251 1252void spx_drft_clear(struct drft_lookup *l) 1253{ 1254 if(l) 1255 { 1256 if(l->trigcache) 1257 speex_free(l->trigcache); 1258 if(l->splitcache) 1259 speex_free(l->splitcache); 1260 } 1261} 1262