unit Unit2;

interface
Procedure FOUR1(var DATA:array of real; NN, ISIGN:integer);
Procedure TWOFFT(DATA1, DATA2:array of real;
                             var FFT1, FFT2:array of real; N:integer);
implementation
Procedure FOUR1(var DATA:array of real; NN, ISIGN:integer);
var
    N,J,JJ,I,M,II,MM,ISTEP,MMAX:integer;
    THETA,WR,WI,WPR,SUM,WTEMP,Y1,Y2,TEMPR,TEMPI,WPI:real;
begin
    N:=2 * NN;
    J:=1;
    For II:=1 To NN do
    begin
        I:= II * 2 - 1;
        If J > I Then
        begin
            TEMPR:=DATA[J];
            TEMPI:=DATA[J + 1];
            DATA[J]:=DATA[I];
            DATA[J + 1]:=DATA[I + 1];
            DATA[I]:=TEMPR;
            DATA[I + 1]:=TEMPI;
        end; 
        M:=N div 2;
        While (M >= 2) And (J > M) do
        begin
            J:=J - M;
            M:=M div 2;
        end;
        J:=J + M;
    end; 
    MMAX:=2;
    While N > MMAX do
    begin
        ISTEP:=2 * MMAX;
        THETA:=6.28318530717959 / (ISIGN * MMAX);
        WPR:=-2 * Sqr(Sin(0.5 * THETA));
        WPI:=Sin(THETA);
        WR:=1;
        WI:=0;
        For II:=1 To (MMAX div 2) do
        begin
            M:= II * 2 - 1;
            For JJ:=0 To ((N-M) div ISTEP) do
            begin
                I:= M + JJ*ISTEP ;
                J:=I + MMAX;
                TEMPR:=WR * DATA[J] - WI * DATA[J + 1];
                TEMPI:=WR * DATA[J + 1] + WI * DATA[J];
                DATA[J]:=DATA[I] - TEMPR;
                DATA[J + 1]:=DATA[I + 1] - TEMPI;
                DATA[I]:=DATA[I] + TEMPR;
                DATA[I + 1]:=DATA[I + 1] + TEMPI;
            end; 
            WTEMP:=WR;
            WR:=WR * WPR - WI * WPI + WR;
            WI:=WI * WPR + WTEMP * WPI + WI;
        end; 
        MMAX:=ISTEP;
    end;
end;

Procedure TWOFFT(DATA1, DATA2:array of real;
                             var FFT1, FFT2:array of real; N:integer);
var
    C1R,C1I,C2R,C2I,CONJR,CONJI,H1R,H1I,H2R,H2I:real;  N2,J,J2:integer;
begin
    C1R:=0.5;
    C1I:=0;
    C2R:=0;
    C2I:=-0.5;
    For J:=1 To N do
    begin
        FFT1[2 * J - 1]:=DATA1[J];
        FFT1[2 * J]:=DATA2[J];
    end;
    FOUR1(FFT1, N, 1);
    FFT2[1]:=FFT1[2];
    FFT2[2]:=0;
    FFT1[2]:=0;
    N2:=2 * (N + 2);
    For J:=2 To ((N div 2) + 1) do
    begin
        J2:=2 * J;
        CONJR:=FFT1[N2 - J2 - 1];
        CONJI:=-FFT1[N2 - J2];
        H1R:=C1R * (FFT1[J2 - 1] + CONJR) - C1I * (FFT1[J2] + CONJI);
        H1I:=C1I * (FFT1[J2 - 1] + CONJR) + C1R * (FFT1[J2] + CONJI);
        H2R:=C2R * (FFT1[J2 - 1] - CONJR) - C2I * (FFT1[J2] - CONJI);
        H2I:=C2I * (FFT1[J2 - 1] - CONJR) + C2R * (FFT1[J2] - CONJI);
        FFT1[J2 - 1]:=H1R;
        FFT1[J2]:=H1I;
        FFT1[N2 - J2 - 1]:=H1R;
        FFT1[N2 - J2]:=-H1I;
        FFT2[J2 - 1]:=H2R;
        FFT2[J2]:=H2I;
        FFT2[N2 - J2 - 1]:=H2R;
        FFT2[N2 - J2]:=-H2I;
    end;
end; 

end.