C @(#)fouri.for 17.1.1.1 (ESO-DMD) 01/25/02 17:40:57 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM FOURI C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program FOURI version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C Fourier transformation C C.PURPOSE C calculate the Fourier transformation of given real + imaginary frame C C.ALGORITHM C use the FFT algorithm as in E.O. Brigham: "The Fast Fourier Transform", C Prentice-Hall, Inc., 1974, page 164 C C.INPUT/OUTPUT C the following keywords are used: C C P1/C/1/60 name of real data frame C P2/C/1/60 name of imaginary data frame C if IN_B is "+", only real input data C P3/C/1/60 name of frame holding real part of FT C P4/C/1/60 name of frame holding imaginary part of FT C C CACTIO/C/1/3 '???' = normal way C 'INV' = inverse FFT C 'POW' = return also power spectrum C 'FRE' = use frequencies C (center shifted to the middle) C 'FPO' = FRE + POW combined C 'FIN' = FRE + INV combined C TMPBUF/C/1/3 (1:1) = spare C (2:2) = I for inverse, P for power stuff C (3:3) will be set to Y, if expansion C to a power2 frame took place or not C C P5/C/1/60 if action = POWER or FPOWER, C name of power spectrum C C.VERSIONS C C 011126 last modif C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,M,N,NAXISA,NAXISB,NN,NOPIX INTEGER IMNOA,IMNOB,IMNOC,IMNOD,IMNOP INTEGER*8 PNTRA,PNTRB,PNTRC,PNTRD,PNTRP INTEGER*8 DI,DR,CT,DT,TS,TC,SWAP1,SWAP2 INTEGER IMNODI,IMNODR,IMNOTC,IMNOTS INTEGER SIZE,DSTAT,STAT INTEGER NLOG(3),NPIXA(3),NPIXB(3),NPIXF(3),EXTFLG,POW2(22) INTEGER LOW(3) INTEGER UNI(1),NULO,MADRID(1) INTEGER EC,EL,ED C CHARACTER*60 WFRAME,FRAMEA,FRAMEB,FRAMEC,FRAMED,FRAMEP CHARACTER IDENT*72,CUNIT*48,ACTION*3,TMP*3 C DOUBLE PRECISION STARTA(3),STARTB(3),STEPA(3),STEPB(3) DOUBLE PRECISION DBUF(6) C REAL CUTS(4),RC C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA FRAMEA /' '/, FRAMEB /' '/ DATA FRAMEC /' '/, FRAMED /' '/ DATA IDENT /' '/, CUNIT /' '/ DATA ACTION /' '/ DATA NPIXA /3*1/, NPIXB /3*1/, NPIXF /3*1/ DATA STARTA /3*0.D0/, STARTB /3*0.D0/ DATA STEPA /3*1.D0/, STEPB /3*1.D0/ DATA EXTFLG /0/ DATA POW2 /1,2,4,8,16,32,64,128,256,512,1024,2048,4096, + 8192,16384,32768,65536,131072,262144,524288, + 1048576,2097152/ DATA RC /0.0/ DATA LOW /3*1/ NLOG /3*0/ DATA CT /1/, DT /1/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get into MIDAS environment CALL STSPRO('FOURI') C C get action CALL STKRDC('CACTIO',1,1,3,IAV,ACTION,UNI,NULO,STAT) CALL UPCAS(ACTION,ACTION) CALL STKRDC('TMPBUF',1,1,3,IAV,TMP,UNI,NULO,STAT) C C get data frames involved CALL STKRDC('P1',1,1,60,IAV,FRAMEA,UNI,NULO,STAT) CALL STKRDC('P2',1,1,60,IAV,FRAMEB,UNI,NULO,STAT) CALL STKRDC('P3',1,1,60,IAV,FRAMEC,UNI,NULO,STAT) CALL STKRDC('P4',1,1,60,IAV,FRAMED,UNI,NULO,STAT) C C map real input frame CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 3,NAXISA,NPIXA,STARTA,STEPA, + IDENT,CUNIT,PNTRA,IMNOA,STAT) C IF (NAXISA.GT.2) + CALL STETER(1,'only 2-dim frames supported currently...') C C get imaginary part IF (FRAMEB(1:1).NE.'+') THEN CALL STIGET(FRAMEB,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 3,NAXISB,NPIXB,STARTB,STEPB, + IDENT,CUNIT,PNTRB,IMNOB,STAT) IF (NAXISB.NE.NAXISA) GOTO 9000 C DO 200, N=1,NAXISA !frames have to be identical IF (NPIXA(N).NE.NPIXB(N)) GOTO 9000 IF (ACTION(1:1).EQ.'F') THEN IF (ABS(STEPA(N)-STEPB(N)).GT.10.D-15) GOTO 9000 IF (ABS(STARTA(N)-STARTB(N)).GT.10.D-15) GOTO 9000 ENDIF 200 CONTINUE ENDIF C C for inverse FFT check descriptor FFT_SPECIAL with original START, STEP C and eventual NPIX before expansion C DBUF(1) = STARTA(1) DBUF(2) = STARTA(2) DBUF(3) = STEPA(1) DBUF(4) = STEPA(2) DBUF(5) = DBLE(NPIXA(1)) DBUF(6) = DBLE(NPIXA(2)) C IF (TMP(2:2).EQ.'I') THEN CALL STECNT('GET',EC,EL,ED) CALL STECNT('PUT',1,0,0) !test, if FFT_SPECIAL is there CALL STDRDD(IMNOA,'FFT_SPECIAL',1,6,IAV,DBUF,UNI,NULO,DSTAT) CALL STECNT('PUT',EC,EL,ED) IF (DSTAT.EQ.0) THEN STARTA(1) = DBUF(1) STARTA(2) = DBUF(2) STEPA(1) = DBUF(3) STEPA(2) = DBUF(4) ENDIF ELSE DSTAT = 1 !indicate error... ENDIF C C calculate log of dimension(s) DO 400, M=1,NAXISA NN = NPIXA(M) NPIXF(M) = NN DO 300, N=1,22 IF (NN.LE.POW2(N)) THEN IF (NN.LT.POW2(N)) THEN NPIXF(M) = POW2(N) EXTFLG = 1 ENDIF NLOG(M) = N - 1 GOTO 400 !jump out of loop ENDIF 300 CONTINUE CALL STETER(3,'x- or y-dimension larger than 2097152 ...') 400 CONTINUE C C calculate size IF (NAXISA.EQ.1) THEN !max. 2-dim frame... SIZE = NPIXF(1) ELSE SIZE = NPIXF(1) * NPIXF(2) CALL STFXMP(2*SIZE,D_R4_FORMAT,CT,STAT) DT = CT + SIZE ENDIF C C see, if we have to expand real input frame IF (EXTFLG.EQ.1) THEN CALL STFCRE('wworrkfr',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + SIZE,IMNOC,STAT) CALL STFMAP(IMNOC,F_X_MODE,1,SIZE,IAV,PNTRC,STAT) CALL CONFIL(RC,MADRID(PNTRC),SIZE) !fill with zero CALL COPWND(MADRID(PNTRA),NPIXA,MADRID(PNTRC),NPIXF, + LOW,LOW,NPIXA) !and insert input frame CALL STFCLO(IMNOA,STAT) IMNOA = IMNOC PNTRA = PNTRC C IF (FRAMEB(1:1).NE.'+') THEN !also expand imaginary file CALL STFCRE('wworrkfi',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + SIZE,IMNOC,STAT) CALL STFMAP(IMNOC,F_X_MODE,1,SIZE,IAV,PNTRC,STAT) CALL CONFIL(RC,MADRID(PNTRC),SIZE) !fill with zero CALL COPWND(MADRID(PNTRB),NPIXA,MADRID(PNTRC),NPIXF, + LOW,LOW,NPIXA) !and insert CALL STFCLO(IMNOB,STAT) IMNOB = IMNOC PNTRB = PNTRC ENDIF ENDIF C C update start, step for forward FFT IF (TMP(2:2).NE.'I') THEN DO 2400, N=1,NAXISA STEPA(N) = ABS( 1.D0 / (DBLE(NPIXF(N)) * STEPA(N)) ) IF (ACTION(1:1).EQ.'F') THEN STARTA(N) = - ( (NPIXF(N)/2) * STEPA(N) ) ELSE STARTA(N) = 0.0D0 ENDIF 2400 CONTINUE ENDIF C C create + map result frames IF (EXTFLG.EQ.1) THEN WFRAME(1:) = 'exp_fftr ' ELSE WFRAME(1:) = FRAMEC(1:) ENDIF CALL STIPUT(WFRAME,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISA,NPIXF,STARTA,STEPA, + IDENT,CUNIT,PNTRC,IMNOC,STAT) IF (EXTFLG.EQ.1) THEN WFRAME(1:) = 'exp_ffti ' ELSE WFRAME(1:) = FRAMED(1:) ENDIF CALL STIPUT(WFRAME,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISA,NPIXF,STARTA,STEPA, + IDENT,CUNIT,PNTRD,IMNOD,STAT) C C allocate memory for working buffer DREAL + DIMAG C as well as for sine + cosine look-up tables (all double prec. data) C NOPIX = MAX(NPIXF(1),NPIXF(2)) CALL STFCRE('wworkkdr',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE, + NOPIX,IMNODR,STAT) CALL STFMAP(IMNODR,F_X_MODE,1,NOPIX,IAV,DR,STAT) CALL STFCRE('wworkkdi',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE, + NOPIX,IMNODI,STAT) CALL STFMAP(IMNODI,F_X_MODE,1,NOPIX,IAV,DI,STAT) C CALL STFCRE('wworkkts',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE, + NOPIX,IMNOTS,STAT) CALL STFMAP(IMNOTS,F_X_MODE,1,NOPIX,IAV,TS,STAT) CALL STFCRE('wworkktc',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE, + NOPIX,IMNOTC,STAT) CALL STFMAP(IMNOTC,F_X_MODE,1,NOPIX,IAV,TC,STAT) C C and allocate space for swap tables CALL STFXMP(NOPIX,D_I4_FORMAT,SWAP1,STAT) SWAP2 = SWAP1 + NOPIX/2 C C now do it in a subroutine IF (FRAMEB(1:1).NE.'+') THEN !imaginary part exists... IF (TMP(2:2).EQ.'I') THEN CALL DFFTB(MADRID(PNTRA),MADRID(PNTRB), + MADRID(PNTRC),MADRID(PNTRD), + NPIXF,ACTION,MADRID(DR),MADRID(DI),NLOG,POW2, + MADRID(TS),MADRID(TC), + MADRID(CT),MADRID(DT), + MADRID(SWAP1),MADRID(SWAP2)) ELSE CALL DFFTF(MADRID(PNTRA),MADRID(PNTRB), + MADRID(PNTRC),MADRID(PNTRD), + NPIXF,ACTION,MADRID(DR),MADRID(DI),NLOG,POW2, + MADRID(TS),MADRID(TC), + MADRID(CT),MADRID(DT), + MADRID(SWAP1),MADRID(SWAP2)) ENDIF CALL STFCLO(IMNOB,STAT) ELSE !just real data IF (TMP(2:2).EQ.'I') THEN CALL DRFFTB(MADRID(PNTRA),MADRID(PNTRC), + MADRID(PNTRD),NPIXF,ACTION, + MADRID(DR),MADRID(DI),NLOG,POW2, + MADRID(TS),MADRID(TC), + MADRID(CT),MADRID(DT), + MADRID(SWAP1),MADRID(SWAP2)) ELSE CALL DRFFTF(MADRID(PNTRA),MADRID(PNTRC), + MADRID(PNTRD),NPIXF,ACTION, + MADRID(DR),MADRID(DI),NLOG,POW2, + MADRID(TS),MADRID(TC), + MADRID(CT),MADRID(DT), + MADRID(SWAP1),MADRID(SWAP2)) ENDIF ENDIF C C get rid of all the allocated temporary space CALL STFCLO(IMNODR,STAT) CALL STFCLO(IMNODI,STAT) CALL STFCLO(IMNOTS,STAT) CALL STFCLO(IMNOTC,STAT) C CALL STFCLO(IMNOA,STAT) CALL STFOPN(FRAMEA,D_OLD_FORMAT,0,F_IMA_TYPE,IMNOA,STAT) CALL DSCUPT(IMNOA,IMNOC,' ',STAT) CALL STDWRD(IMNOC,'FFT_SPECIAL',DBUF,1,6,UNI,STAT) C C if we also want power spectrum (absolute value), do it now... IF (TMP(2:2).EQ.'P') THEN CALL STKRDC('P5',1,1,60,IAV,FRAMEP,UNI,NULO,STAT) IF (EXTFLG.EQ.1) THEN WFRAME(1:) = 'exp_pow ' ELSE WFRAME(1:) = FRAMEP(1:) ENDIF CALL STIPUT(WFRAME,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISA,NPIXF,STARTA,STEPA, + IDENT,CUNIT,DR,IMNODR,STAT) CALL POWER(MADRID(PNTRC),MADRID(PNTRD),MADRID(DR),SIZE) CALL DSCUPT(IMNOA,IMNODR,' ',STAT) CALL STDWRD(IMNODR,'FFT_SPECIAL',DBUF,1,6,UNI,STAT) C IF (EXTFLG.EQ.1) THEN CALL STIPUT(FRAMEP,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISA,NPIXA,STARTA,STEPA, + IDENT,CUNIT,PNTRP,IMNOP,STAT) CALL COPWND(MADRID(DR),NPIXF,MADRID(PNTRP),NPIXA, + LOW,LOW,NPIXA) SIZE = NPIXA(1)*NPIXA(2) CALL MYMX(MADRID(PNTRP),SIZE,CUTS(1)) ELSE CALL MYMX(MADRID(DR),SIZE,CUTS(1)) ENDIF CUTS(3) = CUTS(1) CUTS(4) = CUTS(2) CALL STDWRR(IMNODR,'LHCUTS',CUTS,1,4,UNI,STAT) CALL STFCLO(IMNODR,STAT) IF (EXTFLG.EQ.1) THEN !also update extracted frame CALL DSCUPT(IMNOA,IMNOP,' ',STAT) CALL STDWRR(IMNOP,'LHCUTS',CUTS,1,4,UNI,STAT) CALL STDWRD(IMNOP,'FFT_SPECIAL',DBUF,1,6,UNI,STAT) CALL STFCLO(IMNOP,STAT) ENDIF ENDIF C C for inverse FFT, check DBUF(5,6) to see if we have to extract frames IF ((TMP(2:2).EQ.'I').AND.(DSTAT.EQ.0)) THEN NLOG(1) = NINT(DBUF(5)) NLOG(2) = NINT(DBUF(6)) IF ((NLOG(1).NE.NPIXA(1)) .OR. (NLOG(2).NE.NPIXA(2))) THEN FRAMEC(1:) = 'orig_r ' FRAMED(1:) = 'orig_i ' NPIXA(1) = NLOG(1) NPIXA(2) = NLOG(2) EXTFLG = 2 ENDIf ENDIF C C if necessary, extract result frames IF (EXTFLG.GE.1) THEN CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISA,NPIXA,STARTA,STEPA, + IDENT,CUNIT,DR,IMNODR,STAT) CALL COPWND(MADRID(PNTRC),NPIXF,MADRID(DR),NPIXA, + LOW,LOW,NPIXA) CALL STFCLO(IMNOC,STAT) IMNOC = IMNODR PNTRC = DR CALL DSCUPT(IMNOA,IMNOC,' ',STAT) CALL STDWRD(IMNOC,'FFT_SPECIAL',DBUF,1,6,UNI,STAT) C CALL STIPUT(FRAMED,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISA,NPIXA,STARTA,STEPA, + IDENT,CUNIT,PNTRB,IMNOB,STAT) CALL COPWND(MADRID(PNTRD),NPIXF,MADRID(PNTRB),NPIXA, + LOW,LOW,NPIXA) CALL STFCLO(IMNOD,STAT) IMNOD = IMNOB PNTRD = PNTRB SIZE = NPIXA(1)*NPIXA(2) IF (EXTFLG.EQ.1) THEN TMP(3:3) = 'Y' !tell fft.prg about the extension ELSE TMP(3:3) = 'Z' !tell fft.prg about the extraction ENDIF CALL STKWRC('TMPBUF',1,TMP,1,3,UNI,STAT) ENDIF C C get dynamic range CALL MYMX(MADRID(PNTRC),SIZE,CUTS(1)) CUTS(3) = CUTS(1) CUTS(4) = CUTS(2) CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT) CALL MYMX(MADRID(PNTRD),SIZE,CUTS(1)) CUTS(3) = CUTS(1) CUTS(4) = CUTS(2) CALL STDWRR(IMNOD,'LHCUTS',CUTS,1,4,UNI,STAT) C C That's it folks... CALL STSEPI C C here if real + imaginary input frames don't match 9000 CALL STETER(2,'input frames do not have same dimensions...') END SUBROUTINE DFFTF(A,B,C,D,NPIX,ACT,DREAL,DIMAG,NLOG,POW2, + TSIN,TCOS,CT,DT,SWAPI,SWAPO) C IMPLICIT NONE C INTEGER INDX,N,NDIM,NX,NY,MM INTEGER XFROM,NSIZE(2),NWORK INTEGER KMM,MPOW(21) INTEGER NPIX(2),NLOG(*),POW2(*) INTEGER SWAPI(*),SWAPO(*),SWALIM C REAL A(*),B(*),C(*),D(*),CT(*),DT(*) C DOUBLE PRECISION DREAL(*),DIMAG(*),TSIN(*),TCOS(*) DOUBLE PRECISION DFACT DOUBLE PRECISION DX,DY C CHARACTER*1 ACT C C init NDIM = NPIX(1) DFACT = 1.D0 / DBLE(FLOAT(NDIM)) C C compute array MPOW: 2**(NLOG(1)-1), ..., 2**1, 2**0 MM = NLOG(1) DO 200, N=1,NLOG(1) MPOW(N) = POW2(MM) MM = MM - 1 200 CONTINUE C C generate sine + cosine look-up table in x-dimension CALL DSCTB(TSIN,TCOS,NDIM,NLOG(1),SWAPI,SWAPO,SWALIM) C C here we handle 1-dim FFT C IF (NPIX(2).EQ.1) THEN IF (ACT.NE.'F') THEN DO 400, NX=1,NDIM DREAL(NX) = DBLE(A(NX)) !copy A to DREAL DIMAG(NX) = DBLE(B(NX)) !copy B to DIMAG 400 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM) !do the FFT ELSE C C here for frequency business (1-dim case) C MM = 0 DO 1000, NX=1,NDIM DX = DBLE(A(NX)) DY = DBLE(B(NX)) IF (MM.EQ.0) THEN DREAL(NX) = DX !copy A*(-1)**(NX+NY) to DREAL DIMAG(NX) = DY !copy B*(-1)**(NX+NY) to DIMAG MM = 1 ELSE DREAL(NX) = -DX !copy A*(-1)**(NX+NY) to DREAL DIMAG(NX) = -DY !copy B*(-1)**(NX+NY) to DIMAG MM = 0 ENDIF 1000 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM) !do the FFT ENDIF RETURN ENDIF C C C -------------- here we handle 2-dim FFT ------------------- C NSIZE(1) = 128 NSIZE(2) = 256 C C branch, if frequency scaling is wanted XFROM = 1 IF (ACT.EQ.'F') GOTO 6000 NWORK = NDIM - 1 C C do FFT on each line DO 5800, NY=1,NPIX(2) N = 1 DO 5400, NX = XFROM,XFROM+NWORK DREAL(N) = DBLE(A(NX)) !copy A to DREAL DIMAG(N) = DBLE(B(NX)) !copy B to DIMAG N = N + 1 5400 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C(XFROM),D(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT C XFROM = XFROM + NDIM 5800 CONTINUE C GOTO 7000 C C here for frequency business (2-dim case) 6000 KMM = 0 DO 6800, NY=1,NPIX(2) N = 1 MM = KMM NWORK = XFROM - 1 C DO 6400, NX = 1,NDIM INDX = NWORK + NX DX = DBLE(A(INDX)) DY = DBLE(B(INDX)) IF (MM.EQ.0) THEN DREAL(N) = DX !copy A*(-1)**(NX+NY) to DREAL DIMAG(N) = DY !copy B*(-1)**(NX+NY) to DIMAG MM = 1 ELSE DREAL(N) = -DX !copy A*(-1)**(NX+NY) to DREAL DIMAG(N) = -DY !copy B*(-1)**(NX+NY) to DIMAG MM = 0 ENDIF N = N + 1 6400 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C(XFROM),D(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT C XFROM = XFROM + NDIM KMM = 1 - KMM !alternate between 0 and 1 6800 CONTINUE C C we now do another FFT along the columns 7000 CALL LINCOL(C,NPIX,NSIZE,CT) !change lines -> columns CALL LINCOL(D,NPIX,NSIZE,DT) NDIM = NPIX(2) DFACT = DFACT / DBLE(FLOAT(NDIM)) C C compute array MPOW in y-dim: 2**(NLOG(2)-1), ..., 2**1, 2**0 IF (NLOG(2).NE.NLOG(1)) THEN MM = NLOG(2) DO 7400, N=1,NLOG(2) MPOW(N) = POW2(MM) MM = MM - 1 7400 CONTINUE C C and generate sine + cosine look-up table in y-dimension CALL DSCTB(TSIN,TCOS,NDIM,NLOG(2),SWAPI,SWAPO,SWALIM) ENDIF C XFROM = 1 NWORK = NDIM - 1 C DO 8000, NY=1,NPIX(1) N = 1 DO 7700, NX=XFROM,XFROM+NWORK DREAL(N) = DBLE(CT(NX)) DIMAG(N) = DBLE(DT(NX)) N = N + 1 7700 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW, + TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT C XFROM = XFROM + NDIM 8000 CONTINUE C C finally transpose back again CALL TRAPO(CT,DT,NPIX,NSIZE,C,D) C RETURN END SUBROUTINE DFFTB(A,B,C,D,NPIX,ACT,DREAL,DIMAG,NLOG,POW2, + TSIN,TCOS,CT,DT,SWAPI,SWAPO) C IMPLICIT NONE C INTEGER N,NDIM,NX,NY,MM INTEGER XFROM,NSIZE(2),NWORK INTEGER NPIX(2),NLOG(*),POW2(*) INTEGER KMM,MPOW(21) INTEGER SWAPI(*),SWAPO(*),SWALIM C REAL A(*),B(*),C(*),D(*),CT(*),DT(*) C DOUBLE PRECISION DREAL(*),DIMAG(*),TSIN(*),TCOS(*) DOUBLE PRECISION DFACT C CHARACTER*1 ACT C C init NDIM = NPIX(1) DFACT = 1.D0 / DBLE(FLOAT(NDIM)) C C compute array MPOW: 2**(NLOG(1)-1), ..., 2**1, 2**0 MM = NLOG(1) DO 200, N=1,NLOG(1) MPOW(N) = POW2(MM) MM = MM - 1 200 CONTINUE C C generate sine + cosine look-up table in x-dimension CALL DSCTB(TSIN,TCOS,NDIM,NLOG(1),SWAPI,SWAPO,SWALIM) C C do 1-dim inverse FFT IF (NPIX(2).EQ.1) THEN IF (ACT.NE.'F') THEN DO 400, NX=1,NDIM DREAL(NX) = DBLE(A(NX)) !copy A to DREAL DIMAG(NX) = DBLE(-B(NX)) !copy -B to DIMAG 400 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM) !do the FFT ELSE C C frequency business... DO 600, NX=1,NDIM DREAL(NX) = DBLE(A(NX)) !copy A to DREAL DIMAG(NX) = DBLE(-B(NX)) !copy -B to DIMAG 600 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM) !do the FFT C MM = 0 DO 1000, NX=1,NDIM IF (MM.EQ.0) THEN MM = 1 ELSE C(NX) = - C(NX) D(NX) = - D(NX) MM = 0 ENDIF 1000 CONTINUE ENDIF C RETURN ENDIF C C C ------------------------- do 2-dim inverse FFT ------------------------- C NSIZE(1) = 128 NSIZE(2) = 256 C C handle each line XFROM = 1 NWORK = NDIM - 1 C DO 5800, NY=1,NPIX(2) N = 1 DO 5400, NX=XFROM,XFROM+NWORK DREAL(N) = DBLE(A(NX)) !copy A to DREAL DIMAG(N) = DBLE(-B(NX)) !copy -B to DIMAG N = N + 1 5400 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C(XFROM),D(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT XFROM = XFROM + NDIM 5800 CONTINUE C C handle each column CALL LINCOL(C,NPIX,NSIZE,CT) CALL LINCOL(D,NPIX,NSIZE,DT) NDIM = NPIX(2) CC DFACT = DFACT * DBLE(FLOAT(NDIM)) DFACT = 1.D0 C C compute array MPOW in y-dim: 2**(NLOG(2)-1), ..., 2**1, 2**0 IF (NLOG(2).NE.NLOG(1)) THEN MM = NLOG(2) DO 6000, N=1,NLOG(2) MPOW(N) = POW2(MM) MM = MM - 1 6000 CONTINUE C C and generate sine + cosine look-up table in y-dimension CALL DSCTB(TSIN,TCOS,NDIM,NLOG(2),SWAPI,SWAPO,SWALIM) ENDIF C XFROM = 1 NWORK = NDIM - 1 C IF (ACT.NE.'F') THEN DO 6800, NY=1,NPIX(1) N = 1 DO 6400, NX=XFROM,XFROM+NWORK DREAL(N) = DBLE(CT(NX)) DIMAG(N) = DBLE(DT(NX)) N = N + 1 6400 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW, + TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT XFROM = XFROM + NDIM 6800 CONTINUE C ELSE C C frequency business... KMM = 0 DO 7800, NY=1,NPIX(1) N = 1 DO 7400, NX=XFROM,XFROM+NWORK DREAL(N) = DBLE(CT(NX)) DIMAG(N) = DBLE(DT(NX)) N = N + 1 7400 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW, + TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT C MM = KMM DO 7600, NX=XFROM,XFROM+NWORK IF (MM.EQ.0) THEN MM = 1 ELSE CT(NX) = - CT(NX) DT(NX) = - DT(NX) MM = 0 ENDIF 7600 CONTINUE C XFROM = XFROM + NDIM KMM = 1 - KMM !alternate between 0 and 1 7800 CONTINUE C ENDIF C C finally transpose back again CALL TRAPO(CT,DT,NPIX,NSIZE,C,D) C RETURN END SUBROUTINE DRFFTF(A,C,D,NPIX,ACT,DREAL,DIMAG,NLOG,POW2, + TSIN,TCOS,CT,DT,SWAPI,SWAPO) C IMPLICIT NONE C INTEGER INDX,N,NDIM,NX,NY,MM INTEGER XFROM,NSIZE(2),NWORK INTEGER NPIX(2),NLOG(*),POW2(*) INTEGER KMM,MPOW(21) INTEGER SWAPI(*),SWAPO(*),SWALIM C REAL A(*),C(*),D(*),CT(*),DT(*) C DOUBLE PRECISION DREAL(*),DIMAG(*),TSIN(*),TCOS(*) DOUBLE PRECISION DFACT DOUBLE PRECISION DX C CHARACTER*1 ACT C C init NDIM = NPIX(1) DFACT = 1.D0 / DBLE(FLOAT(NDIM)) C C compute array MPOW: 2**(NLOG(1)-1), ..., 2**1, 2**0 MM = NLOG(1) DO 200, N=1,NLOG(1) MPOW(N) = POW2(MM) MM = MM - 1 200 CONTINUE C C generate sine + cosine look-up table in x-dimension CALL DSCTB(TSIN,TCOS,NDIM,NLOG(1),SWAPI,SWAPO,SWALIM) C C here we handle 1-dim FFT C IF (NPIX(2).EQ.1) THEN IF (ACT.NE.'F') THEN C DO 400, NX=1,NDIM DREAL(NX) = DBLE(A(NX)) !copy A to DREAL DIMAG(NX) = 0.D0 400 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM) !do the FFT ELSE C C here for frequency business MM = 0 DO 1000, NX=1,NDIM DX = DBLE(A(NX)) IF (MM.EQ.0) THEN DREAL(NX) = DX !copy A*(-1)**(NX+NY) to DREAL MM = 1 ELSE DREAL(NX) = -DX !copy A*(-1)**(NX+NY) to DREAL MM = 0 ENDIF DIMAG(NX) = 0.D0 1000 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM) !do the FFT ENDIF C RETURN ENDIF C C C ---------------------- here we handle 2-dim FFT --------------------- C NSIZE(1) = 128 NSIZE(2) = 256 C C branch, if frequency scaling is wanted XFROM = 1 IF (ACT.EQ.'F') GOTO 6000 NWORK = NDIM - 1 C C do FFT on each line DO 5800, NY=1,NPIX(2) N = 1 DO 5400, NX = XFROM,XFROM+NWORK DREAL(N) = DBLE(A(NX)) !copy A to DREAL DIMAG(N) = 0.D0 N = N + 1 5400 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C(XFROM),D(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT XFROM = XFROM + NDIM 5800 CONTINUE C GOTO 7000 C C here for frequency business 6000 KMM = 0 DO 6800, NY=1,NPIX(2) N = 1 MM = KMM NWORK = XFROM - 1 C DO 6400, NX = 1,NDIM INDX = NWORK + NX DX = DBLE(A(INDX)) IF (MM.EQ.0) THEN DREAL(N) = DX !copy A*(-1)**(NX+NY) to DREAL MM = 1 ELSE DREAL(N) = -DX !copy A*(-1)**(NX+NY) to DREAL MM = 0 ENDIF DIMAG(N) = 0.D0 N = N + 1 6400 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C(XFROM),D(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT C KMM = 1 - KMM XFROM = XFROM + NDIM 6800 CONTINUE C C we now do another FFT along the columns 7000 CALL LINCOL(C,NPIX,NSIZE,CT) !change lines -> columns CALL LINCOL(D,NPIX,NSIZE,DT) NDIM = NPIX(2) DFACT = DFACT / DBLE(FLOAT(NDIM)) C C compute array MPOW in y-dim: 2**(NLOG(2)-1), ..., 2**1, 2**0 IF (NLOG(2).NE.NLOG(1)) THEN MM = NLOG(2) DO 7400, N=1,NLOG(2) MPOW(N) = POW2(MM) MM = MM - 1 7400 CONTINUE C C and generate sine + cosine look-up table in y-dimension CALL DSCTB(TSIN,TCOS,NDIM,NLOG(2),SWAPI,SWAPO,SWALIM) ENDIF C XFROM = 1 NWORK = NDIM - 1 C DO 8000, NY=1,NPIX(1) N = 1 DO 7700, NX=XFROM,XFROM+NWORK DREAL(N) = DBLE(CT(NX)) DIMAG(N) = DBLE(DT(NX)) N = N + 1 7700 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW, + TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT XFROM = XFROM + NDIM 8000 CONTINUE C C finally transpose back again CALL TRAPO(CT,DT,NPIX,NSIZE,C,D) C RETURN END SUBROUTINE DRFFTB(A,C,D,NPIX,ACT,DREAL,DIMAG,NLOG,POW2, + TSIN,TCOS,CT,DT,SWAPI,SWAPO) C IMPLICIT NONE C INTEGER N,NDIM,NX,NY,MM INTEGER XFROM,NSIZE(2),NWORK INTEGER NPIX(2),NLOG(*),POW2(*) INTEGER KMM,MPOW(21) INTEGER SWAPI(*),SWAPO(*),SWALIM C REAL A(*),C(*),D(*),CT(*),DT(*) C DOUBLE PRECISION DREAL(*),DIMAG(*),TSIN(*),TCOS(*) DOUBLE PRECISION DFACT C CHARACTER*1 ACT C C init NDIM = NPIX(1) DFACT = 1.D0 / DBLE(FLOAT(NDIM)) C C compute array MPOW: 2**(NLOG(1)-1), ..., 2**1, 2**0 MM = NLOG(1) DO 200, N=1,NLOG(1) MPOW(N) = POW2(MM) MM = MM - 1 200 CONTINUE C C generate sine + cosine look-up table in x-dimension CALL DSCTB(TSIN,TCOS,NDIM,NLOG(1),SWAPI,SWAPO,SWALIM) C C do 1-dim inverse FFT IF (NPIX(2).EQ.1) THEN IF (ACT.NE.'F') THEN DO 400, NX=1,NDIM DREAL(NX) = DBLE(A(NX)) !copy A to DREAL DIMAG(NX) = 0.D0 400 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM) !do the FFT ELSE C C frequency business... DO 600, NX=1,NDIM DREAL(NX) = DBLE(A(NX)) !copy A to DREAL DIMAG(NX) = 0.D0 600 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM) !do the FFT C MM = 0 DO 1000, NX=1,NDIM IF (MM.EQ.0) THEN MM = 1 ELSE C(NX) = - C(NX) D(NX) = - D(NX) MM = 0 ENDIF 1000 CONTINUE C ENDIF RETURN ENDIF C C C ------------------------- do 2-dim inverse FFT ------------------------- C NSIZE(1) = 128 NSIZE(2) = 256 C C handle each line XFROM = 1 NWORK = NDIM - 1 C DO 5800, NY=1,NPIX(2) N = 1 DO 5400, NX=XFROM,XFROM+NWORK DREAL(N) = DBLE(A(NX)) !copy A to DREAL DIMAG(N) = 0.D0 N = N + 1 5400 CONTINUE C CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW, + TSIN,TCOS,DFACT,C(XFROM),D(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT C XFROM = XFROM + NDIM 5800 CONTINUE C C handle each column CALL LINCOL(C,NPIX,NSIZE,CT) CALL LINCOL(D,NPIX,NSIZE,DT) NDIM = NPIX(2) DFACT = DFACT * DBLE(FLOAT(NDIM)) C C compute array MPOW in y-dim: 2**(NLOG(2)-1), ..., 2**1, 2**0 IF (NLOG(2).NE.NLOG(1)) THEN MM = NLOG(2) DO 6000, N=1,NLOG(2) MPOW(N) = POW2(MM) MM = MM - 1 6000 CONTINUE C C and generate sine + cosine look-up table in y-dimension CALL DSCTB(TSIN,TCOS,NDIM,NLOG(2),SWAPI,SWAPO,SWALIM) ENDIF C XFROM = 1 NWORK = NDIM - 1 C IF (ACT.NE.'F') THEN DO 6800, NY=1,NPIX(1) N = 1 DO 6400, NX=XFROM,XFROM+NWORK DREAL(N) = DBLE(CT(NX)) DIMAG(N) = DBLE(DT(NX)) N = N + 1 6400 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW, + TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT C XFROM = XFROM + NDIM 6800 CONTINUE C ELSE C C frequency business... KMM = 0 DO 7800, NY=1,NPIX(1) N = 1 DO 7400, NX=XFROM,XFROM+NWORK DREAL(N) = DBLE(CT(NX)) DIMAG(N) = DBLE(DT(NX)) N = N + 1 7400 CONTINUE C CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW, + TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM), + SWAPI,SWAPO,SWALIM) !do the FFT C MM = KMM DO 7600, NX=XFROM,XFROM+NWORK IF (MM.EQ.0) THEN MM = 1 ELSE CT(NX) = - CT(NX) DT(NX) = - DT(NX) MM = 0 ENDIF 7600 CONTINUE C XFROM = XFROM + NDIM KMM = 1 - KMM 7800 CONTINUE C ENDIF C C finally transpose back again CALL TRAPO(CT,DT,NPIX,NSIZE,C,D) C RETURN END SUBROUTINE DFFT(IFACT,XREAL,XIMAG,N,NU,MPOW,TSIN,TCOS,DFACT, + XR,XI,SWAPI,SWAPO,SWLIM) C C FAST FOURIER TRANSFORM ALGORITHM C K. Banse version 2.10 850625 C 2.20 860217 C 2.40 870109 C 2.50 921022 C C IFACT - factor flag: 1 - yes, multiply by DFACT C 0 - no, do not multiply C XREAL - double precision real part of input array C XIMAG - double precision imaginary part of input array C N - dimension of input array C NU - LOG(N) i.e.: N = 2**NU C MPOW - array with powers of 2 from 1 - N C TSIN - array precomputed sine values C TCOS - array precomputed cosine values C DFACT - double precision constant for multiplication (IFACT = 1) C XR - real part of XREAL C XI - imag part of XIMAG C IMPLICIT NONE C INTEGER IFACT,N,NU INTEGER N1,MPL INTEGER I,K,K1,K1N2,L INTEGER N2,LP INTEGER MPOW(*) INTEGER SWAPI(*),SWAPO(*),SWLIM C DOUBLE PRECISION XREAL(*),XIMAG(*),TSIN(*),TCOS(*),DFACT DOUBLE PRECISION TREAL,TIMAG,TCC,TSS C REAL XR(*),XI(*) REAL TR C C initialize N2 = MPOW(1) ! the same as: N2 = N/2 N1 = N - 1 C C handle vertical arrays DO 500, L=1,NU K = 0 MPL = MPOW(L) C 50 DO 200, I=1,N2 LP = K/MPL + 1 TCC = TCOS(LP) TSS = TSIN(LP) K = K + 1 K1N2 = K + N2 TREAL = XREAL(K1N2)*TCC + XIMAG(K1N2)*TSS TIMAG = XIMAG(K1N2)*TCC - XREAL(K1N2)*TSS XREAL(K1N2) = XREAL(K) - TREAL XIMAG(K1N2) = XIMAG(K) - TIMAG XREAL(K) = XREAL(K) + TREAL XIMAG(K) = XIMAG(K) + TIMAG 200 CONTINUE C K = K + N2 IF (K.LT.N1) GOTO 50 C N2 = N2/2 !set up variables for next array 500 CONTINUE C C convert from double precision to single + multiply (if IFACT = 1) IF (IFACT.EQ.1) THEN DO 1000, K=1,N XR(K) = SNGL(XREAL(K) * DFACT) XI(K) = SNGL(XIMAG(K) * DFACT) 1000 CONTINUE ELSE DO 1200, K=1,N XR(K) = SNGL(XREAL(K)) XI(K) = SNGL(XIMAG(K)) 1200 CONTINUE ENDIF C C unscramble final vector DO 2000, K=1,SWLIM K1 = SWAPI(K) I = SWAPO(K) TR = XR(K1) XR(K1) = XR(I) XR(I) = TR TR = XI(K1) XI(K1) = XI(I) XI(I) = TR 2000 CONTINUE C C that's it folks... RETURN END SUBROUTINE DSCTB(TSIN,TCOS,NDIM,NU,TIN,TOUT,SWLIM) C C subroutine to generate sine and cosine look-up tables C and also to generate swap tables C IMPLICIT NONE C INTEGER IBITR,NDIM,NU,K,M,IBR INTEGER TIN(NDIM),TOUT(NDIM),SWLIM C DOUBLE PRECISION TSIN(*),TCOS(*),ARG,ARG1 C ARG1 = 6.2831852D0 / DBLE(FLOAT(NDIM)) SWLIM = 0 C DO 1000, K=1,NDIM M = K - 1 IBR = IBITR(M,NU) ARG = ARG1 * IBR TSIN(K) = SIN(ARG) TCOS(K) = COS(ARG) C IF (IBR.GT.M) THEN !fill swap tables already SWLIM = SWLIM + 1 TIN(SWLIM) = K TOUT(SWLIM) = IBR + 1 ENDIF 1000 CONTINUE C RETURN END SUBROUTINE TRAPO(A,B,NPIX,NSIZE,C,D) C IMPLICIT NONE C INTEGER NPIX(2),NSIZE(2) INTEGER NPAX(2) C REAL A(*),B(*),C(*),D(*) C NPAX(1) = NPIX(2) NPAX(2) = NPIX(1) CALL LINCOL(A,NPAX,NSIZE,C) CALL LINCOL(B,NPAX,NSIZE,D) C RETURN END FUNCTION IBITR(J,NU) C IMPLICIT NONE C INTEGER IBITR,J,NU INTEGER I,J1,J2 C J1 = J IBITR = 0 DO 1000 I=1,NU J2 = J1/2 IBITR = IBITR*2 + (J1-2*J2) J1 = J2 1000 CONTINUE C RETURN END SUBROUTINE POWER(A,B,C,NDIM) C IMPLICIT NONE C INTEGER N,NDIM C REAL A(NDIM),B(NDIM),C(NDIM) C DO 1000 N=1,NDIM C(N) = SQRT( ( A(N)*A(N) ) + ( B(N)*B(N) ) ) 1000 CONTINUE C RETURN END