sw/rescue/lab8e_goettingen/disk2_11/rkb/paroff/haedof.ft

   1 C       -+-+-+-+-+ \ e  H A E D O F . F T  \ e -+-+-+-+-+
   2 C
   3 C       DO THE FFT
   4         SUBROUTINE DOFFT
   5 C
   6         INCLUDE HAEBUF.FI
   7         INCLUDE HAEGSA.FI
   8         INCLUDE HAEPTI.FI
   9         INCLUDE HAEDEF.FI
  10         INCLUDE HAEI85.FI
  11         INCLUDE HAETTY.FI
  12         INCLUDE HABRK.FI
  13 C
  14         INTEGER TOASCI,NCOND
  15         REAL INTENS,IMAXNT,COND
  16         EXTERNAL TOASCI,INTENS,IMAXNT,NCOND,COND
  17 \fC
  18         INTEGER I,J,K,PNT,BAND,         @ DO LOOP COUNTERS & LIMITS
  19      *  POSIT,POS1,                     @  START OF THE ACTUAL WINDOW (SAMRAT*SPAN)*N
  20      *  I1,I2                           @ INDEX TO COMPLEX SPECTRUM
  21         INTEGER MM,SS,MMM,SSS           @ START OF FFT, END OF FFT IN MINUTES AND SECONDS
  22         REAL X(2050)                    @ BUFFER FOR THE FFT
  23         LOGICAL L1                      @ TEMPORARY STORAGE
  24 C
  25 C       DO THE FFT
  26 C
  27         IF (REC5.LT.1) DEFINE FILE 5(MAXBL5,85,U,REC5) @ OPEN SCRATCH FILE FOR INTERMEDIATE FFT RESULTS
  28 C
  29 C       CLEAR THE OUTPUT FILE UNIT 5 TO ALLOW CNTRL C ABORT
  30 C
  31         REC5=1
  32         REDVAL=0                        @ NO VALID DATA HERE
  33         WRITE (5'REC5) (RCRD0(I),I=1,85) @ CLEAR THE FILE HEADER BLOCK
  34 C
  35 C                                       @ NXTDTA COUNTS FROM 0 TO (ENDS-1)/OVRLAP
  36         CALL MOVE (-85,0,INTE85)        @ CLEAR BUFFER FOR INTERMEDIATE RESULTS AND NXTDTA TOO
  37         DO 50 POS1=BEGIN,ENDS-1,OVRLAP  @ LOOP FOR EACH FFT
  38         POSIT=POS1*SAMRAT
  39         IF (POSIT.GT.ENDS*SAMRAT-INCR) GOTO 50  @ SKIP LAST OVERLAPPING FFT
  40         IF BREAK(11) GOTO 51            @ SW 11 ABORTS TASK
  41         ARTEFK=0                        @ COUNTER FOR OVERRANGE SIGNAL
  42         DO 40 I=1,INCR,1                @ LOOP TO FETCH EACH VALUE FROM INPUT FILE
  43         X(I)=ESAM(POSIT+I-1)
  44 40      ARTEFK=NCOND(ABS(X(I)).GT.509,ARTEFK+1,ARTEFK) @ COUNT NUMBER OF POINTS OUT OF RANGE
  45         CALL FFTC(X,INCR/2,EXPON-1,1.)  @ BLACK MAGIC BOXES
  46         CALL FFTR(X,INCR/2,1.,1.)
  47         CALL HANING (X,INCR/2)          @ SMOOTH REAL PART
  48         CALL HANING (X(INCR/2),INCR/2)  @ SMOOTH IMAG PART OF FFT
  49         CALL FTPOWR(X,INCR)             @ COMPUTE THE POWER SPECTRUM AND THE MINIMUM OF THE SPECTRUM
  50         X(1)=0                          @ CLEAR THE FIRST AND SECOND CHANNEL TO
  51         X(2)=0                          @  STRIP OFF THE DC OFFSET
  52 C
  53 C                                       @ IF SAMPLE RATE 64 HZ THEN COMPUTE THE MAXIMUM OF THE SPECTRUM STARTING AT 1. HZ
  54 C                                       @  ELSE COMPUTE THE MAXIMUM STARTING AT 0. HZ BUT EXCLUDING DC (CHANNEL 0 & 1)
  55         CALL IMAXNT(X,COND(SAMRAT.EQ.64,1.,0.),SAMRAT/2.,K) @ FREQUENCY OF MAXIMUM, INDEX TO INTENSITY X --> K
  56         XMAXI=X(K)                      @ INTENSITY OF MAXIMUM
  57         PNT=MOD(NXTDTA,4)+1             @ POINTER TO NEXT INTE LOCATION WHERE TO INSERT BAND INTENSITY DATA
  58         DO 45 BAND=1,5                  @ FIVE BANDS
  59         L1=SAMRAT/2.GE.FREQU(BAND,2)    @ COMPUTE THE INTENSITY ONLY FROM THOSE FREQUENCY
  60 C                                       @  RANGES WHICH ARE LESS THAN THE HIGHEST FREQUENCY WE MEASURED
  61         INTE(1,BAND,PNT)=0
  62         IF (L1) INTE(1,BAND,PNT)=INTENS(X,FREQU(BAND,1),FREQU(BAND,2)) @ TOTAL INTENSITY OF THIS BAND
  63         IF (L1) INTE(2,BAND,PNT)=IMAXNT(X,FREQU(BAND,1),FREQU(BAND,2),K) @ FREQUENCY OF MAXIMUM [HZ]
  64         INTE(3,BAND,PNT)=POS1            @ TIME OF THIS POWER SPECTRUM
  65         BEGFFT=INTE(3,BAND,PNT)         @ SAVE THE TIME OF THIS POWER SPECTRUM INTO THE FILE HEADER
  66 CX      IF (L1) INTE(3,BAND,PNT)=X(K)   @ INTENSITY OF FREQUENCY MAXIMUM
  67         INTE(4,BAND,PNT)=INTE(1,BAND,PNT)/AMAX1(1.,INTE(1,1,PNT)) @ RELATIVE INTENSITY OF THIS BAND
  68 45      CONTINUE
  69         MM=POS1/60                      @ MINUTES
  70         SSS=POS1+SPAN                   @ SECONDS FFT ENDS
  71         SS=TOASCI(MOD(POS1,60))
  72         MMM=SSS/60
  73         SSS=TOASCI(MOD(SSS,60))
  74         WRITE (3,1) MM,SS,MMM,SSS,((INTE(I,J,PNT),I=1,2),J=1,5),ARTEFK
  75         REC5=STSCAN+NXTDTA/4            @ 2 BLOCK HEADER, 25 BLOCKS POWER SPECTRUM AND THEN WE STORE THE INTE DATA SETS
  76         NXTDTA=NXTDTA+1                 @ WE STORE 4 INTE DATA SET IN ONE BLOCK
  77         WRITE (5'REC5) INTE85           @ SAVE THIS DATA ONTO THE SCRATCH FILE
  78 50      CONTINUE                        @ END OF THE DO LOOP FOR EACH FFT
  79 \f51     CONTINUE                        @ IF USER ABORTS TASK THEN ENTRY HERE
  80         COMP=.TRUE.                     @ COMPUTATIONS DONE FLAG
  81         REC5=1                          @ SAVE THE LAST FFT DONE
  82         REDVAL=6H2DPWLD
  83         WRITE (5'REC5) (RCRD0(I),I=1,85)@ INSERT THE FILE HEADER BLOCK
  84         REDVAL=0                        @ FOR SAVETY ONLY
  85         SETNU5=SETNUM                   @ RECORD # OF THE ACTUAL DATA SET --> HEADER OF UNIT 5
  86         CHANE5=CHANEL                   @ AND WE SAVE THE CHANNEL NUMBER TOO
  87         DATCNT=NXTDTA                   @ INSERT THE NUMBER OF SCANS INTO THE HEADER BLOCK
  88         FILTER=0                        @ NEW DATA, NOTHING FILTERED
  89         CONTNS=0                        @ REC3 TO 16 CONTAINS A EEG POWER SPECTRUM
  90         WRITE (5'REC5) (LCEGSA(I),I=1,LCEGSA(1)) @ WRITE THE DATA HEADER BLOCK ( COMMON CEGESA)
  91         DO 55 J=1,INCR/2+84,85          @ SAVE THE LAST POWER SPECTRUM COMPUTED INTO REC 3-27
  92 55      WRITE (5'REC5) (X(I),I=J,MIN0(1024,J+84)) @ 13 RECORDS TO WRITE
  93         COMP=.TRUE.                     @ SOME COMPUTATIONS DONE
  94         RETURN
  95 C
  96 C
  97 1       FORMAT (I5,1H:,A2,1H-,I3,1H:,A2,5(2X,1P,E8.2,0P,F6.1,4X),I5)
  98         END
  99 \1a