C DSPR SOURCE BP208322 16/11/18 21:16:32 9177 SUBROUTINE DSPR IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-V) IMPLICIT COMPLEX*16 (W-Z) CHARACTER *72 TI CHARACTER*12 MOTX,MOTY CHARACTER*4 MCLE(2) LOGICAL INV C C======================================================================= C = CALCUL DE LA DENSITE SPECTRALE DE PUISSANCE D'UN SIGNAL = C = = C = SYNTAXE : = C = = C = SPEC = DSPR EXP2 EVO1 FMIN V1 FMAX V2 (COUL) ; = C = = C = = C = EXP2 : EXPOSANT DANS NPOINT=2^EXP2, NPOINT ETANT = C = LE NOMBRE DE REELS DANS LISTREEL = C = EVO1 : OBJET DE TYPE EVOLUTIO CONTENANT LE SIGNAL A TRAITER= C = ( UNE COURBE SEULEMENT ) = C = FMIN : MOT-CLE = C = V1 : FREQUENCE MINIE A VISUALISER = C = FMAX : MOT-CLE = C = V2 : FREQUENCE MAXIE A VISUALISER = C = COUL : COULEUR ATTRIBUEE A L'OBJET CREE (FACULTATIF) = C = = C = = C = CREATION : 01/04/87 = C = PROGRAMMEUR : BEAUFILS = C = MODIFICATIO : PEG WARNING LE IRET2.... 22/2/90 = C = MODIFICATIO : PEG FREQUENCES CORRECTES 1/6/90 = C = MODIFICATIO : PEG PERIODOGRAM CORRECT 1/6/90 = C======================================================================= C -INC CCGEOME -INC PPARAM -INC CCOPTIO -INC SMEVOLL -INC SMLREEL C SEGMENT MTRAV1 IMPLIED XX(NCOU),YY(NCOU) ENDSEGMENT C SEGMENT MTRAV2 IMPLIED W(NEXP) ENDSEGMENT C DATA MOTX/'FREQ HZ'/ DATA MOTY/'DSP UN**2/HZ'/ DATA MCLE(1),MCLE(2)/'FMIN','FMAX'/ C C LECTURE DE EXP2 C CALL LIRENT(N2,1,IRET1) IF(IRET1.EQ.0)GOTO 666 C C LECTURE DE L'OBJET EVOLUTIO CONTENANT LE SIGNAL C CALL LIROBJ('EVOLUTIO',IPSIG,1,IRET2) C PEG IF(IRET1.EQ.0)GOTO 666 IF(IRET2.EQ.0)GOTO 666 C C LECTURE DE LA FREQUENCE MINIE C CALL LIRMOT(MCLE(1),1,IRET,0) IF(IRET.EQ.1) THEN CALL LIRREE(FRMI,0,IRET1) IF(IRET1.EQ.0) THEN MOTERR(1:4)=MCLE(1) CALL ERREUR(166) GOTO 666 ENDIF IF(FRMI.LT.(-1.D-20)) FRMI=-FRMI IF(ABS(FRMI).LT.1.D-20) FRMI=-1.D0 ELSE FRMI=-1.D0 ENDIF C C LECTURE DE LA FREQUENCE MAXIE C CALL LIRMOT(MCLE(2),1,IRET,0) IF(IRET.EQ.1) THEN CALL LIRREE(FRMA,0,IRET1) IF(IRET1.EQ.0) THEN MOTERR(1:4)=MCLE(2) CALL ERREUR(166) GOTO 666 ENDIF IF(FRMA.LT.(-1.D-20)) FRMA=-FRMA IF(ABS(FRMA).LT.1.D-20) FRMA=-1.D0 ELSE FRMA=-1.D0 ENDIF C C LECTURE DE LA COULEUR C CALL LIRMOT(NCOUL,NBCOUL,ICOUL,0) IF(ICOUL.EQ.0) ICOUL=IDCOUL+1 ICOUL=ICOUL-1 C IF(IERR.NE.0) GOTO 666 C MEVOL1=IPSIG SEGACT MEVOL1 KEVOL1=MEVOL1.IEVOLL(1) SEGACT KEVOL1 MLREE1=KEVOL1.IPROGX MLREE2=KEVOL1.IPROGY SEGACT MLREE1,MLREE2 C L1=MLREE1.PROG(/1) DT=MLREE1.PROG(2)-MLREE1.PROG(1) SEGDES MLREE1 C NPOINT=2**N2 IF(NPOINT.GT.L1) THEN IF(IIMPI.EQ.1) WRITE(IOIMP,1000) L1,N2,NPOINT 1000 FORMAT(1H ,'LE NOMBRE DE POINTS DANS LISTEMPS ',I6, ' EST ', & 'INFERIEURE @ 2**',I5, & /' ON PRENDRA UNE LONGUEUR DE LISTEMPS DE ',I6,' MOTS ', & /' ET ON COMPLETERA PAR DES ZEROS') ELSE IF(NPOINT.LT.L1) THEN IF(IIMPI.EQ.1) WRITE(IOIMP,1010) N2 1010 FORMAT(1H ,'ON TRONQUE LE SIGNAL A 2**',I5,' MOTS',/) ELSE IF(IIMPI.EQ.1) WRITE(IOIMP,1030)N2,NPOINT 1030 FORMAT(1H ,'LA LONGUEUR DU LISTEMP EST EGALE A 2**',I5, & ' = ',I6,/) ENDIF ENDIF C NCOU=NPOINT SEGINI MTRAV1 C NEXP=NPOINT/2 SEGINI MTRAV2 C C CHARGEMENT DES TABLEAUX DE TRAVAIL C IND1=L1 IF(NPOINT.LT.L1)IND1 = NPOINT DO 10 I=1,IND1 XX(I)=MLREE2.PROG(I) C IF(IIMPI.EQ.1)WRITE(IOIMP,*) ' XX(',I,') = ',XX(I) 10 CONTINUE IF(NPOINT.GT.L1) THEN L2=L1+1 DO 11 I=L2,NPOINT XX(I)=0.D0 11 CONTINUE ENDIF C C CALCUL DE LA FFT C IF(IIMPI.EQ.1) WRITE(IOIMP,*)' CALCUL DE LA FFT DU SIGNAL ' INV=.FALSE. CALL WEXP(INV,NPOINT,W(1)) CALL FFTL(XX(1),YY(1),W(1),NPOINT) C IF(IIMPI.EQ.1) WRITE(IOIMP,*)' FFT CALCULEE ' SEGDES MLREE2 C C CREATION ET CALCUL DES LISTES DE LA DSP C IF(IIMPI.EQ.1) WRITE(IOIMP,*)' CALCUL DE LA DSP DU SIGNAL ' DUREE=DT*DBLE(NPOINT) DFREQ=1.D0/DUREE KHALF=NEXP+1 CPP + KHALF1=KHALF CPP + KMIN=INT(FRMI/DFREQ)+1 KMAX=INT(FRMA/DFREQ)+1 KDEBU=1 IF(KMIN.GT.0) KDEBU=KMIN IF((KMAX.LT.KHALF).AND.(FRMA.GT.0.)) KHALF=KMAX JG=KHALF-KDEBU+1 SEGINI MLREE1 SEGINI MLREE2 CPP ! COEF=2.D0*DFREQ*DT*DT COEF= DFREQ*DT*DT DO 20 I=KDEBU,KHALF FREQ=DBLE(I-1)*DFREQ C IF(IIMPI.EQ.1) WRITE(IOIMP,*)I,FREQ,' XX = ',XX(I) MLREE1.PROG(I-KDEBU+1)=FREQ YCOMP=CONJG(XX(I)) PMODU2=XX(I)*YCOMP cbp: on pourrait aussi ecrire PMODU2=ABS(XX(I))**2 CPP + IF(I.NE.1.AND.I.NE.KHALF1)THEN YCOMP=CONJG(XX(NPOINT-I+2)) PMODU2=PMODU2+XX(NPOINT-I+2)*YCOMP cbp: on pourrait aussi ecrire PMODU2=PMODU2+ABS(XX(NPOINT-I+2))**2 ENDIF CPP + MLREE2.PROG(I-KDEBU+1)=COEF*PMODU2 C IF(IIMPI.EQ.1) WRITE(IOIMP,*)'DSP(',I,')=',MLREE2.PROG(I) 20 CONTINUE C SEGDES MLREE1,MLREE2 C C CREATION DE L'OBJET EVOLUTIO DSP C N=1 SEGINI MEVOLL IPVO=MEVOLL TI(1:72)=TITREE IEVTEX=TI ITYEVO='REEL' SEGINI KEVOLL c KEVTEX=TI KEVTEX=KEVOL1.KEVTEX IEVOLL(1)=KEVOLL TYPX='LISTREEL' TYPY='LISTREEL' C IPROGX=MLREE1 NOMEVX=MOTX(1:12) C IPROGY=MLREE2 NOMEVY=MOTY(1:12) C NUMEVX=ICOUL NUMEVY='REEL' C SEGSUP MTRAV1,MTRAV2 SEGDES KEVOL1 SEGDES MEVOL1 C SEGDES KEVOLL,MEVOLL CALL ECROBJ('EVOLUTIO',IPVO) 666 CONTINUE RETURN END