C TFRINV SOURCE OF166741 25/02/20 21:17:46 12165 SUBROUTINE TFRINV C C======================================================================= C = CALCUL DE LA TRANSFORMEE DE FOURIER INVERSE D'UN SIGNAL C = C = SYNTAXE : FFTR = TFRI EVOL ; C = C = EVOL : OBJET DE TYPE EVOLUTIO CONTENANT L'OBJET EVOLUTION C = = C = CREATION : 26/04/88 = C = F. ROULLIER C======================================================================= C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) CHARACTER *72 TI CHARACTER *4 ITYP1 c RLTOCX : fonction de conversion real*8 -> complex*16 COMPLEX*16 RLTOCX LOGICAL INV -INC PPARAM -INC CCOPTIO -INC SMEVOLL -INC SMLREEL -INC CCREEL C c SEGMENT COMPL c COMPLEX*16 XX(NPT),YY(NPT),W(NEXP) c REAL*8 AR(NPOIN),AI(NPOIN) c REAL*8 DM(NPTD),DP(NPTD) c ENDSEGMENT c segment pour appel a fftpack5.1 segment wfft51 real*8 work(lenwrk) real*8 wsave(lensav) real*8 XX51(NCOU) endsegment C C==== LECTURES ========================================================= c C LECTURE DE L'OBJET EVOLUTION COMPLEXE A TRANSFORMER C CALL LIROBJ('EVOLUTIO',IEV1,1,IRETOU) IF(IERR.NE.0) RETURN CALL ACTOBJ('EVOLUTIO',IEV1,1) C MEVOLL=IEV1 IF (ITYEVO.NE.'COMPLEXE') THEN MOTERR(1:8)='EVOLUTION' CALL ERREUR (302) RETURN ENDIF KEVOL1=IEVOLL(1) KEVOL2=IEVOLL(2) ITYP1=KEVOL1.NUMEVY ICOU1=KEVOL1.NUMEVX MLREE3=KEVOL1.IPROGX MLREE1=KEVOL1.IPROGY MLREE2=KEVOL2.IPROGY C c on recupere le pas en frequence, la dimension et les ordonnees L1=MLREE3.PROG(/1) DF=MLREE3.PROG(2)-MLREE3.PROG(1) C c on verifie le nombre de points NPT=(L1-1)*2 NP=NPT NN=0 DO 1 I=1,1000 NP=NP/2 NN=NN+1 IF (NP.LE.1) GO TO 2 1 CONTINUE 2 NP=2**NN IF (NP.NE.NPT) THEN MOTERR(1:8)='EVOLUTIO' CALL ERREUR (1083) RETURN ENDIF c on verifie -que la 1ere frequence soit 0, IF(MLREE3.PROG(1).NE.0.D0) THEN c La premiere frequence doit etre egale a 0. REAERR(1)=0.d0 CALL ERREUR (1084) RETURN ENDIF C -que la partie imag de la 1ere soit 0 ... IF(ITYP1.EQ.'PREE') THEN DR=MLREE1.PROG(1) DI=MLREE2.PROG(1) DM=SQRT(DR**2+DI**2) ELSE DM = MLREE1.PROG(1) DP = MLREE2.PROG(1) * XPI/180. DI = DM*SIN(DP) ENDIF IF(ABS(DI).GT.(XZPREC*DM)) THEN c La 1 eme valeur imaginaire doit etre egale a 0. INTERR(1)=1 REAERR(1)=0.d0 CALL ERREUR (1086) RETURN ENDIF C c NEXP=NPT/2 c NPOIN=L1 c IF (ITYP1.NE.'PREE') THEN c NPTD=L1 c ELSE c NPTD=0 c ENDIF c SEGINI COMPL c write(*,*) 'TFR en entree de dim:',L1 c write(*,*) 'nombre de point en sortie',NPT,NEXP c write(*,*) 'dim des tableaux de travail',DM(/1),AR(/1),W(/1) C creation du tableau de travail pour FFTPACK5.1 NCOU=NPT lenwrk = 2 * NPT lensav = NPT + int(log ( real (NPT) ) / log (2.0) ) + 4 segini,wfft51 C C==== CHARGEMENT DES TABLEAUX DE TRAVAIL =============================== C c - cas Reel,Imag : IF (ITYP1.EQ.'PREE') THEN c on a forcement le terme constant (f=0Hz) XX51(1) = MLREE1.PROG(1) c la suite par paire ITY=1 DO 10 I=2,L1-1 ITY=ITY+1 XX51(2*I-2)=2. * (-1)**(2*I) * MLREE1.PROG(ITY) XX51(2*I-1)=2. * (-1)**(2*I+1) * MLREE2.PROG(ITY) 10 CONTINUE c dernier terme : on suppose que la TF n'a pas ete tronquee XX51(NPT) = MLREE1.PROG(L1) c - cas Module, Phase : ELSE RAP=XPI/180. c on a forcement le terme constant (f=0Hz) c XX51(1) = MLREE1.PROG(1) DM = MLREE1.PROG(1) DP = MLREE2.PROG(1) * RAP DR = DM*COS(DP) DI = DM*SIN(DP) c write(*,*) '1: DM,DP,DR,DI=',DM,DP,DR,DI XX51(1) = DR c la suite par paire ITY=1 DO 11 I=2,L1-1 ITY=ITY+1 c MODULE , PHASE -> REEL, IMAG DM = MLREE1.PROG(ITY) DP = MLREE2.PROG(ITY) * RAP DR = DM*COS(DP) DI = DM*SIN(DP) c write(*,*) ITY,I,': DM,DP,DR,DI=',DM,DP,DR,DI XX51(2*I-2)=2.0 * (-1)**(2*I) * DR XX51(2*I-1)=2.0 * (-1)**(2*I+1) * DI 11 CONTINUE c dernier terme : on suppose que la TF n'a pas ete tronquee c XX51(NPT) = 0.25d0 * MLREE1.PROG(L1) DM = MLREE1.PROG(L1) DP = MLREE2.PROG(L1) * RAP DR = DM*COS(DP) DI = DM*SIN(DP) XX51(NPT) = DR c c si la TF est fournie sous la forme (module,phase), conversion c DO 8 I=1,L1 c DM(I)=MLREE1.PROG(I) c DP(I)=MLREE2.PROG(I) c 8 CONTINUE c CALL CONVCP(AR,AI,DM,DP,L1,-1) c ELSE c c sinon on recupere directement la TF (Re,Im) c DO 9 I=1,L1 c AR(I)=MLREE1.PROG(I) c AI(I)=MLREE2.PROG(I) c 9 CONTINUE ENDIF C==== CALCUL DE LA FFT inverse ========================================= c +++ via FFTPACK5.1 +++ c Initialize the WSAVE array. call rfft1i (NPT,wsave,lensav,ier) IF (IERR.ne.0) RETURN c c Compute inverse FFT of coefficients. Should get back the c original data. inc = 1 lenr = NPT call rfft1b(NPT,inc,XX51,lenr,wsave,lensav,work,lenwrk,ier) IF (IERR.ne.0) RETURN 19 IF(IIMPI.EQ.1) WRITE(IOIMP,*)' FFT INVERSE CALCULEE ' C C==== CREATION ET CALCUL DES LISTES DE LA FFT ========================== C JG=NPT SEGINI MLREE1 SEGINI MLREE2 FMAX=DF*REAL(NPT) DTEMP=1.D0/FMAX C DO 20 I=1,NPT TEMP=REAL(I-1)*DTEMP MLREE1.PROG(I)=TEMP MLREE2.PROG(I)=DF*XX51(I) 20 CONTINUE SEGSUP wfft51 C C CREATION DE L'OBJET EVOLUTION REEL FFT INVERSE C N=1 SEGINI MEVOLL IPVO=MEVOLL IEVTEX=TITREE ITYEVO='REEL ' C SEGINI KEVOLL c KEVTEX=TITREE KEVTEX='TF^{-1}' TYPX='LISTREEL' TYPY='LISTREEL' IEVOLL(1)=KEVOLL IPROGX=MLREE1 NOMEVX='TEMPS SEC' IPROGY=MLREE2 NOMEVY='TF INVERSE ' NUMEVX=ICOU1 NUMEVY='REEL' C CALL ACTOBJ('EVOLUTIO',IPVO,1) CALL ECROBJ('EVOLUTIO',IPVO) 666 CONTINUE END