DECLARE FUNCTION atan2! (r!, i!)
DECLARE SUB fft (xr!(), yi!(), n%, inverse%)
DECLARE SUB FFT2DCRCalc (xreal!(), yimag!(), c%, r%, flag%)
DECLARE SUB FFT2DRCCalc (xreal!(), yimag!(), c%, r%, flag%)
DECLARE SUB fft2 (xr!(), y!(), n%, m%, ks%)
DECLARE SUB reorder (xr!(), y!(), n%, m%, ks%, reel%)
DECLARE SUB trans (xr!(), y!(), n%, eval%)
DECLARE SUB rfft2 (xr!(), y!(), n%, m%, ks%)
DECLARE SUB WindowData (x!(), numdat%, wind%)
DECLARE SUB FFTCalc (xreal(), yimag(), numdat%)
DECLARE SUB FFTInvCalc (xreal(), yimag(), numdat%)
DECLARE FUNCTION Parzen! (n%, j%)
DECLARE FUNCTION Hanning! (n%, j%)
DECLARE FUNCTION Hamming! (n%, j%)
DECLARE FUNCTION ExactBlackman! (n%, j%)
DECLARE FUNCTION Welch! (n%, j%)
DECLARE SUB FFTSolve (xreal(), yimag(), numdat%, flag%)
DECLARE FUNCTION SquareAndSum! (a!, b!)
DECLARE FUNCTION FFTMagnitude! (xr(), yi(), n%, i%)
DECLARE FUNCTION FFTPhase! (xr(), yi(), n%, i%)
DECLARE FUNCTION FFTFrequency! (n%, samplefreq, i%)

CONST pi = 3.141592653589793#

FUNCTION atan2 (r, i)
CONST NearZero = 1E-20
DIM TempAngle
    IF ABS(r) < NearZero THEN
      IF ABS(i) < NearZero THEN TempAngle = 0!
      IF i < -(NearZero) THEN TempAngle = 3! * pi / 2!
      IF i > NearZero THEN TempAngle = pi / 2!
      atan2 = TempAngle
      EXIT FUNCTION
    ELSE
      IF ABS(i) < NearZero THEN
	IF (r < -(NearZero)) THEN TempAngle = pi
	IF r > NearZero THEN TempAngle = 0
	atan2 = TempAngle
	EXIT FUNCTION
      ELSE
	IF (r > NearZero) AND (i > NearZero) THEN TempAngle = ATN(i / r)
	IF (r < -NearZero) AND (i > NearZero) THEN TempAngle = pi - ATN(-i / r)
	IF (r < -NearZero) AND (i < -NearZero) THEN TempAngle = pi + ATN(i / r)
	IF (r > NearZero) AND (i < -NearZero) THEN TempAngle = 2! * pi - ATN(-i / r)
      END IF
    END IF
    atan2 = TempAngle
END FUNCTION 'CAngle

SUB Convolve (filtcoef(), x(), y(), filtLen%, aLen%)
DIM n%, m%
DIM summ

  FOR n% = 0 TO (aLen% + filtLen% - 2)
    summ = 0!
    FOR m% = 0 TO filtLen% - 1
	IF ((n% - m%) >= 0) AND ((n% - m%) <= (aLen% - 1)) THEN
	  summ = summ + filtcoef(m%) * x(n% - m%)
	END IF
    NEXT
    y(n%) = summ
  NEXT
END SUB

FUNCTION ExactBlackman (n%, j%)
  ExactBlackman = .42 - .5 * COS(2! * pi * j% / (n% - 1)) + .08 * COS(4! * pi * j% / (n% - 1))
END FUNCTION

SUB fft (xr(), yi(), n%, inverse%)
DIM nn%, m%, j%
DIM pp, qq

      m% = -1
      nn% = n%

      WHILE (nn% > 0)
	nn% = nn% \ 2
	m% = m% + 1
      WEND

      CALL fft2(xr(), yi(), n%, m%, n%)
      CALL reorder(xr(), yi(), n%, m%, n%, 0)

      IF (inverse% = 1) THEN
	pp = 1! / (n%)
	qq = -pp
	FOR j% = 0 TO n% - 1
	  xr(j%) = xr(j%) * pp
	  yi(j%) = yi(j%) * qq
	NEXT
      ELSE
	FOR j% = 0 TO n% - 1
	  yi(j%) = -yi(j%)
	NEXT
      END IF
END SUB

SUB fft2 (xr(), y(), n%, m%, ks%)
DIM cc%(m%)
DIM indx%, k0%, k1%, k2%, k3%, span%, j%, k%, kb%, kn%, fp%
DIM mm%, mk%, breakf%, bf2%
DIM rad, c1, c2, c3, s1, s2, s3, ck, sk, sq
DIM x0, x1, x2, x3, y0, y1, y2, y3

      sq = .707106781187#
      sk = .382683432366#
      ck = .92387953251#

      cc%(m%) = ks%
      mm% = (m% \ 2) * 2
      kn% = 0
      FOR k% = m% - 1 TO 0 STEP -1
	cc%(k%) = cc%(k% + 1) \ 2
      NEXT
      rad = 6.28318530718# / ((cc%(0)) * (ks%))
      mk% = m% - 5
      breakf% = 0
      WHILE ((kn% < n%) AND (breakf% = 0))

999     breakf% = 0
	kb% = kn%
	kn% = kn% + ks%

	IF (mm% <> m%) THEN
	  k2% = kn%
	  k0% = cc%(mm%) + kb%
	  WHILE (k0% > kb%)
	    k2% = k2% - 1
	    k0% = k0% - 1
	    x0 = xr(k2%)
	    y0 = y(k2%)
	    xr(k2%) = xr(k0%) - x0
	    xr(k0%) = xr(k0%) + x0
	    y(k2%) = y(k0%) - y0
	    y(k0%) = y(k0%) + y0
	  WEND
	END IF

	c1 = 1!
	s1 = 0!
	indx% = 0
	k% = mm% - 2
	j% = 3

	IF (k% >= 0) THEN
	   GOTO 1002
	ELSE
	   IF (kn% < n%) THEN
	     GOTO 999
	   ELSE
	     breakf% = 1
	     GOTO 1008
	   END IF
	END IF
1001    bf2% = 1
	WHILE (bf2% = 1)
	  IF (cc%(j%) <= indx%) THEN
	    indx% = indx% - cc%(j%)
	    j% = j% - 1
	    IF (cc%(j%) <= indx%) THEN
	      indx% = indx% - cc%(j%)
	      j% = j% - 1
	      k% = k% + 2
	    END IF
	  ELSE
	    bf2% = 0
	  END IF
	WEND
	indx% = cc%(j%) + indx%
	j% = 3

1002    span% = cc%(k%)
	fp% = 0
	IF (indx% <> 0) THEN
	  c2 = (indx%) * (span%) * rad
	  c1 = COS(c2)
	  s1 = SIN(c2)
	END IF
1005    IF ((fp% = 1) OR (indx% <> 0)) THEN
	   c2 = c1 * c1 - s1 * s1
	   s2 = 2! * s1 * c1
	   c3 = c2 * c1 - s2 * s1
	   s3 = c2 * s1 + s2 * c1
	END IF
	FOR k0% = kb% + span% - 1 TO kb% STEP -1
	  k1% = k0% + span%
	  k2% = k1% + span%
	  k3% = k2% + span%
	  x0 = xr(k0%)
	  y0 = y(k0%)
	  IF (s1 = 0!) THEN
	    x1 = xr(k1%)
	    y1 = y(k1%)
	    x2 = xr(k2%)
	    y2 = y(k2%)
	    x3 = xr(k3%)
	    y3 = y(k3%)
	  ELSE
	    x1 = xr(k1%) * c1 - y(k1%) * s1
	    y1 = xr(k1%) * s1 + y(k1%) * c1
	    x2 = xr(k2%) * c2 - y(k2%) * s2
	    y2 = xr(k2%) * s2 + y(k2%) * c2
	    x3 = xr(k3%) * c3 - y(k3%) * s3
	    y3 = xr(k3%) * s3 + y(k3%) * c3
	  END IF
	  xr(k0%) = x0 + x2 + x1 + x3
	  y(k0%) = y0 + y2 + y1 + y3
	  xr(k1%) = x0 + x2 - x1 - x3
	  y(k1%) = y0 + y2 - y1 - y3
	  xr(k2%) = x0 - x2 - y1 + y3
	  y(k2%) = y0 - y2 + x1 - x3
	  xr(k3%) = x0 - x2 + y1 - y3
	  y(k3%) = y0 - y2 - x1 + x3
	NEXT
	IF (k% > 0) THEN
	  k% = k% - 2
	  GOTO 1002
	END IF
	kb% = k3% + span%
	IF (kb% < kn%) THEN
	  IF (j% = 0) THEN
	    k% = 2
	    j% = mk%
	    GOTO 1001
	  END IF
	  j% = j% - 1
	  c2 = c1
	  IF (j% = 1) THEN
	    c1 = c1 * ck + s1 * sk
	    s1 = s1 * ck - c2 * sk
	  ELSE
	    c1 = (c1 - s1) * sq
	    s1 = (s1 + c2) * sq
	  END IF
	  fp% = 1
	  GOTO 1005
	END IF

1008  ' CONTINUE

      WEND

END SUB

SUB FFT2DCalc (xreal(), yimag(), c%, r%, flag%)
 IF (flag% = 0) THEN
   CALL FFT2DRCCalc(xreal(), yimag(), c%, r%, flag%)
 ELSE
   CALL FFT2DCRCalc(xreal(), yimag(), c%, r%, flag%)
 END IF
END SUB

SUB FFT2DCRCalc (xreal(), yimag(), c%, r%, flag%)
DIM i%, a%, j%
DIM xVector(c% - 1)
DIM yVector(c% - 1)

      FOR j% = 0 TO r% - 1
	   FOR a% = 0 TO c% - 1
		xVector(a%) = xreal(j%, a%)
		yVector(a%) = yimag(j%, a%)
	   NEXT
	   IF (flag% = 0) THEN
		CALL FFTCalc(xVector(), yVector(), c%)
	   ELSE
		CALL FFTInvCalc(xVector(), yVector(), c%)
	   END IF
	   FOR a% = 0 TO c% - 1
	       xreal(j%, a%) = xVector(a%)
	       yimag(j%, a%) = yVector(a%)
	   NEXT
      NEXT

      FOR i% = 0 TO c% - 1
	   FOR a% = 0 TO r% - 1
		xVector(a%) = xreal(a%, i%)
		yVector(a%) = yimag(a%, i%)
	   NEXT
	   IF (flag% = 0) THEN
		CALL FFTCalc(xVector(), yVector(), r%)
	   ELSE
		CALL FFTInvCalc(xVector(), yVector(), r%)
	   END IF
	   FOR a% = 0 TO r% - 1
		xreal(a%, i%) = xVector(a%)
		yimag(a%, i%) = yVector(a%)
	   NEXT
      NEXT
END SUB

SUB FFT2DRCCalc (xreal(), yimag(), c%, r%, flag%)
DIM i%, a%, j%
DIM xVector(r% - 1), yVector(r% - 1)

      FOR j% = 0 TO c% - 1
	   FOR a% = 0 TO r% - 1
		xVector(a%) = xreal(a%, j%)
		yVector(a%) = yimag(a%, j%)
	   NEXT
	   IF (flag% = 0) THEN
		CALL FFTCalc(xVector(), yVector(), c%)
	   ELSE
		CALL FFTInvCalc(xVector(), yVector(), c%)
	   END IF
	   FOR a% = 0 TO r% - 1
	       xreal(a%, j%) = xVector(a%)
	       yimag(a%, j%) = yVector(a%)
	   NEXT
      NEXT

      FOR i% = 0 TO r% - 1
	   FOR a% = 0 TO c% - 1
		xVector(a%) = xreal(i%, a%)
		yVector(a%) = yimag(i%, a%)
	   NEXT
	   IF (flag% = 0) THEN
		CALL FFTCalc(xVector(), yVector(), r%)
	   ELSE
		CALL FFTInvCalc(xVector(), yVector(), r%)
	   END IF
	   FOR a% = 0 TO c% - 1
		xreal(i%, a%) = xVector(a%)
		yimag(i%, a%) = yVector(a%)
	   NEXT
      NEXT
END SUB

SUB FFTCalc (xreal(), yimag(), numdat%)
   CALL fft(xreal(), yimag(), numdat%, 0)
END SUB

FUNCTION FFTFrequency (n%, samplefreq, i%)
DIM in

      in = i% / n%
      IF (i% <= n% / 2) THEN
	result = samplefreq * in
      ELSE
	i% = n% - i%
	result = -samplefreq * in
      END IF
      FFTFrequency = result
END FUNCTION

SUB FFTInvCalc (xreal(), yimag(), numdat%)
   CALL fft(xreal(), yimag(), numdat%, 1)
END SUB

FUNCTION FFTMagnitude (xr(), yi(), n%, i%)
DIM result, nr, rr, ii

      nr = n%
      IF (i% > n% / 2) THEN
	result = 0!
      ELSEIF (i% = 0) THEN
	 result = SQR(xr(0) * xr(0) + yi(0) * yi(0)) / nr
      ELSE
	  rr = ABS(xr(i%)) + ABS(xr(n% - i%))
	  ii = ABS(yi(i%)) + ABS(yi(n% - i%))
	  result = SQR(rr * rr + ii * ii) / nr
      END IF
      FFTMagnitude = result
END FUNCTION

FUNCTION FFTPhase (xr(), yi(), n%, i%)
DIM result
      result = atan2(xr(i%), yi(i%))
      FFTPhase = result
END FUNCTION

SUB FFTSwap (x, y)
DIM temp

      temp = x
      x = y
      y = temp
END SUB

SUB FIRFreqSample (filtcoef(), fp, n%, dc%, filtType%, win%)
DIM np%, i%, j%, m%, m1%, n2%
DIM lowbin, highbin, xt, am, q
DIM a(n% - 1)

      m% = (n% + 1) \ 2
      am = (n% + 1) / 2!
      m1% = n% \ 2 + 1
      q = 2! * pi / n%
      n2% = n% \ 2
      IF (dc% = 1) THEN
	  np% = INT(n% * fp + .5#)
      ELSE
	  np% = INT(n% * fp + 1#)
      END IF
      IF (filtType% = 1) THEN
	  lowbin = 0#
	  highbin = 1#
      ELSE
	  lowbin = 1#
	  highbin = 0#
      END IF

      FOR j% = 0 TO np%
	  a(j%) = lowbin
      NEXT
      FOR j% = np% + 1 TO m1%
	  a(j%) = highbin
      NEXT

      IF (dc% = 0) THEN
	   FOR j% = 1 TO m%
	       xt = a(1) / 2!
	       FOR i% = 2 TO m%
		    xt = xt + a(i%) * COS(q * (am - j%) * (i% - 1))
	       NEXT
	       filtcoef(j% - 1) = 2! * xt / n%
	       filtcoef(n% - j%) = filtcoef(j% - 1)
	   NEXT
      ELSE
	   FOR j% = 1 TO m%
		xt = 0!
		FOR i% = 1 TO n2%
		    xt = xt + a(i%) * COS(q * (am - j%) * (i% - .5))
		NEXT
		IF (am <> m%) THEN xt = xt + a(m%) * COS(pi * (am - j%)) / 2!
		filtcoef(j% - 1) = 2! * xt / (n%)
		filtcoef(n% - j%) = filtcoef(j% - 1)
	   NEXT
      END IF
      IF (win% > 0) THEN CALL WindowData(filtcoef(), n%, win%)
END SUB

SUB FreqResponse (filtcoef(), a(), n%, k%)
DIM i%, j%, m%, n2%
DIM att, am, q

      q = pi / k%
      am = (n% + 1) / 2!
      m% = (n% + 1) \ 2
      n2% = n% \ 2
      FOR j% = 1 TO k%
	  att = 0!
	  IF (am = m%) THEN att = filtcoef(m% - 1) / 2!
	  FOR i% = 1 TO n2%
		att = att + filtcoef(i% - 1) * COS(q * (am - i%) * (j% - 1))
	  NEXT
	  a(j% - 1) = 2! * att
      NEXT
END SUB

FUNCTION Hamming (n%, j%)
  Hamming = .54 - .46 * COS(2! * pi * j% / (n% - 1))
END FUNCTION

FUNCTION Hanning (n%, j%)
  Hanning = .5 * (1 - COS(2! * pi * j% / (n% - 1!)))
END FUNCTION

FUNCTION Parzen (n%, j%)
  Parzen = 1 - ABS((j% - .5 * (n% - 1)) / (.5 * (n% - 1)))
END FUNCTION

SUB PowerSpectrumCalc (xreal(), yimag(), numdat%, delta)
DIM numdatr, normal, timespan
DIM i%

   numdatr = numdat% * 1!
   normal = 1! / numdatr ^ 2

   CALL fft(xreal(), yimag(), numdat%, 0)

   xreal(0) = SquareAndSum(xreal(0), yimag(0)) * normal
   FOR i% = 1 TO (numdat% \ 2) - 1
     xreal(i%) = normal * (SquareAndSum(xreal(i%), yimag(i%)) + SquareAndSum(xreal(numdat% - i%), yimag(numdat% - i%)))
   NEXT
   i% = numdat% \ 2
   xreal(i%) = SquareAndSum(xreal(i%), yimag(i%)) * normal
   timespan = numdat% * delta
   FOR i% = 0 TO (numdat% \ 2)
     yimag(i%) = INT(i%) / timespan
   NEXT
END SUB

SUB RealFFT2 (xr(), sinc(), n%, inverse%)
DIM nn%, j%, m%
DIM pp

      m% = -1
      nn% = n% \ 2

      ' determine power of 2
      WHILE (nn% > 0)
	nn% = nn% \ 2
	m% = m% + 1
      WEND
      nn% = n% \ 2

      IF (inverse% = 1) THEN
	pp = .5 / (nn%)
	CALL trans(xr(), sinc(), nn%, 1)
	CALL fft2(xr(), sinc(), nn%, m%, nn%)
	FOR j% = 0 TO nn% - 1
	  xr(j%) = xr(j%) * pp
	  sinc(j%) = sinc(j%) * (-pp)
	NEXT
	CALL reorder(xr(), sinc(), nn%, m%, nn%, 1)
	FOR j% = 0 TO nn% - 1     ' combine
	  xr(j% + nn%) = sinc(j%)
	NEXT
      ELSE
	FOR j% = 0 TO nn% - 1           ' split array in 2
	  sinc(j%) = xr(j% + nn%)
	NEXT
	CALL reorder(xr(), sinc(), nn%, m%, nn%, 1)
	CALL rfft2(xr(), sinc(), nn%, m%, 1)
	FOR j% = 0 TO nn% - 1
	  xr(j%) = xr(j%) * .5
	  sinc(j%) = sinc(j%) * .5
	NEXT
	CALL trans(xr(), sinc(), nn%, 0)
	FOR j% = 0 TO nn% - 1
	  sinc(j%) = -sinc(j%)
	NEXT
	FOR j% = 1 TO nn% - 1
	  xr(j% + nn%) = xr(nn% - j%)
	  sinc(j% + nn%) = -sinc(nn% - j%)
	NEXT
      END IF
END SUB

SUB reorder (xr(), y(), n%, m%, ks%, reel%)

      DIM i%, j%, indx%, k%, kk%, kb%, k2%, KU%, lim%, p%
      DIM breakf%
      DIM temp
      DIM cc%(m%)
      DIM lst(m%)

      cc%(m%) = ks%
      FOR k% = m% TO 1 STEP -1
	cc%(k% - 1) = cc%(k%) \ 2
      NEXT
      p% = m% - 1
      j% = m% - 1
      i% = 0
      kb% = 0
      IF (reel% <> 0) THEN
	KU% = n% - 2
	FOR k% = 0 TO KU% STEP 2
	  temp = xr(k% + 1)
	  xr(k% + 1) = y(k%)
	  y(k%) = temp
	NEXT
      ELSE
	  m% = m% - 1
      END IF
      lim% = (m% + 2) \ 2
      breakf% = 0
      IF (p% > 0) THEN
	WHILE (breakf% = 0)

3000      KU% = cc%(j%) + kb%
	  k2% = KU%
	  indx% = cc%(m% - j%)
	  kk% = kb% + indx%

	  WHILE (kk% < KU%)
	     k% = kk% + indx%
	     WHILE (kk% < k%)
	       temp = xr(kk%)
	       xr(kk%) = xr(k2%)
	       xr(k2%) = temp
	       temp = y(kk%)
	       y(kk%) = y(k2%)
	       y(k2%) = temp
	       kk% = kk% + 1
	       k2% = k2% + 1
	     WEND 'kk% < k%

	     kk% = kk% + indx%
	     k2% = k2% + indx%
	  WEND ' KK < KU

	  IF (j% > lim%) THEN
	    j% = j% - 1
	    i% = i% + 1
	    lst(i%) = j%
	    GOTO 3000
	  END IF
	  kb% = k2%
	  IF (i% > 0) THEN
	    j% = lst(i%)
	    i% = i% - 1
	    GOTO 3000
	  END IF
	  IF (kb% < n%) THEN
	    j% = p%
	  ELSE
	   breakf% = 1
	   GOTO 3010
	  END IF
3010    ' CONTINUE
	WEND ' True

      END IF
END SUB

SUB rfft2 (xr(), y(), n%, m%, ks%)

DIM k0%, k1%, k2%, k3%, k4%, span%, j%, indx%, k%, kb%, nt%, kn%, mk%
DIM Breakf1%, fp%
DIM rad, c1, c2, c3, s1, s2, s3, ck, sk, sq
DIM x0, x1, x2, x3, y0, y1, y2, y3, re, im
DIM cc%(m%)

      sq = .707106781187#
      sk = .382683432366#
      ck = .92387953251#

      cc%(0) = ks%
      kn% = 0
      k4% = 4 * ks%
      mk% = m% - 4

      FOR k% = 1 TO m%
	ks% = ks% + ks%
	cc%(k%) = ks%
      NEXT
      rad = 3.14159265359# / ((cc%(0)) * (ks%))

      WHILE (kn% < n%)
	kb% = kn% + 4
	kn% = kn% + ks%
	IF (m% <> 1) THEN
	  k% = 0
	  indx% = 0
	  j% = mk%
	  nt% = 3
	  c1 = 1!
	  s1 = 0!
	  Breakf1% = 0
	  WHILE (Breakf1% = 0)
2000        span% = cc%(k%)
	    fp% = 0
	    IF (indx% <> 0) THEN
	       c2 = (indx%) * (span%) * rad
	       c1 = COS(c2)
	       s1 = SIN(c2)
	    END IF
2001        IF ((fp% = 1) OR (indx% <> 0)) THEN
	      c2 = c1 * c1 - s1 * s1
	      s2 = 2! * s1 * c1
	      c3 = c2 * c1 - s2 * s1
	      s3 = c2 * s1 + s2 * c1
	    ELSE
	      s1 = 0!
	    END IF
	    k3% = kb% - span%
	    WHILE (k3% < kb%)
	      k2% = k3% - span%
	      k1% = k2% - span%
	      k0% = k1% - span%

	      x0 = xr(k0%)
	      y0 = y(k0%)
	      x1 = xr(k1%)
	      y1 = y(k1%)
	      x2 = xr(k2%)
	      y2 = y(k2%)
	      x3 = xr(k3%)
	      y3 = y(k3%)

	      xr(k0%) = x0 + x2 + x1 + x3
	      y(k0%) = y0 + y2 + y1 + y3
	      IF (s1 = 0!) THEN
		xr(k1%) = x0 - x1 - y2 + y3
		y(k1%) = y0 - y1 + x2 - x3
		xr(k2%) = x0 + x1 - x2 - x3
		y(k2%) = y0 + y1 - y2 - y3
		xr(k3%) = x0 - x1 + y2 - y3
		y(k3%) = y0 - y1 - x2 + x3
	      ELSE
		re = x0 - x1 - y2 + y3
		im = y0 - y1 + x2 - x3
		xr(k1%) = re * c1 - im * s1
		y(k1%) = re * s1 + im * c1
		re = x0 + x1 - x2 - x3
		im = y0 + y1 - y2 - y3
		xr(k2%) = re * c2 - im * s2
		y(k2%) = re * s2 + im * c2
		re = x0 - x1 + y2 - y3
		im = y0 - y1 - x2 + x3
		xr(k3%) = re * c3 - im * s3
		y(k3%) = re * s3 + im * c3
	      END IF 's1 = 0
	      k3% = k3% + 1
	    WEND  'while (k3% < kb%)
	    nt% = nt% - 1
	    IF (nt% >= 0) THEN
	      c2 = c1
	      IF (nt% = 1) THEN
		c1 = c1 * ck + s1 * sk
		s1 = s1 * ck - c2 * sk
	      ELSE
		c1 = (c1 - s1) * sq
		s1 = (s1 + c2) * sq
	      END IF
	      kb% = kb% + k4%
	      IF (kb% <= kn%) THEN
		fp% = 1
		GOTO 2001
	      ELSE
		Breakf1% = 1
		GOTO 2005
	      END IF
	    END IF  'nt% >= 0
	    IF (nt% = -1) THEN
	      k% = 2
	      GOTO 2000
	    END IF
	    IF (cc%(j%) <= indx%) THEN
	      indx% = indx% - cc%(j%)
	      j% = j% - 1
	      IF (cc%(j%) <= indx%) THEN
		indx% = indx% - cc%(j%)
		j% = j% - 1
		k% = k% + 2
	      ELSE
		indx% = cc%(j%) + indx%
		j% = mk%
	      END IF
	    ELSE
	      indx% = cc%(j%) + indx%
	      j% = mk%
	    END IF 'cc%[j%] <= indx%
	    IF (j% < mk%) GOTO 2000
	    k% = 0
	    nt% = 3
	    kb% = kb% + k4%
	    IF (kb% > kn%) THEN Breakf1 = 1
2005        ' CONTINUE
	  WEND  ' while true
	END IF  'm% <> 1
	k% = (m% \ 2) * 2

	IF (k% <> m%) THEN
	  k2% = kn%
	  k0% = kn% - cc%(k%)
	  j% = k0%
	  WHILE (k2% > j%)
	    k2% = k2% - 1
	    k0% = k0% - 1
	    x0 = xr(k2%)
	    y0 = y(k2%)
	    xr(k2%) = xr(k0%) - x0
	    xr(k0%) = xr(k0%) + x0
	    y(k2%) = y(k0%) - y0
	    y(k0%) = y(k0%) + y0
	  WEND  'while (k2% > j%)
	END IF    ' k% <> m%
      WEND         'while (kn% < n%)
END SUB

FUNCTION SquareAndSum (a, b)
  SquareAndSum = a ^ 2 + b ^ 2
END FUNCTION

SUB trans (xr(), y(), n%, eval%)

DIM k%, nk%, nn%
DIM x1, ab, ba, y1, re, im, ck, sk, dc, ds, r

      nn% = n% \ 2
      r = 3.14159265359# / (n%)
      ds = SIN(r)
      r = 2! * SIN(.5 * r)
      r = -r * r
      dc = -.5 * r
      ck = 1!
      sk = 0!

      IF (eval% = 0) THEN
	xr(n%) = xr(0)
	y(n%) = y(0)
      END IF
      FOR k% = 0 TO nn%
	nk% = n% - k%
	x1 = xr(k%) + xr(nk%)
	ab = xr(k%) - xr(nk%)
	ba = y(k%) + y(nk%)
	y1 = y(k%) - y(nk%)

	re = ck * ba + sk * ab
	im = sk * ba - ck * ab
	y(nk%) = im - y1
	y(k%) = im + y1
	xr(nk%) = x1 - re
	xr(k%) = x1 + re
	dc = r * ck + dc
	ck = ck + dc
	ds = r * sk + ds
	sk = sk + ds
      NEXT
END SUB

FUNCTION Welch (n%, j%)
  Welch = 1 - ((j% - .5 * (n% - 1)) / (.5 * (n% + 1))) ^ 2
END FUNCTION

SUB WindowData (x(), numdat%, wind%)
DIM multiplier, i%

  FOR i% = 0 TO numdat% - 1
    SELECT CASE wind%
       CASE 0
	     multiplier = 1
       CASE 1
	     multiplier = Parzen(numdat%, i%)
       CASE 2
	     multiplier = Hanning(numdat%, i%)
       CASE 3
	      multiplier = Welch(numdat%, i%)
       CASE 4
	      multiplier = Hamming(numdat%, i%)
       CASE 5
	      multiplier = ExactBlackman(numdat%, i%)
       CASE ELSE
	      multiplier = 1!
   END SELECT
    x(i%) = multiplier * x(i%)
  NEXT
END SUB

SUB WindowFFTData (xreal(), yimag(), numdat%, wind%)
  CALL WindowData(xreal(), numdat%, wind%)
  CALL WindowData(yimag(), numdat%, wind%)
END SUB