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