PROGRAM FOUANM ! VERS AUG81 C C This program computes the phase and amplitude C of the fundamental frequency (k=1) using the C complex FOURIER analysis method. C C Amplitudes are normalized to 10000...65000, C phases to 0....2*pi*10000 C The phase at max. amplitude is set to +pi/2. C The pgm assumes volume-curve of LV in SA 6, Curve A. C C Revised: B.Warshaw C Hadassah University Hospital C Jerusalem, Israel. C DIMENSION RADMIN(128),IADMIN(256),B(3,42),IP(75) DIMENSION IDBLK1(4),JADMIN(256),SADMIN(128) DIMENSION IMAT(256),NM1(32),IM4(16),ICV(256) DIMENSION PEXP(24),CMFL(512) VIRTUAL PHASES(4096),AMPLIS(4096),ICURVE(256,40) C C COMPLEX PEXP,ZERO,CO,CMFL C LOGICAL*1 MADMIN(512),FILNAM(12),TITLE(6) LOGICAL*1 LADMIN(512),ELV(2),THR,BLK,NAM(12) C COMMON /TABS/ LSTP,PEXP,CMFL C EQUIVALENCE (LADMIN,IADMIN,RADMIN),(IDBLK1,JADMIN) EQUIVALENCE (MADMIN,JADMIN,SADMIN) C DATA IM4/0,256,512,768,1024,1280,1536,1792,2048,2304, 1 2560,2816,3072,3328,3584,3840/,ELV/'1','1'/ C DATA FILNAM/'S','V','A','S','V','A','R','0','6','S','Y','S'/ 1,THR/'3'/,BLK/"0/ DATA TITLE/'F','O','U','A','N','M'/ C DATA PI/3.1415927/,PI2/1.5707963/,PITWO/6.283185/ DATA CMFL/512*(0.,0.)/ C CALL CNTXSW CALL GAMSTP(JADMIN,TITLE) MAXIND=4096 DO 5 I=1,MAXIND AMPLIS(I)=0. 5 PHASES(I)=-PI ZERO=(0.,0.) CALL CMTB ! prepare CMFL variable C CALL ASSIGN (1,'SY0:FOUR.DAT',12,'NEW') CALL IRAD50(12,FILNAM,IDBLK1) ICHN1=IGETC() IF(ICHN1 .LT. 0) STOP 'NO CHNL' IF(LOOKUP(ICHN1,IDBLK1) .LT. 0) STOP 'No file SVAR06 ' IF(IREADW(256,JADMIN,0,ICHN1) .LT. 0) STOP 'Can"t read JADMIN' N1=5 IF(IREADW(256,RADMIN,N1,ICHN1) .LT. 0) STOP 'Can"t read RADMIN' CALL CLOSEC(ICHN1) C C Display curve and mark last EVEN point for analysis C WRITE(5,2) 2 FORMAT(1X,'Mark End-Diastolic Frame !') CALL FGPLOT(JADMIN,RADMIN) CALL FGPICK(NEND,Y) C WRITE(5,4) 4 FORMAT(1X,'Mark End-Systolic Frame !') CALL FGPICK(NENS,Y) C WRITE(5,10) 10 FORMAT(1X,'Mark last EVEN numbered Frame !') 25 CALL FGPICK(LSTP,Y) C IF(LSTP .LE. 44 .AND. (LSTP/2)*2 .EQ. LSTP) GO TO 15 TYPE *,'Point marked',LSTP,'is too big or ODD' GO TO 25 C 15 TYPE *,'PROCESSING TAKES APPROX 2 MIN' CALL R50ASC(12,JADMIN,NAM) WRITE(1,13) NAM 13 FORMAT(12A1) WRITE(1,12) NEND WRITE(1,12) NENS WRITE(1,12) LSTP 12 FORMAT(I5) CALL CLOSE(1) C ICHN1=IGETC() IF(ICHN1 .LT. 0) STOP 'NO CHNL' IF(LOOKUP(ICHN1,IDBLK1) .LT. 0) STOP 'No file IDBLK1 ' IF(IREADW(256,IADMIN,0,ICHN1) .LT. 0) STOP 'Can"t read IADMIN' CALL FGADM1(LADMIN,IP,B) NTYP=IP(12) NN=64 DO 20 I=1,LSTP II=I 20 NM1(II)=IFGFRM(IP,II)-1 FN=2./FLOAT(LSTP) NF2=LSTP/2 CALL ANGTB ! prepare complex PEXP C AMPMAX=0. IMAX=1 C C Loop over the number of blocks per frame C assuming type 4 matrix only C IF(NTYP .NE. 4) STOP'? FOUN : Wrong matrix type ?' C INDX=16 ! number of blocks per frame DO 100 I1=1,INDX ! loop over blocks/frame II=I1-1 IL=IM4(I1) C C NM1(i) = Pointer to first word of frame i. C DO 35 I2=1,LSTP ! loop over frames M1=NM1(I2)+II IF(IREADW(256,ICV,M1,ICHN1) .LT. 0) 1 STOP 'Read Error in CRV?' DO 36 I = 1,256 36 ICURVE(I,I2) = ICV(I) 35 CONTINUE C C Begin complex fourier analysis C NBLK=256 DO 150 I2=1,NBLK I4=IL+I2 CO=ZERO C C Indices: IL = Multiples of 256 C I4 = Running index (1-4096) C C In the following DO-LOOP, the index of CMFL(i) is C computed by the difference in the countrates of C two points separated by N/2 units. C To insure a positive non zero index 256 is added. C DO 30 I=1,NF2 30 CO=CO+PEXP(I)*CMFL(ICURVE(I2,I)-ICURVE(I2,I+NF2)+256) SI=AIMAG(CO) CI=REAL(CO) AMPLI=CABS(CO) C C What follows is similiar to previous version. C 31 IF(CI .EQ. 0.) GO TO 40 PHASE=ATAN(SI/CI) IF(CI .GT. 0.) GO TO 55 ! ALL OK PHASE=PHASE-PI IF(SI .LE. 0.) PHASE=PHASE+PITWO GO TO 55 40 PHASE=-PI IF(SI .NE. 0.) PHASE = PI2*(SI/ABS(SI)) 55 CONTINUE IF(I4 .LT. 131 .OR. I4 .GT. 3966) GO TO 45 IF(AMPLI.LE.AMPMAX) GO TO 45 ! find max amplitude AMPMAX=AMPLI IMAX=I4 45 PHASES(I4)=PHASE AMPLIS(I4)=AMPLI 150 CONTINUE 100 CONTINUE AMPMAX=AMPMAX+AMPMAX C C Normalize phases and amplitudes. C get factor for amplitudes C FACTOR=1. FMAX=650. AMPSUM=0. DO 200 I=1,3 IF(AMPMAX .GT. FMAX) GO TO 210 FMAX=FMAX/10. FACTOR=FACTOR*10. 200 CONTINUE FACTOR=2.*FACTOR C 210 PHSUM=0. DO 250 I=1,MAXIND ! phase at max amplitude = PI/2 PHASE=PHASES(I)+PI IF(PHASE .LT. 0) PHASE=PHASE+PITWO IF(PHASE .GT. PITWO) PHASE=PHASE-PITWO PHSUM=PHSUM+PHASE PHASES(I)=PHASE*1000. AMPLI=AMPLIS(I)*FACTOR AMPSUM=AMPSUM+AMPLI AMPLIS(I)=AMPLI 250 CONTINUE PHMAX=2000.*PI PHASES(4096)=PHMAX ! Uniform display up to 2 PI AMPMAX=AMPMAX*FACTOR C C Write amplitudes in SA 11, C Phases in 13 C Update SA-descriptor block C MADMIN(165)="0 ! LT=0 JADMIN(7)=1 ! matrix data JADMIN(8)=11 ! SA-number JADMIN(69)=0 ! STATIC IMAGE JADMIN(84)=MAXIND ! NBR OF WORDS IN MATRIX JADMIN(85)=0 ! WORD MODE JADMIN(87)=ISPR(AMPMAX) ! max cts JADMIN(88)=0 ! min cts RADMIN(45)=RDPR(AMPSUM) ! tot cts JADMIN(91)=INT(AMPSUM/4096.) ! av cts JADMIN(125)=MAXIND ! NMBR OF WORDS IN MATRIX JADMIN(184)=0 ! nmbr of rois=0 DO 310 I=1,10 I1=I+224 MADMIN(I1)=BLK 310 CONTINUE C ICHN1=IGETC() IF(ICHN1 .LT. 0) STOP 'NO CHNL' ENCODE(2,325,FILNAM(8)) ELV 325 FORMAT(2A1) CALL IRAD50(12,FILNAM,IDBLK1) IF(LOOKUP(ICHN1,IDBLK1) .LT. 0) STOP 'No file SVAR11 ' IF(IWRITW(256,JADMIN,0,ICHN1) .LT. 0) STOP 'Can"t write JADMIN0' C DO 330 I1=1,INDX I3=IM4(I1) DO 320 I2=1,256 I4=I3+I2 IMAT(I2)=ISPR(AMPLIS(I4)) IF(IMAT(I2).GE.0)GO TO 320 IMAT(I2)=0 320 CONTINUE IF=I1 321 IF(IWRITW(256,IMAT,IF,ICHN1) .LT. 0) STOP 'Can"t write IMAT1' 330 CONTINUE CALL CLOSEC(ICHN1) C C Write phases in SA 13 C ENCODE(1,335,FILNAM(9)) THR 335 FORMAT(A1) ICHN1=IGETC() IF(ICHN1 .LT. 0) STOP 'NO CHNL' CALL IRAD50(12,FILNAM,IDBLK1) IF(LOOKUP(ICHN1,IDBLK1) .LT. 0) STOP 'No file SVAR13 ' C JADMIN(8)=13 JADMIN(87)=ISPR(PHMAX) SADMIN(45)=RDPR(PHASUM) JADMIN(91)=INT(AMPSUM/4096.) IF(IWRITW(256,JADMIN,0,ICHN1) .LT. 0) STOP 'Can"t write JADMIN1' C DO 410 I=1,10 I1=I+224 MADMIN(I1)=BLK 410 CONTINUE DO 430 I1=1,INDX I3=IM4(I1) DO 420 I2=1,256 I4=I3+I2 PHASE=PHASES(I4) IF(PHASE .LT. 1.) PHASE=0. ! make sure phases be 0 if needed IMAT(I2)=ISPR(PHASE) 420 CONTINUE IF=I1 IF(IWRITW(256,IMAT,IF,ICHN1) .LT. 0) STOP 'Can"t write MADMIN' 430 CONTINUE CALL CLOSEC(ICHN1) C CALL GAMCA(JADMIN,-1) END