c--lpfiltd.f c------------------------------------------------------------------------------ c Exchange alogorithm for designing low pass finite length c response digital filters with linear phase. c c*** This is a modified version *** D. Crossley, October 1989 *** c modifications : c all routines made double precision c substituting DFOURT for Singleton's FFT c c J.H. McClellan and T.W. Parks---Rice University---December 11, 1971 c The filter length (NF0, the passband cutoff frequency (FP), c the stopband freqnecy (FS) and the weighting factor (AA) are c inputs to the alogrithm. Then the stopband deviation (DS) is c minimised and thus the passband deviation (DP=AA*DS) is also c minimised. The program seeks to determine the best location of c the peaks of the error curve in order to mimimise the deviations. c c Reference--Parks, T.W. and J.H. McClellan, "Chebychev approximation c of non-recursive digital filters with linear phase", IEEE Trans. c on Circuit Theory, Vol CT-19, March 1972. c c The program appears in--Parks, T.W. and J.H. McClellan, "A program c for the design of linear phase finite impulse response digital filters" c IEEE Trans. on Audio and Electroacoustics, Vol. AU-20, No. 3, August 1972. c c The program is set up in a unending loop so that any number of cases c can be run at one time. The loop is terminated by entering 0 for c the filter length (NF). c c The input data format for an individual case is as follows : c there are five cards necessary for each case. c Card one--the filter length is entered in the first 3 columns c in an integer format. The filter length must be an odd integer, c if not the program will make it an odd integer. c Card two-- the passband is entered in the first ten columns in a real c format. The passband frequency is normalised so the Nyquist c frequency is equal to one. Thus the maximum frequency is 0.5. c Card three--stopband cutoff frequency same as card two. c Card four--Weighting factor is entered in a real format in the first c ten columns. There are no resrictions on the value of the c weighting factor except that it be non-zero. c Card five--Impulse response option flag. If the impulse response c is desired enter a 1 in col. 1 otherwise leave it blank. c A sample case setup c 021 c 0.2 c 0.3 c 10.0 c 1 c The foregoing case is for NF=21, FP=0.2, FS=0.3, weighting factor=10.0 c and the impulse response will be obtained. c c......................................................................... implicit real*8(a-h,o-z) common/calc/x(1025),scale,af(8300),ad(1025),y(1025) common/ppf/rho,delf,nf common/cheb/n,ngrid,kgrid,kfp,kfs,ifr(1025) common/chek/nl,ng,ka,kb dimension fopt(1025),b(2049) character*60 title character*20 dunit,tunit,fmt character*1 ans data dunit,tunit,fmt/'amplitude','sec','(4e15.7)'/, & delt/5.0/ c c initialisation of the grid of 16(n-1) points, the constant lgrid c can be changed if a different grid density is desired c c--initialise the files with fixed names open(1,file='lpfiltd.dat') open(2,file='lpfiltd.out') open(3,file='lpfiltd.in') c--read in time unit read(3,'(a20)') tunit c--read in original sampling interval read(3,*) delt1 write(*,'(a,f8.5,2x,a)') 'original sampling ',delt1,tunit c lgrid=16 c lgrid=8 68 continue c-- c write(*,300) 300 format(' Enter filter length (odd integer <= 2049) --> '$) read (3,*) nf if(nf.lt.3.or.nf.gt.2049) then write(*,*) 'filter length too small or too big' stop endif c-- c write(*,301) c 301 format(' Enter passband cutoff (FP) --> '$) c read (*,302) fp c 302 format(f20.3) c c--read new sampling interval in original time units read(3,*) delt2 fs = delt1/(2*delt2) write(*,'(a,f10.7)') 'stopband frequency ',fs c c--read passband cutoff period in original time units read(3,*) fp fp = delt1/fp write(*,'(a,f10.7)') 'passband frequency ',fp c-- c write(*,303) c 303 format(' Enter stopband cutoff (FS) --> '$) c read (*,302) fs c write(*,304) c 304 format(' Enter weighting factor (AA) --> '$) c read (*,302) aa c c--read weighting factor c--(size of ripples in passband / size of ripples in stopband) read(3,*) aa c-- n=nf/2+1 nf=2*n-1 pi2=8.d0*datan(1.d0) kcount=0 tw=fs-fp ngrid=lgrid*(n-1) kgrid=2 delf=0.5/ngrid af(1)=0.d0 do 2 j=1,ngrid 2 af(j+1)=j*delf kfp=fp/delf-1 kfs=fs/delf c-- 3 continue if(fp.lt.af(kfp+1)) goto 4 kfp=kfp+1 goto 3 c-- 4 continue if(fs.le.af(kfs)) goto 5 kfs=kfs+1 goto 4 c-- 5 continue af(kfp)=fp af(kfs)=fs ngrid=ngrid+1 do 11 j=1,ngrid 11 af(j)=dcos(pi2*af(j)) c c initial guess for the n+1 optimal frequencies c kp=(kfp-1)/2 kp=2*kp+1 ks=kfs/2 ks=2*ks+1 xpp=kp-1 xss=ngrid-ks xaa=xpp+xss pt=xaa/(n-1) if(pt.gt.3) goto 101 write(2,106) 106 format(' Error in input data') stop c-- 101 kr=xpp/pt if(kr.eq.0) kr=kr+1 pt=xpp/kr ifr(1)=1 do 6 j=1,kr xtw=j*pt ntw=xtw+1 ntw=ntw/2 6 ifr(j+1)=2*ntw+1 kr=kr+1 ifr(kr)=kp ifr(kr+1)=ks ntw=n-kr if(ntw.ge.0) goto 102 c-- write(2,106) stop c-- 102 if(ntw.eq.0) goto 8 pt=xss/ntw do 7 j=1,ntw xtw=j*pt+ks-1.d0 nb=xtw nb=nb/2 7 ifr(kr+1+j)=2*nb-1 8 continue c c Main calculation routine c calculate rho on the optimal set c The first 4 passes are made on a grid of 8(n-1) points then the c finer grid of 16(n-1) points is used c A maximum of 20 iterations is allowed c 100 continue do 12 j=1,n+1 12 x(j)=af(ifr(j)) m=(n+1)/2 do 13 j=1,n+1 13 ad(j)=d(j,n+1,m) r1=0.d0 do 14 j=1,kr 14 r1=r1+ad(j) r2=0.d0 k=1 do 15 j=1,kr r2=r2+ad(j)*aa*k 15 k=-k do 16 j=kr+1,n+1 r2=r2+k*ad(j) 16 k=-k rho=r1/r2 write(2,23) rho nu=1 if(rho.gt.0.d0) nu=-1 rho=dabs(rho) k=nu xray=x(kr+1) do 17 j=1,kr ad(j)=ad(j)*(xray-x(j)) y(j)=1.d0+aa*k*rho 17 k=-k k=-k do 18 j=kr+1,n x(j)=x(j+1) ad(j)=ad(j+1)*(xray-x(j)) y(j)=k*rho 18 k=-k y(n+1)=k*rho nl=0 ng=0 test=gee(ifr(kr+1)) ro=dabs(test-rho)/rho if(ro.lt.0.1.or.kcount.lt.4) goto 801 c-- write(2,802) 802 format(' Interpolation error') stop c-- 801 continue c c search for new extremal frequencies c ntot=0 npt=ifr(1) ng=ifr(2) ka=1 kb=1 ksign=nu xt=y(1) nfir=ifr(1) nold=ifr(n+1) if(kr.eq.1) goto 400 if(npt-1)202,200,202 202 zt=gee(npt-kgrid) if(ksign*(zt-xt))200,200,201 201 kpt=npt-kgrid 220 kpt=kpt-kgrid xt=zt if(kpt-1)210,211,211 211 zt=gee(kpt) if(ksign*(zt-xt))210,210,220 210 ifr(1)=kpt+kgrid ntot=ntot+1 ksign=-ksign goto 260 200 zt=gee(npt+kgrid) if(ksign*(zt-xt))230,230,231 231 kpt=npt 240 kpt=kpt+kgrid xt=zt zt=gee(kpt+kgrid) if(ksign*(zt-xt))250,250,240 250 ifr(1)=kpt ntot=ntot+1 230 ksign=-ksign 260 xfir=dabs(1.d0-xt)/aa ka=2 440 if(ka.gt.(n+1)) goto 441 nl=npt ng=ngrid+1 if(ka.le.n) ng=ifr(ka+1) npt=ifr(ka) xt=y(ka)*kb if(ka.eq.kr) goto 400 if(ka.eq.(kr+1)) goto 401 goto 403 400 if((kfp-npt).lt.kgrid) goto 407 npt=kfp ntot=1 goto 407 401 kb=-1 if((npt-kfs).lt.kgrid) goto 407 npt=kfs ntot=1 goto 407 403 kpt=npt zt=gee(kpt-kgrid) if(ksign*(zt-xt))404,404,405 405 xt=zt kpt=kpt-kgrid zt=gee(kpt-kgrid) if(ksign*(zt-xt))406,406,405 406 ifr(ka)=kpt ksign=-ksign ntot=ntot+1 ka=ka+1 goto 440 404 if((kpt+kgrid).gt.ngrid) goto 407 zt=gee(kpt+kgrid) if(ksign*(zt-xt))407,407,408 408 xt=zt kpt=kpt+kgrid if((kgrid+kpt).gt.ngrid) goto 409 zt=gee(kpt+kgrid) if(ksign*(zt-xt))409,409,408 409 ifr(ka)=kpt ksign=-ksign ntot=ntot+1 ka=ka+1 goto 440 407 ifr(ka)=npt ksign=-ksign ka=ka+1 goto 440 c-- 441 continue xold=dabs(xt)*aa nl=0 ng=0 if(ifr(1).eq.1) goto 450 if(ifr(n+1).eq.ngrid) goto 451 write(2,452) 452 format(' Error-neither endpoint is an extremum') stop c-- 451 if(nfir.eq.1) goto 444 test=gee(1) test=-nu*(test-1.d0)*xold if(test.lt.0.d0) goto 444 ntot=1 do 445 j=1,n 445 ifr(n+2-j)=ifr(n+1-j) ifr(1)=1 kr=kr+1 goto 444 450 if(ifr(n+1).eq.ngrid) goto 444 if(nold.eq.ngrid) goto 444 test=gee(ngrid) test=ksign*test-xfir if(test.lt.0.d0) goto 444 ntot=1 do 446 j=1,n 446 ifr(j)=ifr(j+1) ifr(n+1)=ngrid kr=kr-1 c-- 444 continue if(kr.gt.0) goto 454 write(2,453) 453 format(' Error') stop c-- 454 continue write(2,23) rho 23 format(' rho ',d20.12) kcount=kcount+1 if(kcount.gt.50) goto 448 if(kcount.eq.4.and.kgrid.eq.2) goto 447 if(ntot.ne.0) goto 100 if(kgrid.eq.2) goto 447 goto 448 447 kgrid=1 goto 100 c-- 448 continue c c output of the filter design c if the impulse is desired use the following 5 statements c to call the subroutine pref which in turn calls the FFT c jkl=0 c write(*,83) c 83 format(' For the filter response enter a 1 --> '$) c read 81,jkl c 81 format(i1) read(3,*) jkl if(jkl.eq.1) call pref(b) c-- c write(*,'(a$)') ' title --> ' c read(*,'(a)') title read(3,'(a60)') title write(1,'(a)') title c write(*,'(a$)') ' delt --> ' c read(*,*) delt write(1,'(i10,g12.5)') nf,delt1 write(1,'(3a20)') fmt,dunit,tunit write(1,fmt) (b(j),j=1,nf) c-- do 72 j=1,n+1 72 fopt(j)=delf*(ifr(j)-1) fopt(kr)=fp fopt(kr+1)=fs c-- write(2,21) (fopt(j),j=1,n+1) 21 format(' Extremal frequencies'/(f10.7)) krp=n-kr r1=aa*rho c-- write(2,22) nf,kr,krp,tw,aa,rho,r1,kcount 22 format(' Filter length = ',i5/' Passband ripples = ',i5/ 1' Stopband ripples = ',i5/' Transition width = ',f10.7/ 2' Weighting factor = ',f12.7/' Deviation in stopband = ',f14.11/ 3' Deviation in passband = ',f14.11/' Number of iterations = ',i5) write(2,95) fp,fs 95 format(' Passband edge = ',f10.7/' Stopband edge = ',f10.7) write(*,'(a$)') ' another filter ? [y/n] --> ' read(*,'(a1)') ans if(ans.eq.'y') goto 68 stop c-- end c---------------------------------------------------------------------------- function d(k,n,m) c c The subroutine d(,,,) calculates the Lagrangian interpolation c coefficients for use in the barycentric interpolation formula c c......................................................................... implicit real*8(a-h,o-z) common/calc/x(1025),scale,af(8300),ad(1025),y(1025) c-- d=1.d0 q=x(k) do 1 j=1,m if(j-k)2,3,2 2 d=2.d0*d*(x(j)-q) 3 if(n+1-j-k)4,1,4 4 d=2.d0*d*(x(n+1-j)-q) 1 continue if(n-2*m)5,6,5 5 if(m+1-k)7,6,7 7 d=2.d0*d*(x(m+1)-q) 6 d=1.d0/d c-- return end c--------------------------------------------------------------------------- subroutine pref(b) c c Any FFT can be used, this uses Singleton's FFT c ***modified 23 October 1989 to use DFOURT (Brennan) *** c The subroutine pref samples the frequency response at 2**m points c so that the FFT (R.C.Singleton) can be called to inverse transform c the data to obtain the impulse response c c......................................................................... implicit real*8(a-h,o-z) complex*16 cdata(1025),work(1025) common/calc/x(1025),scale,af(8300),ad(1025),y(1025) common/ppf/rho,delf,nf common/cheb/n,ngrid,kgrid,kfp,kfs,ifr(1025) dimension a(2049),b(2049),rf(2049) equivalence (a,cdata) c-- fsh=1.0d-08 pi2=8.d0*datan(1.d0) nx=2 1 if(nx.ge.nf) goto 2 nx=2*nx goto 1 2 nn=nx nx=nx/2 xnn=nn df=1.d0/xnn if(ifr(1).eq.1) goto 3 p=goo(1.d0) goto 4 3 p=y(1) 4 a(1)=p if(ifr(n+1).eq.ngrid) goto 5 p=goo(-1.d0) goto 6 5 p=y(n) 6 a(nx+1)=p l=1 do 7 j=1,nx-1 at=df*j at=dcos(pi2*at) 10 as=x(l) if(at.gt.as) goto 8 if((as-at).lt.fsh) goto 9 l=l+1 goto 10 9 a(j+1)=y(l) goto 11 8 if((at-as).lt.fsh) goto 9 a(j+1)=goo(at) 11 a(nn+1-j)=a(j+1) b(j+1)=0.d0 if(l.gt.1) l=l-1 7 b(nn+1-j)=0.d0 b(1)=0.d0 b(nx+1)=0.d0 c call fft(a,b,nn,nn,nn,1) c do 12 j=1,nx c 12 b(j)=a(nx+1-j)/nn c-- nyq=nn/2+1 write(2,'(1x,3i5)') nf,nyq,n write(2,100) (a(i),i=1,nyq) 100 format(1x,f16.8) call dfourt(a,nn,1,-1,0,work,1025) c c--find coefficients c do 110 i=1,n rf(i)=real(cdata(i))/nn 110 continue c c--transfer to two-sided array b and find sum c do 12 j=1,n b(j) = rf(n-j+1) if(j.eq.n) goto 12 b(j+n) = rf(j+1) 12 continue c c--find sum of coefficients c sum=0. do 115 i=1,nf sum=sum+b(i) 115 continue write(2,120) sum 120 format(' sum of coefficients = ',g12.5) c c--write out impulse response c write(2,130) (i,b(i),i=1,nf) 130 format(/' impulse response'/(i10,f16.8)) c return end c------------------------------------------------------------------- function gee(k) c c The subroutine gee evaluates the freqnecy response using the c interpolation coefficients which are calculated in the routine d(,,,) c c......................................................................... implicit real*8(a-h,o-z) common/calc/x(1025),scale,af(8300),ad(1025),y(1025) common/ppf/rho,delf,nf common/cheb/n,ngrid,kgrid,kfp,kfs,ifr(1025) common/chek/nl,ng,ka,kb c-- if(k-nl)3,2,3 3 if(k-ng)4,2,4 4 p=0.d0 xf=af(k) d=0.d0 do 1 j=1,n c=ad(j)/(xf-x(j)) d=d+c 1 p=p+c*y(j) gee=p/d c-- return 2 gee=y(ka+kb) c-- return end c------------------------------------------------------------------- function goo(f) c c The subroutine goo is the same as gee except that it is called by c a frequency value rather than by a grid index value c c......................................................................... implicit real*8(a-h,o-z) common/calc/x(1025),scale,af(8300),ad(1025),y(1025) common/ppf/rho,delf,nf common/cheb/n,ngrid,kgrid,kfp,kfs,ifr(1025) c-- xf=f p=0.d0 d=0.d0 do 1 j=1,n c=ad(j)/(xf-x(j)) d=d+c 1 p=p+c*y(j) goo=p/d c-- return end c---------------------------------------------------------------------------- SUBROUTINE DFOURT(DATA,N,NDIM,ISIGN,IFORM,WORK,NWORK) c c--Brenner's version of FOURT in double precision c C COOLEY-TUKEY FAST FOURIER TRANSFORM IN USASI BASIC FORTRAN. C MULTI-DIMENSIONAL TRANSFORM, DIMENSIONS OF ARBITRARY SIZE, C COMPLEX OR REAL DATA. N POINTS CAN BE TRANSFORMED IN TIME C PROPORTIONAL TO N*LOG(N), WHEREAS OTHER METHODS TAKE N**2 TIME. C FURTHERMORE, LESS ERROR IS BUILT UP. WRITTEN BY NORMAN BRENNER C OF MIT LINCOLN LABORATORY, JUNE 1968. C C DIMENSION DATA(N(1),N(2),...),TRANSFORM(N(1),N(2),...),N(NDIM) C TRANSFORM(K1,K2,...) = SUM(DATA(J1,J2,...)*EXP(ISIGN*2*PI*SQRT(-1) C *((J1-1)*(K1-1)/N(1)+(J2-1)*(K2-1)/N(2)+...))), SUMMED FOR ALL C J1 AND K1 FROM 1 TO N(1), J2 AND K2 FROM 1 TO N(2), ETC. FOR ALL C NDIM SUBSCRIPTS. NDIM MUST BE POSITIVE AND EACH N(IDIM) MAY BE C ANY INTEGER. ISIGN IS +1 OR -1. LET NTOT = N(1)*N(2)... C ...*N(NDIM). THEN A -1 TRANSFORM FOLLOWED BY A +1 ONE C (OR VICE VERSA) RETURNS NTOT TIMES THE ORIGINAL DATA. C IFORM = 1, 0 OR -1, AS DATA IS COMPLEX, REAL OR THE C FIRST HALF OF A COMPLEX ARRAY. TRANSFORM VALUES ARE C RETURNED TO ARRAY DATA. THEY ARE COMPLEX, REAL OR C THE FIRST HALF OF A COMPLEX ARRAY, AS IFORM = 1, -1 OR 0. C THE TRANSFORM OF A REAL ARRAY (IFORM = 0) DIMENSIONED N(1) BY N(2) C BY ... WILL BE RETURNED IN THE SAME ARRAY, NOW CONSIDERED TO C BE COMPLEX OF DIMENSIONS N(1)/2+1 BY N(2) BY .... NOTE THAT IF C IFORM = 0 OR -1, N(1) MUST BE EVEN, AND ENOUGH ROOM MUST BE C RESERVED. THE MISSING VALUES MAY BE OBTAINED BY COMPLEX CONJU- C GATION. THE REVERSE TRANSFORMATION, OF A HALF COMPLEX ARRAY C DIMENSIONED N(1)/2+1 BY N(2) BY ..., IS ACCOMPLISHED SETTING IFORM C TO -1. IN THE N ARRAY, N(1) MUST BE THE TRUE N(1), NOT N(1)/2+1. C THE TRANSFORM WILL BE REAL AND RETURNED TO THE INPUT ARRAY. C WORK IS A ONE-DIMENSIONAL COMPLEX ARRAY USED FOR WORKING STORAGE. C ITS LENGTH, NWORK, NEED NEVER BE LARGER THAN THE LARGEST N(IDIM) C AND FREQUENTLY MAY BE MUCH SMALLER. FOURT COMPUTES THE MINIMUM C LENGTH WORKING STORAGE REQUIRED AND CHECKS THAT NWORK IS AT LEAST C AS LONG. THIS MINIMUM LENGTH IS COMPUTED AS SHOWN BELOW. C C FOR EXAMPLE-- C DIMENSION DATA(1960),WORK(10) C COMPLEX DATA,WORK C CALL FOURT(DATA,1960,1,-1,+1,WORK,10) C C THE MULTI-DIMENSIONAL TRANSFORM IS BROKEN DOWN INTO ONE-DIMEN- C SIONAL TRANSFORMS OF LENGTH N(IDIM). THESE ARE FURTHER BROKEN C DOWN INTO TRANSFORMS OF LENGTH IFACT(IF), WHERE THESE ARE THE C PRIME FACTORS OF N(IDIM). FOR EXAMPLE, N(1) = 1960, IFACT(IF) = C 2, 2, 2, 5, 7 AND 7. THE RUNNING TIME IS PROPORTIONAL TO NTOT * C SUM(IFACT(IF)), THOUGH FACTORS OF TWO AND THREE WILL RUN ESPE- C CIALLY FAST. NAIVE TRANSFORM PROGRAMS WILL RUN IN TIME NTOT**2. C ARRAYS WHOSE SIZE NTOT IS PRIME WILL RUN MUCH SLOWER THAN THOSE C WITH COMPOSITE NTOT. FOR EXAMPLE, NTOT = N(1) = 1951 (A PRIME), C RUNNING TIME WILL BE 1951*1951, WHILE FOR NTOT = 1960, IT WILL C BE 1960*(2+2+2+5+7+7), A SPEEDUP OF EIGHTY TIMES. NAIVE CALCUL- C ATION WILL RUN BOTH IN THE SLOWER TIME. IF AN ARRAY IS OF C INCONVENIENT LENGTH, SIMPLY ADD ZEROES TO PAD IT OUT. THE RESULTS C WILL BE INTERPOLATED ACCORDING TO THE NEW LENGTH (SEE BELOW). C C A FOURIER TRANSFORM OF LENGTH IFACT(IF) REQUIRES A WORK ARRAY C OF THAT LENGTH. THEREFORE, NWORK MUST BE AS BIG AS THE LARGEST C PRIME FACTOR. FURTHER, WORK IS NEEDED FOR DIGIT REVERSAL-- C EACH N(IDIM) (BUT N(1)/2 IF IFORM = 0 OR -1) IS FACTORED SYMMETRI- C CALLY, AND NWORK MUST BE AS BIG AS THE CENTER FACTOR. (TO FACTOR C SYMMETRICALLY, SEPARATE PAIRS OF IDENTICAL FACTORS TO THE FLANKS, C COMBINING ALL LEFTOVERS IN THE CENTER.) FOR EXAMPLE, N(1) = 1960 C =2*2*2*5*7*7=2*7*10*7*2, SO NWORK MUST AT LEAST MAX(7,10) = 10. C C AN UPPER BOUND FOR THE RMS RELATIVE ERROR IS GIVEN BY GENTLEMAN C AND SANDE (3)-- 3 * 2**(-B) * SUM(F**1.5), WHERE 2**(-B) IS THE C SMALLEST BIT IN THE FLOATING POINT FRACTION AND THE SUM IS OVER C THE PRIME FACTORS OF NTOT. C C IF THE INPUT DATA ARE A TIME SERIES, WITH INDEX J REPRESENTING C A TIME (J-1)*DELTAT, THEN THE CORRESPONDING INDEX K IN THE C TRANSFORM REPRESENTS THE FREQUENCY (K-1)*2*PI/(N*DELTAT), WHICH C BY PERIODICITY, IS THE SAME AS FREQUENCY -(N-K+1)*2*PI/(N*DELTAT). C THIS IS TRUE FOR N = EACH N(IDIM) INDEPENDENTLY. C C REFERENCES-- C 1. COOLEY, J.W. AND TUKEY, J.W., AN ALGORITHM FOR THE MACHINE C CALCULATION OF COMPLEX FOURIER SERIES. MATH. COMP., 19, 90, C (APRIL 1967), 297-301. C 2. RADER, C., ET AL., WHAT IS THE FAST FOURIER TRANSFORM, IEEE C TRANSACTIONS ON AUDIO AND ELECTROACOUSTICS, AU-15, 2 (JUNE 1967). C (SPECIAL ISSUE ON THE FAST FOURIER TRANSFORM AND ITS APPLICATIONS) C 3. GENTLEMAN, W.M. AND SANDE, G., FAST FOURIER TRANSFORMS-- C FOR FUN AND PROFIT. 1966 FALL JOINT COMP. CONF., SPARTAN BOOKS, C WASHINGTON, 1966. C 4. GOERTZEL, G., AN ALGORITHM FOR THE EVALUATION OF FINITE C TRIGONOMETRIC SERIES. AM. MATH. MO., 65, (1958), 34-35. C 5. SINGLETON, R.C., A METHOD FOR COMPUTING THE FAST FOURIER C TRANSFORM WITH AUXILIARY MEMORY AND LIMITED HIGH-SPEED STORAGE. C IN (2). c c...................................................................... c implicit real*8(a-h,o-z) DIMENSION DATA(1), N(1), WORK(1), IFSYM(32), IFCNT(10), IFACT(32) c IF (IFORM) 10,10,40 10 IF (N(1)-2*(N(1)/2)) 20,40,20 20 WRITE (6,30) IFORM,(N(IDIM),IDIM=1,NDIM) 30 FORMAT (26H0ERROR IN FOURT. IFORM = ,I2,23H (REAL OR HALF-COMPLEX $)23H, BUT N(1) IS NOT EVEN./14H DIMENSIONS = 20I5) RETURN 40 NTOT=1 DO 50 IDIM=1,NDIM 50 NTOT=NTOT*N(IDIM) NREM=NTOT IF (IFORM) 60,70,70 60 NREM=1 NTOT=(NTOT/N(1))*(N(1)/2+1) C LOOP OVER ALL DIMENSIONS. 70 DO 230 JDIM=1,NDIM IF (IFORM) 80,90,90 80 IDIM=NDIM+1-JDIM GO TO 100 90 IDIM=JDIM NREM=NREM/N(IDIM) 100 NCURR=N(IDIM) IF (IDIM-1) 110,110,140 110 IF (IFORM) 120,130,140 120 CALL DFIXRL (DATA,N(1),NREM,ISIGN,IFORM) NTOT=(NTOT/(N(1)/2+1))*N(1) 130 NCURR=NCURR/2 140 IF (NCURR-1) 190,190,150 C FACTOR N(IDIM), THE LENGTH OF THIS DIMENSION. 150 CALL DFACTR (NCURR,IFACT,NFACT) IFMAX=IFACT(NFACT) C ARRANGE THE FACTORS SYMMETRICALLY FOR SIMPLER DIGIT REVERSAL. CALL DSMFAC (IFACT,NFACT,ISYM,IFSYM,NFSYM,ICENT,IFCNT,NFCNT) IFMAX=MAX0(IFMAX,ICENT) IF (IFMAX-NWORK) 180,180,160 160 WRITE (6,170) NWORK,IDIM,NCURR,ICENT,(IFACT(IF),IF=1,NFACT) 170 FORMAT (26H0ERROR IN FOURT. NWORK = ,I4,20H IS TOO SMALL FOR N(, $I1,4H) = ,I5,17H, WHOSE CENTER = ,I4,31H, AND WHOSE PRIME FACTORS $ARE--/(1X,20I5)) RETURN 180 NPREV=NTOT/(N(IDIM)*NREM) C DIGIT REVERSE ON SYMMETRIC FACTORS, FOR EXAMPLE 2*7*6*7*2. CALL DSYMRV (DATA,NPREV,NCURR,NREM,IFSYM,NFSYM) C DIGIT REVERSE THE ASYMMETRIC CENTER, FOR EXAMPLE, ON 6 = 2*3. CALL DASMRV (DATA,NPREV*ISYM,ICENT,ISYM*NREM,IFCNT,NFCNT,WORK) C FOURIER TRANSFORM ON EACH FACTOR, FOR EXAMPLE, ON 2,7,2,3,7 AND 2. CALL DCOOL (DATA,NPREV,NCURR,NREM,ISIGN,IFACT,WORK) 190 IF (IFORM) 200,210,230 200 NREM=NREM*N(IDIM) GO TO 230 210 IF (IDIM-1) 220,220,230 220 CALL DFIXRL (DATA,N(1),NREM,ISIGN,IFORM) NTOT=NTOT/N(1)*(N(1)/2+1) 230 CONTINUE RETURN END c------------------------------------------------------------------------- SUBROUTINE DFACTR (N,IFACT,NFACT) c C FACTOR N INTO ITS PRIME FACTORS, NFACT IN NUMBER. FOR EXAMPLE, C FOR N = 1960, NFACT = 6 AND IFACT(IF) = 2, 2, 2, 5, 7 AND 7. c implicit real*8(a-h,o-z) DIMENSION IFACT(1) IF=0 NPART=N DO 50 ID=1,N,2 IDIV=ID IF (ID-1) 10,10,20 10 IDIV=2 20 IQUOT=NPART/IDIV IF (NPART-IDIV*IQUOT) 40,30,40 30 IF=IF+1 IFACT(IF)=IDIV NPART=IQUOT GO TO 20 40 IF (IQUOT-IDIV) 60,60,50 50 CONTINUE 60 IF (NPART-1) 80,80,70 70 IF=IF+1 IFACT(IF)=NPART 80 NFACT=IF RETURN END c------------------------------------------------------------------------- SUBROUTINE DSMFAC(IFACT,NFACT,ISYM,IFSYM,NFSYM,ICENT,IFCNT,NFCNT) c C REARRANGE THE PRIME FACTORS OF N INTO A SQUARE AND A NON- C SQUARE. N = ISYM*ICENT*ISYM, WHERE ICENT IS SQUARE-FREE. C ISYM = IFSYM(1)*...*IFSYM(NFSYM), EACH A PRIME FACTOR. C ICENT = IFCNT(1)*...*IFCNT(NFCNT), EACH A PRIME FACTOR. C FOR EXAMPLE, N = 1960 = 14*10*14. THEN ISYM = 14, ICENT = 10, C NFSYM = 2, NFCNT = 2, NFACT = 6, IFSYM(IFS) = 2, 7, IFCNT(IFC) = C 2, 5 AND IFACT(IF) = 2, 7, 2, 5, 7, 2. c implicit real*8(a-h,o-z) DIMENSION IFSYM(1), IFCNT(1), IFACT(1) ISYM=1 ICENT=1 IFS=0 IFC=0 IF=1 10 IF (IF-NFACT) 20,40,50 20 IF (IFACT(IF)-IFACT(IF+1)) 40,30,40 30 IFS=IFS+1 IFSYM(IFS)=IFACT(IF) ISYM=IFACT(IF)*ISYM IF=IF+2 GO TO 10 40 IFC=IFC+1 IFCNT(IFC)=IFACT(IF) ICENT=IFACT(IF)*ICENT IF=IF+1 GO TO 10 50 NFSYM=IFS NFCNT=IFC NFSM2=2*NFSYM NFACT=2*NFSYM+NFCNT IF (NFCNT) 80,80,60 60 NFSM2=NFSM2+1 IFSYM(NFSYM+1)=ICENT DO 70 IFC=1,NFCNT IF=NFSYM+IFC 70 IFACT(IF)=IFCNT(IFC) 80 IF (NFSYM) 110,110,90 90 DO 100 IFS=1,NFSYM IFSCJ=NFSM2+1-IFS IFSYM(IFSCJ)=IFSYM(IFS) IFACT(IFS)=IFSYM(IFS) IFCNJ=NFACT+1-IFS 100 IFACT(IFCNJ)=IFSYM(IFS) 110 NFSYM=NFSM2 RETURN END c----------------------------------------------------------------------- SUBROUTINE DSYMRV(DATA,NPREV,N,NREM,IFACT,NFACT) c C SHUFFLE THE DATA ARRAY BY REVERSING THE DIGITS OF ONE INDEX. C DIMENSION DATA(NPREV,N,NREM) C REPLACE DATA(I1,I2,I3) BY DATA(I1,I2REV,I3) FOR ALL I1 FROM 1 TO C NPREV, I2 FROM 1 TO N AND I3 FROM 1 TO NREM. I2REV-1 IS THE C INTEGER WHOSE DIGIT REPRESENTATION IN THE MULTI-RADIX NOTATION C OF FACTORS IFACT(IF) IS THE REVERSE OF THE REPRESENTATION OF I2-1. C FOR EXAMPLE, IF ALL IFACT(IF) = 2, I2-1 = 11001, I2REV-1 = 10011. C THE FACTORS MUST BE SYMMETRICALLY ARRANGED, I.E., IFACT(IF) = C IFACT(NFACT+1-IF). c implicit real*8(a-h,o-z) DIMENSION DATA(1), IFACT(1) c IF (NFACT-1) 80,80,10 10 IP0=2 IP1=IP0*NPREV IP4=IP1*N IP5=IP4*NREM I4REV=1 DO 70 I4=1,IP4,IP1 IF (I4-I4REV) 20,40,40 20 I1MAX=I4+IP1-IP0 DO 30 I1=I4,I1MAX,IP0 DO 30 I5=I1,IP5,IP4 I5REV=I4REV+I5-I4 TEMPR=DATA(I5) TEMPI=DATA(I5+1) DATA(I5)=DATA(I5REV) DATA(I5+1)=DATA(I5REV+1) DATA(I5REV)=TEMPR 30 DATA(I5REV+1)=TEMPI 40 IP3=IP4 DO 60 IF=1,NFACT IP2=IP3/IFACT(IF) I4REV=I4REV+IP2 IF (I4REV-IP3) 70,70,50 50 I4REV=I4REV-IP3 60 IP3=IP2 70 CONTINUE 80 RETURN END c------------------------------------------------------------------------- SUBROUTINE DASMRV(DATA,NPREV,N,NREM,IFACT,NFACT,WORK) c C SHUFFLE THE DATA ARRAY BY REVERSING THE DIGITS OF ONE INDEX. C THE OPERATION IS THE SAME AS IN SYMRV, EXCEPT THAT THE FACTORS C NEED NOT BE SYMMETRICALLY ARRANGED, I.E., GENERALLY IFACT(IF) NOT= C IFACT(NFACT+1-IF). CONSEQUENTLY, A WORK ARRAY OF LENGTH N IS C NEEDED. c implicit real*8(a-h,o-z) DIMENSION DATA(1), WORK(1), IFACT(1) IF (NFACT-1) 60,60,10 10 IP0=2 IP1=IP0*NPREV IP4=IP1*N IP5=IP4*NREM DO 50 I1=1,IP1,IP0 DO 50 I5=I1,IP5,IP4 IWORK=1 I4REV=I5 I4MAX=I5+IP4-IP1 DO 40 I4=I5,I4MAX,IP1 WORK(IWORK)=DATA(I4REV) WORK(IWORK+1)=DATA(I4REV+1) IP3=IP4 DO 30 IF=1,NFACT IP2=IP3/IFACT(IF) I4REV=I4REV+IP2 IF (I4REV-IP3-I5) 40,20,20 20 I4REV=I4REV-IP3 30 IP3=IP2 40 IWORK=IWORK+IP0 IWORK=1 DO 50 I4=I5,I4MAX,IP1 DATA(I4)=WORK(IWORK) DATA(I4+1)=WORK(IWORK+1) 50 IWORK=IWORK+IP0 60 RETURN END c------------------------------------------------------------------------- SUBROUTINE DCOOL (DATA,NPREV,N,NREM,ISIGN,IFACT,WORK) c C FOURIER TRANSFORM OF LENGTH N. IN PLACE COOLEY-TUKEY METHOD, C DIGIT-REVERSED TO NORMAL ORDER, SANDE-TUKEY FACTORING (2). C DIMENSION DATA(NPREV,N,NREM) C COMPLEX DATA C DATA(I1,J2,I3) = SUM(DATA(I1,I2,I3)*EXP(ISIGN*2*PI*I*((I2-1)* C (J2-1)/N))), SUMMED OVER I2 = 1 TO N FOR ALL I1 FROM 1 TO NPREV, C J2 FROM 1 TO N AND I3 FROM 1 TO NREM. THE FACTORS OF N ARE GIVEN C IN ANY ORDER IN ARRAY IFACT. FACTORS OF TWO ARE DONE IN PAIRS C AS MUCH AS POSSIBLE (FOURIER TRANSFORM OF LENGTH FOUR), FACTORS OF C THREE ARE DONE SEPARATELY, AND ALL FACTORS FIVE OR HIGHER C ARE DONE BY GOERTZEL'S ALGORITHM (4). c implicit real*8(a-h,o-z) DIMENSION DATA(1), WORK(1), IFACT(1) pi = 4.0d0*datan(1.d0) c TWOPI=6.283185307*FLOAT(ISIGN) TWOPI=2.d0*pi*DFLOAT(ISIGN) IP0=2 IP1=IP0*NPREV IP4=IP1*N IP5=IP4*NREM IF=0 IP2=IP1 10 IF (IP2-IP4) 20,240,240 20 IF=IF+1 IFCUR=IFACT(IF) IF (IFCUR-2) 60,30,60 30 IF (4*IP2-IP4) 40,40,60 40 IF (IFACT(IF+1)-2) 60,50,60 50 IF=IF+1 IFCUR=4 60 IP3=IP2*IFCUR THETA=TWOPI/DFLOAT(IFCUR) SINTH=dsin(THETA/2.d0) ROOTR=-2.d0*SINTH*SINTH C COS(THETA)-1, FOR ACCURACY. ROOTI=dsin(THETA) THETA=TWOPI/DFLOAT(IP3/IP1) SINTH=dsin(THETA/2.d0) WSTPR=-2.d0*SINTH*SINTH WSTPI=dsin(THETA) WR=1.d0 WI=0.d0 DO 230 I2=1,IP2,IP1 IF (IFCUR-4) 70,70,210 70 IF ((I2-1)*(IFCUR-2)) 240,90,80 80 W2R=WR*WR-WI*WI W2I=2.d0*WR*WI W3R=W2R*WR-W2I*WI W3I=W2R*WI+W2I*WR 90 I1MAX=I2+IP1-IP0 DO 200 I1=I2,I1MAX,IP0 DO 200 I5=I1,IP5,IP3 J0=I5 J1=J0+IP2 J2=J1+IP2 J3=J2+IP2 IF (I2-1) 140,140,100 100 IF (IFCUR-3) 130,120,110 C APPLY THE PHASE SHIFT FACTORS 110 TEMPR=DATA(J3) DATA(J3)=W3R*TEMPR-W3I*DATA(J3+1) DATA(J3+1)=W3R*DATA(J3+1)+W3I*TEMPR TEMPR=DATA(J2) DATA(J2)=WR*TEMPR-WI*DATA(J2+1) DATA(J2+1)=WR*DATA(J2+1)+WI*TEMPR TEMPR=DATA(J1) DATA(J1)=W2R*TEMPR-W2I*DATA(J1+1) DATA(J1+1)=W2R*DATA(J1+1)+W2I*TEMPR GO TO 140 120 TEMPR=DATA(J2) DATA(J2)=W2R*TEMPR-W2I*DATA(J2+1) DATA(J2+1)=W2R*DATA(J2+1)+W2I*TEMPR 130 TEMPR=DATA(J1) DATA(J1)=WR*TEMPR-WI*DATA(J1+1) DATA(J1+1)=WR*DATA(J1+1)+WI*TEMPR 140 IF (IFCUR-3) 150,160,170 C DO A FOURIER TRANSFORM OF LENGTH TWO 150 TEMPR=DATA(J1) TEMPI=DATA(J1+1) DATA(J1)=DATA(J0)-TEMPR DATA(J1+1)=DATA(J0+1)-TEMPI DATA(J0)=DATA(J0)+TEMPR DATA(J0+1)=DATA(J0+1)+TEMPI GO TO 200 C DO A FOURIER TRANSFORM OF LENGTH THREE 160 SUMR=DATA(J1)+DATA(J2) SUMI=DATA(J1+1)+DATA(J2+1) TEMPR=DATA(J0)-0.5d0*SUMR TEMPI=DATA(J0+1)-.5d0*SUMI DATA(J0)=DATA(J0)+SUMR DATA(J0+1)=DATA(J0+1)+SUMI DIFR=ROOTI*(DATA(J2+1)-DATA(J1+1)) DIFI=ROOTI*(DATA(J1)-DATA(J2)) DATA(J1)=TEMPR+DIFR DATA(J1+1)=TEMPI+DIFI DATA(J2)=TEMPR-DIFR DATA(J2+1)=TEMPI-DIFI GO TO 200 C DO A FOURIER TRANSFORM OF LENGTH FOUR (FROM BIT REVERSED ORDER) 170 T0R=DATA(J0)+DATA(J1) T0I=DATA(J0+1)+DATA(J1+1) T1R=DATA(J0)-DATA(J1) T1I=DATA(J0+1)-DATA(J1+1) T2R=DATA(J2)+DATA(J3) T2I=DATA(J2+1)+DATA(J3+1) T3R=DATA(J2)-DATA(J3) T3I=DATA(J2+1)-DATA(J3+1) DATA(J0)=T0R+T2R DATA(J0+1)=T0I+T2I DATA(J2)=T0R-T2R DATA(J2+1)=T0I-T2I IF (ISIGN) 180,180,190 180 T3R=-T3R T3I=-T3I 190 DATA(J1)=T1R-T3I DATA(J1+1)=T1I+T3R DATA(J3)=T1R+T3I DATA(J3+1)=T1I-T3R 200 CONTINUE GO TO 220 C DO A FOURIER TRANSFORM OF LENGTH FIVE OR MORE 210 CALL DGOERT(DATA(I2),NPREV,IP2/IP1,IFCUR,IP5/IP3,WORK,WR,WI,ROOTR, $ROOTI) 220 TEMPR=WR WR=WSTPR*TEMPR-WSTPI*WI+TEMPR 230 WI=WSTPR*WI+WSTPI*TEMPR+WI IP2=IP3 GO TO 10 240 RETURN END c------------------------------------------------------------------------- SUBROUTINE DGOERT(DATA,NPREV,IPROD,IFACT,IREM,WORK,WMINR,WMINI, $ ROOTR,ROOTI) c C PHASE-SHIFTED FOURIER TRANSFORM OF LENGTH IFACT BY THE GOERTZEL C ALGORITHM (4). IFACT MUST BE ODD AND AT LEAST 5. FURTHER SPEED C IS GAINED BY COMPUTING TWO TRANSFORM VALUES AT THE SAME TIME. C DIMENSION DATA(NPREV,IPROD,IFACT,IREM) C DATA(I1,1,J3,I5) = SUM(DATA(I1,1,I3,I5) * W**(I3-1)), SUMMED C OVER I3 = 1 TO IFACT FOR ALL I1 FROM 1 TO NPREV, J3 FROM 1 TO C IFACT AND I5 FROM 1 TO IREM. C W = WMIN * EXP(ISIGN*2*PI*I*(J3-1)/IFACT). c implicit real*8(a-h,o-z) DIMENSION DATA(1), WORK(1) IP0=2 IP1=IP0*NPREV IP2=IP1*IPROD IP3=IP2*IFACT IP5=IP3*IREM IF (WMINI) 10,40,10 C APPLY THE PHASE SHIFT FACTORS 10 WR=WMINR WI=WMINI I3MIN=1+IP2 DO 30 I3=I3MIN,IP3,IP2 I1MAX=I3+IP1-IP0 DO 20 I1=I3,I1MAX,IP0 DO 20 I5=I1,IP5,IP3 TEMPR=DATA(I5) DATA(I5)=WR*TEMPR-WI*DATA(I5+1) 20 DATA(I5+1)=WR*DATA(I5+1)+WI*TEMPR TEMPR=WR WR=WMINR*TEMPR-WMINI*WI 30 WI=WMINR*WI+WMINI*TEMPR 40 DO 90 I1=1,IP1,IP0 DO 90 I5=I1,IP5,IP3 C STRAIGHT SUMMATION FOR THE FIRST TERM SUMR=0.d0 SUMI=0.d0 I3MAX=I5+IP3-IP2 DO 50 I3=I5,I3MAX,IP2 SUMR=SUMR+DATA(I3) 50 SUMI=SUMI+DATA(I3+1) WORK(1)=SUMR WORK(2)=SUMI WR=ROOTR+1.d0 WI=ROOTI IWMIN=1+IP0 IWMAX=IP0*((IFACT+1)/2)-1 DO 80 IWORK=IWMIN,IWMAX,IP0 TWOWR=WR+WR I3=I3MAX OLDSR=0.d0 OLDSI=0.d0 SUMR=DATA(I3) SUMI=DATA(I3+1) I3=I3-IP2 60 TEMPR=SUMR TEMPI=SUMI SUMR=TWOWR*SUMR-OLDSR+DATA(I3) SUMI=TWOWR*SUMI-OLDSI+DATA(I3+1) OLDSR=TEMPR OLDSI=TEMPI I3=I3-IP2 IF (I3-I5) 70,70,60 C IN A FOURIER TRANSFORM THE W CORRESPONDING TO THE POINT AT K C IS THE CONJUGATE OF THAT AT IFACT-K (THAT IS, EXP(TWOPI*I* C K/IFACT) = CONJ(EXP(TWOPI*I*(IFACT-K)/IFACT))). SINCE THE C MAIN LOOP OF GOERTZELS ALGORITHM IS INDIFFERENT TO THE IMAGINARY C PART OF W, IT NEED BE SUPPLIED ONLY AT THE END. 70 TEMPR=-WI*SUMI TEMPI=WI*SUMR SUMR=WR*SUMR-OLDSR+DATA(I3) SUMI=WR*SUMI-OLDSI+DATA(I3+1) WORK(IWORK)=SUMR+TEMPR WORK(IWORK+1)=SUMI+TEMPI IWCNJ=IP0*(IFACT+1)-IWORK WORK(IWCNJ)=SUMR-TEMPR WORK(IWCNJ+1)=SUMI-TEMPI C SINGLETON'S RECURSION, FOR ACCURACY AND SPEED (5). TEMPR=WR WR=WR*ROOTR-WI*ROOTI+WR 80 WI=TEMPR*ROOTI+WI*ROOTR+WI IWORK=1 DO 90 I3=I5,I3MAX,IP2 DATA(I3)=WORK(IWORK) DATA(I3+1)=WORK(IWORK+1) 90 IWORK=IWORK+IP0 RETURN END c------------------------------------------------------------------------- SUBROUTINE DFIXRL(DATA,N,NREM,ISIGN,IFORM) c C FOR IFORM = 0, CONVERT THE TRANSFORM OF A DOUBLED-UP REAL ARRAY, C CONSIDERED COMPLEX, INTO ITS TRUE TRANSFORM. SUPPLY ONLY THE C FIRST HALF OF THE COMPLEX TRANSFORM, AS THE SECOND HALF HAS C CONJUGATE SYMMETRY. FOR IFORM = -1, CONVERT THE FIRST HALF C OF THE TRUE TRANSFORM INTO THE TRANSFORM OF A DOUBLED-UP REAL C ARRAY. N MUST BE EVEN. C USING COMPLEX NOTATION AND SUBSCRIPTS STARTING AT ZERO, THE C TRANSFORMATION IS-- C DIMENSION DATA(N,NREM) C ZSTP = EXP(ISIGN*2*PI*I/N) C DO 10 I2=0,NREM-1 C DATA(0,I2) = CONJ(DATA(0,I2))*(1+I) C DO 10 I1=1,N/4 C Z = (1+(2*IFORM+1)*I*ZSTP**I1)/2 C I1CNJ = N/2-I1 C DIF = DATA(I1,I2)-CONJ(DATA(I1CNJ,I2)) C TEMP = Z*DIF C DATA(I1,I2) = (DATA(I1,I2)-TEMP)*(1-IFORM) C 10 DATA(I1CNJ,I2) = (DATA(I1CNJ,I2)+CONJ(TEMP))*(1-IFORM) C IF I1=I1CNJ, THE CALCULATION FOR THAT VALUE COLLAPSES INTO C A SIMPLE CONJUGATION OF DATA(I1,I2). c implicit real*8(a-h,o-z) DIMENSION DATA(1) pi = 4.d0*datan(1.d0) c TWOPI=6.283185307*FLOAT(ISIGN) TWOPI=2.d0*pi*DFLOAT(ISIGN) IP0=2 IP1=IP0*(N/2) IP2=IP1*NREM IF (IFORM) 10,70,70 C PACK THE REAL INPUT VALUES (TWO PER COLUMN) 10 J1=IP1+1 DATA(2)=DATA(J1) IF (NREM-1) 70,70,20 20 J1=J1+IP0 I2MIN=IP1+1 DO 60 I2=I2MIN,IP2,IP1 DATA(I2)=DATA(J1) J1=J1+IP0 IF (N-2) 50,50,30 30 I1MIN=I2+IP0 I1MAX=I2+IP1-IP0 DO 40 I1=I1MIN,I1MAX,IP0 DATA(I1)=DATA(J1) DATA(I1+1)=DATA(J1+1) 40 J1=J1+IP0 50 DATA(I2+1)=DATA(J1) 60 J1=J1+IP0 70 DO 80 I2=1,IP2,IP1 TEMPR=DATA(I2) DATA(I2)=DATA(I2)+DATA(I2+1) 80 DATA(I2+1)=TEMPR-DATA(I2+1) IF (N-2) 200,200,90 90 THETA=TWOPI/DFLOAT(N) SINTH=dsin(THETA/2.d0) ZSTPR=-2.d0*SINTH*SINTH ZSTPI=dsin(THETA) ZR=(1.d0-ZSTPI)/2.d0 ZI=(1.d0+ZSTPR)/2.d0 IF (IFORM) 100,110,110 100 ZR=1.d0-ZR ZI=-ZI 110 I1MIN=IP0+1 I1MAX=IP0*(N/4)+1 DO 190 I1=I1MIN,I1MAX,IP0 DO 180 I2=I1,IP2,IP1 I2CNJ=IP0*(N/2+1)-2*I1+I2 IF (I2-I2CNJ) 150,120,120 120 IF (ISIGN*(2*IFORM+1)) 130,140,140 130 DATA(I2+1)=-DATA(I2+1) 140 IF (IFORM) 170,180,180 150 DIFR=DATA(I2)-DATA(I2CNJ) DIFI=DATA(I2+1)+DATA(I2CNJ+1) TEMPR=DIFR*ZR-DIFI*ZI TEMPI=DIFR*ZI+DIFI*ZR DATA(I2)=DATA(I2)-TEMPR DATA(I2+1)=DATA(I2+1)-TEMPI DATA(I2CNJ)=DATA(I2CNJ)+TEMPR DATA(I2CNJ+1)=DATA(I2CNJ+1)-TEMPI IF (IFORM) 160,180,180 160 DATA(I2CNJ)=DATA(I2CNJ)+DATA(I2CNJ) DATA(I2CNJ+1)=DATA(I2CNJ+1)+DATA(I2CNJ+1) 170 DATA(I2)=DATA(I2)+DATA(I2) DATA(I2+1)=DATA(I2+1)+DATA(I2+1) 180 CONTINUE TEMPR=ZR-.5d0 ZR=ZSTPR*TEMPR-ZSTPI*ZI+ZR 190 ZI=ZSTPR*ZI+ZSTPI*TEMPR+ZI C RECURSION SAVES TIME, AT A SLIGHT LOSS IN ACCURACY. IF AVAILABLE, C USE DOUBLE PRECISION TO COMPUTE ZR AND ZI. 200 IF (IFORM) 270,210,210 C UNPACK THE REAL TRANSFORM VALUES (TWO PER COLUMN) 210 I2=IP2+1 I1=I2 J1=IP0*(N/2+1)*NREM+1 GO TO 250 220 DATA(J1)=DATA(I1) DATA(J1+1)=DATA(I1+1) I1=I1-IP0 J1=J1-IP0 230 IF (I2-I1) 220,240,240 240 DATA(J1)=DATA(I1) DATA(J1+1)=0.d0 250 I2=I2-IP1 J1=J1-IP0 DATA(J1)=DATA(I2+1) DATA(J1+1)=0.d0 I1=I1-IP0 J1=J1-IP0 IF (I2-1) 260,260,230 260 DATA(2)=0.d0 270 RETURN END