c--period.f c--simple periodogram using fourt c------------------------------------------------------------------- program PERIOD implicit real*8(a-h,o-y),integer*4(i-n),logical(z) c c--change the dimensions between Unix and DOS c--Unix : ndim = 1100000, ndim2 = 2100000 (2^21 = 2,097,152) c--DOS ndim = 66000, ndim2 = 132000 (2^17 = 131,072) parameter(ndim=1100000,ndim2=2100000) complex*16 cdata(ndim),cwork(16384),ctemp dimension data(ndim2) character*60 title character*40 fname,ftemp character*20 fmt1,fmt2,dunit,tunit,funit character*1 ch equivalence (cdata,data) common // data,cwork c c--initialise data fmt2/'(2g15.7)'/,nwork/16384/ pi = 4.d0*datan(1.d0) rdeg = 180.d0/pi do 2 i = 1,21 np = 2**i if(np.le.ndim2) imax = i 2 continue 1 write(*,'(//a)') ' PERIOD - (mag,phase) FFT of a real data set' write(*,'(a)') ' Default responses are [ ]' 5 write(*,'(a)') & ' ----------------------------------------------------' c c--read data from file write(*,'(a$)') ' Input filename --> ' read(*,'(a)') fname inquire(file = fname,exist=zexist) if(.not.zexist) then write(*,'(a$)') ' file not found, retry? [y]/n --> ' read(*,'(a1)') ch if(ch.ne.'n') goto 5 stop endif open(1,file = fname) read(1,'(a)') title read(1,*) n,delt read(1,'(3a20)') fmt1,dunit,tunit write(*,'(2a)') ' title : ',title write(*,'(a,i10)') ' # points : ',n write(*,'(a,g12.5)') ' delta t : ',delt write(*,'(2a)') ' data units: ',dunit write(*,'(2a)') ' time units: ',tunit write(*,'(2a)') ' format : ',fmt1 read(1,fmt1) (data(i),i = 1,n) close(1) c c--energy and power of original signal in time domain sumtpo = 0.d0 do 6 i = 1,n sumtpo = sumtpo+data(i)*data(i)*delt 6 continue sumtpo = sumtpo/(n*delt) c c--remove mean if requested write(*,'(a$)') ' Remove mean? [y]/n --> ' read(*,'(a1)') ch if(ch.ne.'n') call rmean(data,n,dmean,0) c c--data window call window(data,n,fmt1,delt,tunit,iw,lcb) c c--energy and power of windowed signal in time domain sumtpw = 0.d0 do 26 i = 1,n sumtpw = sumtpw+data(i)*data(i)*delt 26 continue sumtpw = sumtpw/(n*delt) c c--check for next power of 2 do 10 index = 1,21 np = 2**index if(np.ge.n) goto 12 10 continue c c--pad with zeros and adjust n to np 12 write(*,'(2(a,i7),a,i3)')' Data length = ',n,'; np = ',np, & ' = 2**',index write(*,'(a,i2,a$)') & ' new index (0 = no pad; max = ',imax,')--> ' read(*,*) index index = min0(index,imax) if(index.eq.0) then np = n else np = 2**index do 13 i = n+1,np data(i) = 0.d0 13 continue write(*,'(a,i10,a,i10)') ' data padded from ',n,' to',np endif nyq = (np/2)+1 c c--call fourt for FFT c write(*,'(a$)') ' calling fourt ... ' call dfourt(data,np,1,-1,0,cwork,nwork) c write(*,'(a)') ' finished. ' c c--estimate power in frequency domain from squared amplitudes c--if zero or Nyquist frequency, count term only once, otherwise c--double the term at all other positive frequencies sumfe = 0.d0 do 27 k = 1,nyq delta = 2.0 if(k.eq.1.or.k.eq.nyq) delta = 1.0 s = cdata(k)*conjg(cdata(k)) sumfe = sumfe+s*delta 27 continue sumfe = sumfe*delt/np sumfp = sumfe/(n*delt) write(*,'(a,g13.5)') & ' Original time domain power = ',sumtpo, & ' Windowed time domain power = ',sumtpw, & ' Frequency domain power = ',sumfp c c--output file write(*,'(a$)') ' Normalise output amplitudes? y/[n] --> ' read(*,'(a1)') ch inorm = 0 if(ch.eq.'y') inorm = 1 call freq(tunit,funit) delf = 1.d0/(np*delt) c c--get output filename write(*,'(a$)') ' Output filename ( -> input.fft) --> ' read(*,'(a)') ftemp lft = lstr(ftemp,40) if(lft.gt.1) fname = ftemp call fnext(fname,'fft') write(*,'(2a)') ' writing ',fname c c--write output open(2,file = fname) write(2,'(a)') title write(2,'(3i10,g12.5,2i10)') inorm,n,nyq,delf,iw,lcb write(2,'(3a20)') fmt2,dunit,funit c c--calculate and write Fourier amplitudes and phases do 30 k = 1,nyq ctemp = cdata(k) delta = 2.0 if(k.eq.1.or.k.eq.nyq) delta = 1.0 c--adjust normalisation if required if(inorm.eq.1) ctemp = ctemp*delta/float(n) d1 = real(ctemp) d2 = dimag(ctemp) amp = cdabs(ctemp) if(dabs(d1).gt.1.d-32) then phase = rdeg*datan2(d2,d1) else phase = 180.d0 endif if(fname.ne.'nul') write(2,fmt2) amp,phase 30 continue if(fname.ne.'nul') close(2) write(*,'(a$)') ' Another file? y/[n] --> ' read(*,'(a1)') ch if(ch.ne.'y') stop goto 5 c end c-------------------------------------------------------------------------- subroutine freq(tunit,funit) c c--creates frequency units out of time units c--adds 'cycles/' to beginning of tunit implicit real*8(a-h,o-y),integer*4(i-n),logical(z) character*20 tunit,funit,temp1,temp2 character*1 dum1(20),dum2(20) equivalence (temp1,dum1),(temp2,dum2) data temp2/'cycles/ '/ c temp1 = tunit l1 = lstr(temp1,-20) do 10 i = 1,l1 dum2(7+i) = dum1(i) 10 continue funit = temp2 c return end c--------------------------------------------------------------------------- subroutine window(data,n,fmt1,delt,tunit,iw,lcb) c c apply a data window c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) dimension data(1) character*40 fname character*20 fmt1,fmt,dunit,tunit character*1 ans c c--initialise data fmt,dunit/'(g15.7)','amplitude'/ pi = 4.d0*datan(1.d0) c c--select window write(*,'(a$)') &' Data window : 1 box, 2 cosine, 3 Hann --> ' read(*,*) iw c--quit now if boxcar if(iw.eq.1) return c--find percentage length of cosine bell if(iw.eq.2) then write(*,'(a$)') ' cosine bell length (e.g. 10) in % --> ' read(*,*) lcb endif c c--save options zans1 = .false. write(*,'(a$)') ' save data window? y/[n] --> ' read(*,'(a1)') ans if(ans.eq.'y'.or.ans.eq.'Y') zans1 = .true. zans2 = .false. write(*,'(a$)') ' save windowed data? y/[n] --> ' read(*,'(a1)') ans if(ans.eq.'y'.or.ans.eq.'Y') zans2 = .true. c c--open window files if required if(zans1) then if(iw.eq.2) then open(9,file = 'cbell.dat') write(9,'(a,i3,A)') ' cosine bell, length ',lcb,'%' else open(9,file = 'hanning.dat') write(9,'(a)') ' Hanning Window' endif write(9,'(i10,g12.5)') n,delt write(9,'(3a20)') fmt,dunit,tunit endif c c--for the cosine and hanning windows, first remove mean from the data call rmean(data,n,dmean,0) c c--now do the windowing if(iw.eq.2) then c--cosine bell lcos = nint(float(n*lcb)/200.0) pil = pi/lcos do 40 i = 1,n if(i.gt.(lcos+1)) goto 10 w = 0.5d0*(1.d0-dcos(pil*(i-1))) goto 20 10 w = 1.d0 if(i.le.(n-lcos+1)) goto 20 w = 0.5d0*(1.d0-dcos(pil*(n+1-i))) 20 data(i) = data(i)*w if(zans1) write(9,fmt) w 40 continue else c--Hanning do 50 i = 1,n t = float(2*i-2-n)/float(n) w = 0.5d0*(1.d0+dcos(pi*t)) data(i) = data(i)*w if(zans1) write(9,fmt) w 50 continue endif close(9) c c--compensate for energy loss due to window write(*,'(a$)') ' compensate for window? [y]/n --> ' read(*,'(a1)') ans if(ans.ne.'n') then if(iw.eq.2) then factor = 1.d0 - float(lcb)/160.d0 else factor = 0.375 endif factor = dsqrt(1.d0/factor) f2 = factor*factor write(*,'(a,g13.5)') & ' window correction factor (power) = ',f2 do 60 i = 1,n data(i) = factor*data(i) 60 continue endif c c--restore data mean for cosine and hanning windows call rmean(data,n,dmean,1) c c--maybe write windowed data if(zans2) then write(*,'(a$)') ' filename for windowed data --> ' read(*,'(a)') fname open(9,file = fname) write(9,'(a)') ' Windowed Data' write(9,'(i10,g12.5)') n,delt write(9,'(3a20)') fmt1,dunit,tunit write(9,fmt1) (data(i),i=1,n) close(9) endif c 999 return end c---------------------------------------------------------------------------- subroutine rmean(y,ly,ymean,k) c c rmean removes (k=0) and restores (k=1) mean of an array y(i) c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) dimension y(1) c c--if k=0, first calculate ymean if(k.eq.0) then sum = 0.0 do 10 i = 1,ly sum = sum+y(i) 10 continue ymean = sum/float(ly) fact = -1.0 else fact = 1.0 endif c c--subtract (k=0)/add (k=1) the mean as required do 20 i= 1,ly y(i) = y(i)+fact*ymean 20 continue c return end c--------------------------------------------------------------------------- subroutine fnext(fname,ext) c c changes filename extension to '.ext' c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) character*40 fname,fnamec character*3 ext,extc character*1 dum1(40),dum3(3) equivalence (fnamec,dum1),(extc,dum3) c c original length of filename is 'len0' c find length of extension 'len3' by searching last 4 characters c fnamec = fname extc = ext len0 = lstr(fnamec,-40) len3 = 0 ltmp = 0 do 5 i = 1,4 l = len0-i+1 if(l.eq.0) goto 6 if(dum1(l).eq.'.') then len3 = ltmp+1 goto 6 endif ltmp = ltmp+1 5 continue c c find length of main segment of the filename, back to a '\' c if it exists as part of the filename c 6 len2 = 0 do 10 i = 1,8 l = len0-len3-i+1 if(l.eq.0.or.dum1(l).eq.'\\') goto 15 len2 = len2+1 10 continue c c find length of first segment (before any '\') c new overall length of filename is len1+(new)len2+'.ext' c 15 len1 = len0-len3-len2 len = len1+len2+4 c c now assign last 4 characters of new filename c dum1(len-3) = '.' dum1(len-2) = dum3(1) dum1(len-1) = dum3(2) dum1(len) = dum3(3) fname = fnamec c return end c--------------------------------------------------------------- function lstr(name,n) c c Function to return the length of a string. On input the string c 'name' has a nominal length n; if n is negative, blank characters c at the beginning of 'name' are first truncated. c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) character*1 name(1) c c--zero length case first if(n.eq.0) then lstr = 0 return endif k = 0 ld = n c--ld is negative, we truncate leading blanks if(ld.lt.0) then ld = -ld j = 0 10 j = j+1 if(name(j).eq.' ') goto 10 nb = j-1 c--there are nb blanks to be erased in front of 'name' if(nb.gt.0) then do 20 i = nb+1,ld name(i-nb) = name(i) name(i) = ' ' 20 continue c--reduce length of string by number of blanks ld = ld-nb endif endif c--ld is non-zero, find length of string c--first set true length equal to nominal length k = ld+1 c--and decrease k until a non-blank character appears or beginning c--of string is reached 40 k = k-1 if(k.gt.0.and.name(k).eq.' ') goto 40 lstr = k 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-y),integer*4(i-n),logical(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-y),integer*4(i-n),logical(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-y),integer*4(i-n),logical(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-y),integer*4(i-n),logical(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-y),integer*4(i-n),logical(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-y),integer*4(i-n),logical(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-y),integer*4(i-n),logical(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),integer*4(i-n) 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