/*  procedures for the Fast Fourier Transform,
  taken from CACM Algorithm 338 by Richard C. Singleton.
  See CACM Vol 11, No 11, Nov 1968, page 773 - 776 for documentation.
  To satisfy the 7 - character limitation for external names, the names
  "COMPLEXTRANSFORM" and "REALTRANSFORM" have been changed to
  "DFTCOMPLEX" and "DFTREAL".  */
#include <math.h>

FFT2 (A, B, N, M, KS)
float   A[],
        B[];
int     N,
        M,
        KS;
{
    int     K0,
            K1,
            K2,
            K3,
            SPAN,
            J,
            JJ,
            K,
            KB,
            KN,
            MM,
            MK;
    float   RAD,
            C1,
            C2,
            C3,
            S1,
            S2,
            S3,
            CK,
            SK,
            SQ;
    float   A0,
            A1,
            A2,
            A3,
            B0,
            B1,
            B2,
            B3;
    int     CC[30];
    SQ = 0.707106781187;
    SK = 0.382683432366;
    CK = 0.92387953251;
    CC[M] = KS;
    MM = (M / 2) * 2;
    KN = 0;
    for (K = M - 1; K >= 0; --K)
	CC[K] = CC[K + 1] / 2;
    RAD = 6.28318530718 / (CC[0] * KS);
    MK = M - 5;
L:  KB = KN;
    KN = KN + KS;
    if (MM != M)
    {
	K2 = KN;
	K0 = CC[MM] + KB;
L2: 	K2 = K2 - 1;
	K0 = K0 - 1;
	A0 = A[K2];
	B0 = B[K2];
	A[K2] = A[K0] - A0;
	A[K0] = A[K0] + A0;
	B[K2] = B[K0] - B0;
	B[K0] = B[K0] + B0;
	if (K0 > KB)
	    goto L2;
    }
    C1 = 1.0;
    S1 = 0.0;
    JJ = 0;
    K = MM - 2;
    J = 3;
    if (K >= 0)
	goto L4;
    else
	goto L6;
L3: if (CC[J] <= JJ)
    {
	JJ = JJ - CC[J];
	J = J - 1;
	if (CC[J] <= JJ)
	{
	    JJ = JJ - CC[J];
	    J = J - 1;
	    K = K + 2;
	    goto L3;
	}
    }
    JJ = CC[J] + JJ;
    J = 3;
L4: SPAN = CC[K];
    if (JJ != 0)
    {
	C2 = JJ * SPAN * RAD;
	C1 = cos ((double) C2);
	S1 = sin ((double) C2);
L5: 	C2 = C1 * C1 - S1 * S1;
	S2 = 2.0 * C1 * S1;
	C3 = C2 * C1 - S2 * S1;
	S3 = C2 * S1 + S2 * C1;
    }
    for (K0 = KB + SPAN - 1; K0 >= KB; --K0)
    {
	K1 = K0 + SPAN;
	K2 = K1 + SPAN;
	K3 = K2 + SPAN;
	A0 = A[K0];
	B0 = B[K0];
	if (S1 == 0.0)
	{
	    A1 = A[K1];
	    B1 = B[K1];
	    A2 = A[K2];
	    B2 = B[K2];
	    A3 = A[K3];
	    B3 = B[K3];
	}
	else
	{
	    A1 = A[K1] * C1 - B[K1] * S1;
	    B1 = A[K1] * S1 + B[K1] * C1;
	    A2 = A[K2] * C2 - B[K2] * S2;
	    B2 = A[K2] * S2 + B[K2] * C2;
	    A3 = A[K3] * C3 - B[K3] * S3;
	    B3 = A[K3] * S3 + B[K3] * C3;
	}
	A[K0] = A0 + A2 + A1 + A3;
	B[K0] = B0 + B2 + B1 + B3;
	A[K1] = A0 + A2 - A1 - A3;
	B[K1] = B0 + B2 - B1 - B3;
	A[K2] = A0 - A2 - B1 + B3;
	B[K2] = B0 - B2 + A1 - A3;
	A[K3] = A0 - A2 + B1 - B3;
	B[K3] = B0 - B2 - A1 + A3;
    }
    if (K > 0)
    {
	K = K - 2;
	goto L4;
    }
    KB = K3 + SPAN;
    if (KB < KN)
    {
	if (J == 0)
	{
	    K = 2;
	    J = MK;
	    goto L3;
	}
	J = J - 1;
	C2 = C1;
	if (J == 1)
	{
	    C1 = C1 * CK + S1 * SK;
	    S1 = S1 * CK - C2 * SK;
	}
	else
	{
	    C1 = (C1 - S1) * SQ;
	    S1 = (C2 + S1) * SQ;
	}
	goto L5;
    }
L6: if (KN < N)
	goto L;

}

REVFFT2 (A, B, N, M, KS)
float   A[],
        B[];
int     N,
        M,
        KS;
{
    int     K0,
            K1,
            K2,
            K3,
            K4,
            SPAN,
            NN,
            J,
            JJ,
            K,
            KB,
            NT,
            KN,
            MK;
    float   RAD,
            C1,
            C2,
            C3,
            S1,
            S2,
            S3,
            CK,
            SK,
            SQ;
    float   A0,
            A1,
            A2,
            A3,
            B0,
            B1,
            B2,
            B3,
            RE,
            IM;
    int     CC[30];

    SQ = 0.707106781187;
    SK = 0.382683432366;
    CK = 0.92387953251;
    CC[0] = KS;
    KN = 0;
    K4 = 4 * KS;
    MK = M - 4;
    for (K = 1; K <= M; ++K)
    {
	CC[K] = KS = KS + KS;
    }
    RAD = 3.14159265359 / (CC[0] * KS);
L:  KB = KN + K4;
    KN = KN + KS;
    if (M <= 1)
	goto L5;
    K = JJ = 0;
    J = MK;
    NT = 3;
    C1 = 1.0;
    S1 = 0.0;
L2: SPAN = CC[K];
    if (JJ != 0)
    {
	C2 = JJ * SPAN * RAD;
	C1 = cos ((double) C2);
	S1 = sin ((double) C2);
L3: 	C2 = C1 * C1 - S1 * S1;
	S2 = 2 * C1 * S1;
	C3 = C2 * C1 - S2 * S1;
	S3 = C2 * S1 + S2 * C1;
    }
    else
	S1 = 0.0;

    K3 = KB - SPAN;
L4: K2 = K3 - SPAN;
    K1 = K2 - SPAN;
    K0 = K1 - SPAN;
    A0 = A[K0];
    B0 = B[K0];
    A1 = A[K1];
    B1 = B[K1];
    A2 = A[K2];
    B2 = B[K2];
    A3 = A[K3];
    B3 = B[K3];
    A[K0] = A0 + A1 + A2 + A3;
    B[K0] = B0 + B1 + B2 + B3;

    if (S1 == 0.0)
    {
	A[K1] = A0 - A1 - B2 + B3;
	B[K1] = B0 - B1 + A2 - A3;
	A[K2] = A0 + A1 - A2 - A3;
	B[K2] = B0 + B1 - B2 - B3;
	A[K3] = A0 - A1 + B2 - B3;
	B[K3] = B0 - B1 - A2 + A3;
    }
    else
    {
	RE = A0 - A1 - B2 + B3;
	IM = B0 - B1 + A2 - A3;
	A[K1] = RE * C1 - IM * S1;
	B[K1] = RE * S1 + IM * C1;
	RE = A0 + A1 - A2 - A3;
	IM = B0 + B1 - B2 - B3;
	A[K2] = RE * C2 - IM * S2;
	B[K2] = RE * S2 + IM * C2;
	RE = A0 - A1 + B2 - B3;
	IM = B0 - B1 - A2 + A3;
	A[K3] = RE * C3 - IM * S3;
	B[K3] = RE * S3 + IM * C3;
    }

    K3 = K3 + 1;
    if (K3 < KB)
	goto L4;
    NT = NT - 1;
    if (NT >= 0)
    {
	C2 = C1;
	if (NT == 1)
	{
	    C1 = C1 * CK + S1 * SK;
	    S1 = S1 * CK - C2 * SK;
	}
	else
	{
	    C1 = (C1 - S1) * SQ;
	    S1 = (C2 + S1) * SQ;
	}
	KB = KB + K4;
	if (KB <= KN)
	    goto L3;
	else
	    goto L5;
    }

    if (NT == -1)
    {
	K = 2;
	goto L2;
    }
    if (CC[J] <= JJ)
    {
	JJ = JJ - CC[J];
	J = J - 1;
	if (CC[J] <= JJ)
	{
	    JJ = JJ - CC[J];
	    J = J - 1;
	    K = K + 2;
	}
	else
	{
	    JJ = CC[J] + JJ;
	    J = MK;
	}
    }
    else
    {
	JJ = CC[J] + JJ;
	J = MK;
    }

    if (J < MK)
	goto L2;
    K = 0;
    NT = 3;
    KB = KB + K4;
    if (KB <= KN)
	goto L2;
L5: K = (M / 2) * 2;

    if (K != M)
    {
	K2 = KN;
	K0 = J = KN - CC[K];
L6: 	K2 = K2 - 1;
	K0 = K0 - 1;
	A0 = A[K2];
	B0 = B[K2];
	A[K2] = A[K0] - A0;
	A[K0] = A[K0] + A0;
	B[K2] = B[K0] - B0;
	B[K0] = B[K0] + B0;
	if (K2 > J)
	    goto L6;
    }
    if (KN < N)
	goto L;

}

REORDER (A, B, N, M, KS, REEL)
float   A[],
        B[];
int     N,
        M,
        KS,
        REEL;
{
    int     I,
            J,
            JJ,
            K,
            KK,
            KB,
            K2,
            KU,
            LIM,
            P;
    float   T;
    int     CC[30],
            LST[30];

    CC[M] = KS;
    for (K = M; K >= 1; --K)
	CC[K - 1] = CC[K] / 2;
    P = J = M - 1;
    I = KB = 0;

    if (REEL)
    {
	KU = N - 2;
	for (K = 0; K <= KU; K = K + 2)
	{
	    T = A[K + 1];
	    A[K + 1] = B[K];
	    B[K] = T;
	}
    }
    else
	M = M - 1;

    LIM = (M + 2) / 2;
    if (P <= 0)
	goto L4;
L:  KU = K2 = CC[J] + KB;
    JJ = CC[M - J];
    KK = KB + JJ;
L2: K = KK + JJ;
L3: T = A[KK];
    A[KK] = A[K2];
    A[K2] = T;
    T = B[KK];
    B[KK] = B[K2];
    B[K2] = T;
    KK = KK + 1;
    K2 = K2 + 1;
    if (KK < K)
	goto L3;
    KK = KK + JJ;
    K2 = K2 + JJ;
    if (KK < KU)
	goto L2;

    if (J > LIM)
    {
	J = J - 1;
	I = I + 1;
	LST[I] = J;
	goto L;
    }

    KB = K2;
    if (I > 0)
    {
	J = LST[I];
	I = I - 1;
	goto L;
    }

    if (KB < N)
    {
	J = P;
	goto L;
    }
L4: ;

}

REALTRAN (A, B, N, EVALUATE)
float   A[],
        B[];
int     N,
        EVALUATE;
{
    int     K,
            NK,
            NH;
    float   AA,
            AB,
            BA,
            BB,
            RE,
            IM,
            CK,
            SK,
            DC,
            DS,
            R;

    NH = N / 2;
    R = 3.14159265359 / N;
    DS = sin ((double) R);
    R = (2 * sin (0.5 * R));
    R = -R * R;
    DC = -0.5 * R;
    CK = 1.0;
    SK = 0.0;

    if (EVALUATE)
    {
	CK = -1.0;
	DC = -DC;
    }
    else
    {
	A[N] = A[0];
	B[N] = B[0];
    }

    for (K = 0; K <= NH; ++K)
    {
	NK = N - K;
	AA = A[K] + A[NK];
	AB = A[K] - A[NK];
	BA = B[K] + B[NK];
	BB = B[K] - B[NK];
	RE = CK * BA + SK * AB;
	IM = SK * BA - CK * AB;
	B[NK] = IM - BB;
	B[K] = IM + BB;
	A[NK] = AA - RE;
	A[K] = AA + RE;
	DC = R * CK + DC;
	CK = CK + DC;
	DS = R * SK + DS;
	SK = SK + DS;
    }
}

DFTREAL (A, B, M, INVERSE)
float   A[],
        B[];
int     M,
        INVERSE;
{
    int     N,
            J;
    float   P;

    N = pow ((double) 2, (double) M);

    if (INVERSE)
    {
	REALTRAN (A, B, N, 1);
	for (J = N - 1; J >= 0; --J)
	    B[J] = -B[J];
	FFT2 (A, B, N, M, N);
	for (J = N - 1; J >= 0; --J)
	{
	    A[J] = 0.5 * A[J];
	    B[J] = -0.5 * B[J];
	}
	REORDER (A, B, N, M, N, 1);
    }
    else
    {
	REORDER (A, B, N, M, N, 1);
	REVFFT2 (A, B, N, M, 1);
	P = 0.5 / N;
	for (J = N - 1; J >= 0; --J)
	{
	    A[J] = P * A[J];
	    B[J] = P * B[J];
	}
	REALTRAN (A, B, N, 0);
    }
    return;
}

DFTCOMPLEX (A, B, M, INVERSE)
float   A[],
        B[];
int     M,
        INVERSE;
{
    int     N,
            J;
    float   P,
            Q;
    N = pow ((double) 2, (double) M);
    P = Q = 1 / sqrt ((double) N);
    if (INVERSE)
    {
	Q = -Q;
	for (J = N - 1; J >= 0; --J)
	    B[J] = -B[J];
    }
    FFT2 (A, B, N, M, N);
    REORDER (A, B, N, M, N, 0);
    for (J = N - 1; J >= 0; --J)
    {
	A[J] = A[J] * P;
	B[J] = B[J] * Q;
    }
    return;
}