C	subroutines mainly from the book by Silvia & Robinson

C	note that the vectors (because of Fortran) go from index 1 (ej 0)
C	cf. the middle of p. 185 in the book





	SUBROUTINE ZERO (LX, X)
C	zero is to store zero in an array
	REAL X(*)
	IF (LX) 30,30,10
10	DO 20 I=1,LX
20	X(I)=0.0
30	RETURN
	END





	SUBROUTINE DOT (L,X,Y, ANS)
C	dot is dot product
	REAL X(*),Y(*)
	ANS=0.0
	IF (L) 30,30,10
10	DO 20 I=1,L
20	ANS=ANS+X(I)*Y(I)
30	RETURN
	END





	SUBROUTINE FOLD (LA,A,LB,B, LC,C)
C	fold is folding (convolution)
C	note: lc and c are output, cf. foldny below
	REAL A(*),B(*),C(*)
	LC=LA+LB-1
	CALL ZERO(LC,C)
	DO 10 I=1,LA
	 DO 10 J=1,LB
	  K=I+J-1
10	C(K)=A(I)*B(J)+C(K)
	RETURN
	END





	SUBROUTINE FOLDNY (LA,A,LB,B,LC, C)
C	fold is folding (convolution)
C	note: only c is output, cf. fold above
	REAL A(*),B(*),C(*)
	CALL ZERO(LC,C)
	DO 10 I=1,LA
	 DO 10 J=1,MIN(LB,LC+1-I)
	  K=I+J-1
10	C(K)=A(I)*B(J)+C(K)
	RETURN
	END





	SUBROUTINE CROSS (LX,X,LY,Y,LC, C)
C	cross is cross product (cross-correlation)
	REAL X(*),Y(*),C(*)
	DO 10 I=1,LC
10	CALL DOT(MIN0(LY+I-1,LX)-I+1,X(I),Y,C(I))
	RETURN
	END





	SUBROUTINE POLYDV (N,DVS,M,DVD,L, Q)
C	perform division of power series of the form
C	(dvd(1)+x*dvd(2)+x*x*dvd(3)+...)/
C	(dvs(1)+x*dvs(2)+x*x*dvs(3)+...)
C	the result is = q(1)+x*q(2)+x*x*q(3)+...
C	n is the number of coefficients in the divisor polynomial
C	m is the number of coefficients in the dividend polynomial
C	l is the number of coefficients in the quotient polynomial
C	dvs is the input divisor polynomial coefficients
C	dvd is the input dividend polynomial coefficients
C	q is the output quotient polynomial coefficients
C	dvs(1) must be nonzero
C	this could be programmed to allow equivalence (dvd,q)
	REAL DVS(*),DVD(*),Q(*)
	NM=N-1
	DO 8 I=1,L
8	Q(I)=0.0
C	move the used portion of dvd to q
	MML=MIN0(M,L)
	DO 10 I=1,MML
10	Q(I)=DVD(I)
	DO 50 I=1,L
	 Q(I)=Q(I)/DVS(1)
	 IF (I-L) 30,20,30
20	 RETURN
30	 K=I
	 ISUB=MIN0(NM,L-I)
	 DO 40 J=1,ISUB
	  K=K+1
	  Q(K)=Q(K)-Q(I)*DVS(J+1)
40	 CONTINUE
50	CONTINUE
C	program never gets here
	RETURN
	END





	SUBROUTINE PEO (LR,R, A)
C	peo is the auxiliary Toeplitz recursion
C	it gives the prediction-error operator a(i), a(1)=1.0
	REAL R(LR),A(LR)
	V=R(1)
	A(1)=1.0
	IF (LR.EQ.1) RETURN
	DO 6 L=2,LR
	 D=0.0
	 L3=L-1
	 DO 1 J=1,L3
	  K=L-J+1
1	 D=D+A(J)*R(K)
	 C=D/V
	 IF (L.EQ.2) GOTO 4
	 L1=(L-2)/2
	 L2=L1+1
	 IF (L2.LT.2) GOTO 3
	 DO 2 J=2,L2
	  HOLD=A(J)
	  K=L-J+1
	  A(J)=A(J)-C*A(K)
2	 A(K)=A(K)-C*HOLD
	 IF (2*L1.EQ.L-2) GOTO 4
3	 LT3=L2+1
	 A(LT3)=A(LT3)-C*A(LT3)
4	 A(L)=-C
6	V=V-C*D
	RETURN
	END





	SUBROUTINE FFT (X,Y, TABLE,M,LL,ISN)
C	fft is in place dft computation using Sande algorithm
C	   and Markel pruning modification
C	x = array of length 2**m used to hold real part of complex input
C	y = array of length 2**m used to hold imaginary part of complex input
C	table = array of length n/4+1 (n=2**m) , contains
C	   quarter-length cosine table
C	m = integer exponent to two which determines size (2**m) of 
C	   fft to be performed  (bit reverse table is set for a
C	   maximum of n=2**16=65536)
C	ll = integer exponent to two which determines the number (2**ll) 
C	   actual data points, number of stages in which no pruning is 
C	   allowable
C	isn :  set to -1 for forward fft (dft)
C	       set to +1 for inverse fft (idft)
C	x = real part of complex output
C	y = imaginary part of complex output
C	
C	observera : vektorerna gaar ju fraan 1 (ej 0) och konventionen med
C	   shift (se oeverst i denna fil) bibehaallen, saa vid dft blir 
C	   & "in-x(1)" vaerde vid tiden 0
C	   & "in-y(1)" vaerde vid tiden 0
C	   & "ut-x(1)" 0-frekvens (dvs dc) komponent
C	   & "ut-y(1)" 0-frekvens (dvs dc) komponent


C	this is exactly what happens to the x,y arrays:

C	for k=1,2,...,2**m
C	  x(out)(k) + i*y(out)(k)   :=   SUM (n=1,2**m)  
C	             (x(in)(n)+i*y(in)(n)) * exp(isn*2*pi*i*(n-1)*(k-1)/2**m)
C	end
C	if (isn.gt.0) then all elements of the x,y arrays are finally
C	                   divided by 2**m


	REAL X(*),Y(*),TABLE(*)
	INTEGER L(16)
	EQUIVALENCE (L16,L(1)),(L15,L(2)),(L14,L(3)),(L13,L(4)),
     $	(L12,L(5)),(L11,L(6)),(L10,L(7)),(L9,L(8)),(L8,L(9)),
     $	(L7,L(10)),(L6,L(11)),(L5,L(12)),(L4,L(13)),(L3,L(14)),
     $	(L2,L(15)),(L1,L(16))
	IF (M.GT.16) THEN
	  WRITE(*,*) ' %%% FFT reverse table set too small'
	  STOP 1616
	END IF
	N=2**M
	ND4=N/4
	ND4P1=ND4+1
	ND4P2=ND4P1+1
	ND2P2=ND4+ND4P2
	LLL=2**LL
	DO 8 LO=1,M
	 LMX=2**(M-LO)
	 LMM=LMX
	 LIX=2*LMX
	 ISCL=N/LIX
C	 test for pruning
	 IF (LO-M+LL) 1,2,2
1	 LMM=LLL
2	 DO 8 LM=1,LMM
	  IARG=(LM-1)*ISCL+1
	  IF (IARG.LE.ND4P1) GOTO 4
	   K1=ND2P2-IARG
	   C=-TABLE(K1)
	   K3=IARG-ND4
	   S=ISN*TABLE(K3)
	   GOTO 6
4	  C=TABLE(IARG)
	   K2=ND4P2-IARG
	   S=ISN*TABLE(K2)
6	  CONTINUE
	  DO 8 LI=LIX,N,LIX
	   J1=LI-LIX+LM
	   J2=J1+LMX
	   T1=X(J1)-X(J2)
	   T2=Y(J1)-Y(J2)
	   X(J1)=X(J1)+X(J2)
	   Y(J1)=Y(J1)+Y(J2)
	   X(J2)=C*T1-S*T2
	   Y(J2)=C*T2+S*T1
8	CONTINUE

C	perform bit reversal

	DO 40 J=1,16
	 L(J)=1
	 IF (J-M) 31,31,40
31	 L(J)=2**(M+1-J)
40	CONTINUE
	JN=1
	DO 60 J1=1,L1
	 DO 60 J2=J1,L2,L1
	  DO 60 J3=J2,L3,L2
	   DO 60 J4=J3,L4,L3
	    DO 60 J5=J4,L5,L4
	     DO 60 J6=J5,L6,L5
	      DO 60 J7=J6,L7,L6
	       DO 60 J8=J7,L8,L7
	        DO 60 J9=J8,L9,L8
	         DO 60 J10=J9,L10,L9
	          DO 60 J11=J10,L11,L10
	           DO 60 J12=J11,L12,L11
	            DO 60 J13=J12,L13,L12
	             DO 60 J14=J13,L14,L13
	              DO 60 J15=J14,L15,L14
	               DO 60 JR=J15,L16,L15
	                IF (JN-JR) 51,51,54
51	                R=X(JN)
	                X(JN)=X(JR)
	                X(JR)=R
	                FI=Y(JN)
	                Y(JN)=Y(JR)
	                Y(JR)=FI
54	                IF (ISN) 53,53,52
52	                X(JR)=X(JR)/FLOAT(N)
	                Y(JR)=Y(JR)/FLOAT(N)
53	                JN=JN+1
60	CONTINUE
	RETURN
	END





	SUBROUTINE COSQT (M, TABLE)
C	cosqt generates quarter-length cosine table
	REAL TABLE(*)
	N=2**M
	ND4P1=N/4+1
	SCL=6.283185307/FLOAT(N)
	DO 10 I=1,ND4P1
	 ARG=FLOAT(I-1)*SCL
10	TABLE(I)=COS(ARG)
	RETURN
	END





	SUBROUTINE CENV (X,Y, TABLE,M)

C	cenv aer rutin foer diskret "complex envelope" utnyttjande fft-
C	   rutinen ovan, jmf Farnbach i BSSA 65:4 (1975) sid. 952

C	x = array of length 2**m used to hold real part of complex input
C	y = array of length 2**m used to hold imaginary part of complex input
C	table = array of length n/4+1 (n=2**m) , contains
C	   quarter-length cosine table
C	m = integer exponent to two which determines size (2**m) of 
C	   cenv to be performed
C	x = real part of complex output
C	y = imaginary part of complex output

C	note: the accuracy may be increased by extending "input x&y" with 
C	zeroes, since the frequency sampling of the DFT will then be 
C	more dense

	REAL X(*),Y(*),TABLE(*)
	CALL FFT(X,Y,TABLE,M,M,-1)
	N=2**M
	ND2=N/2
	ND2P2=ND2+2
	DO K=2,ND2
	 X(K)=2.0*X(K)
	 Y(K)=2.0*Y(K)
	END DO
	DO K=ND2P2,N
	 X(K)=0.0
	 Y(K)=0.0
	END DO
	CALL FFT(X,Y,TABLE,M,M,+1)
	RETURN
	END











	SUBROUTINE CFFT (CXY, TABLE,M,LL,ISN)
C	cfft is in place dft computation using Sande algorithm
C	   and Markel pruning modification
C	   (complex input/output)
C	cxy = array of length 2**m used to hold complex input
C	table = array of length n/4+1 (n=2**m) , contains
C	   quarter-length cosine table
C	m = integer exponent to two which determines size (2**m) of 
C	   cfft to be performed  (bit reverse table is set for a
C	   maximum of n=2**16=65536)
C	ll = integer exponent to two which determines the number (2**ll) 
C	   actual data points, number of stages in which no pruning is 
C	   allowable
C	isn :  set to -1 for forward cfft (dft)
C	       set to +1 for inverse cfft (idft)
C	cxy = complex output
	
C	observera : vektorerna gaar ju fraan 1 (ej 0) och konventionen med
C	   shift (se oeverst i denna fil) bibehaallen, saa vid dft blir 
C	   & "in-cxy(1)" vaerde vid tiden 0
C	   & "ut-cxy(1)" 0-frekvens (dvs dc) komponent

C	this is exactly what happens to the cxy array:

C	for k=1,2,...,2**m
C	  cxy(out)(k) :=   SUM (n=1,2**m)  
C	             cxy(in)(n) * exp(isn*2*pi*i*(n-1)*(k-1)/2**m)
C	end
C	if (isn.gt.0) then all elements of the cxy array are finally
C	                   divided by 2**m


	COMPLEX CXY(*)
	REAL TABLE(*)
	INTEGER L(16)
	COMPLEX CT12,CR
	EQUIVALENCE (L16,L(1)),(L15,L(2)),(L14,L(3)),(L13,L(4)),
     $	(L12,L(5)),(L11,L(6)),(L10,L(7)),(L9,L(8)),(L8,L(9)),
     $	(L7,L(10)),(L6,L(11)),(L5,L(12)),(L4,L(13)),(L3,L(14)),
     $	(L2,L(15)),(L1,L(16))
	IF (M.GT.16) THEN
	  WRITE(*,*) ' %%% CFFT reverse table set too small'
	  STOP 1616
	END IF
	N=2**M
	ND4=N/4
	ND4P1=ND4+1
	ND4P2=ND4P1+1
	ND2P2=ND4+ND4P2
	LLL=2**LL
	DO 8 LO=1,M
	 LMX=2**(M-LO)
	 LMM=LMX
	 LIX=2*LMX
	 ISCL=N/LIX
C	 test for pruning
	 IF (LO-M+LL) 1,2,2
1	 LMM=LLL
2	 DO 8 LM=1,LMM
	  IARG=(LM-1)*ISCL+1
	  IF (IARG.LE.ND4P1) GOTO 4
	   K1=ND2P2-IARG
	   C=-TABLE(K1)
	   K3=IARG-ND4
	   S=ISN*TABLE(K3)
	   GOTO 6
4	  C=TABLE(IARG)
	   K2=ND4P2-IARG
	   S=ISN*TABLE(K2)
6	  CONTINUE
	  DO 8 LI=LIX,N,LIX
	   J1=LI-LIX+LM
	   J2=J1+LMX
	   CT12=CXY(J1)-CXY(J2)
	   CXY(J1)=CXY(J1)+CXY(J2)
	   CXY(J2)=CMPLX(C,S)*CT12
8	CONTINUE

C	perform bit reversal

	DO 40 J=1,16
	 L(J)=1
	 IF (J-M) 31,31,40
31	 L(J)=2**(M+1-J)
40	CONTINUE
	JN=1
	DO 60 J1=1,L1
	 DO 60 J2=J1,L2,L1
	  DO 60 J3=J2,L3,L2
	   DO 60 J4=J3,L4,L3
	    DO 60 J5=J4,L5,L4
	     DO 60 J6=J5,L6,L5
	      DO 60 J7=J6,L7,L6
	       DO 60 J8=J7,L8,L7
	        DO 60 J9=J8,L9,L8
	         DO 60 J10=J9,L10,L9
	          DO 60 J11=J10,L11,L10
	           DO 60 J12=J11,L12,L11
	            DO 60 J13=J12,L13,L12
	             DO 60 J14=J13,L14,L13
	              DO 60 J15=J14,L15,L14
	               DO 60 JR=J15,L16,L15
	                IF (JN-JR) 51,51,54
51	                CR=CXY(JN)
	                CXY(JN)=CXY(JR)
	                CXY(JR)=CR
54	                IF (ISN) 53,53,52
52	                CXY(JR)=CXY(JR)/N
53	                JN=JN+1
60	CONTINUE
	RETURN
	END