	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
