unit2.pas

《Delphi常用数值算法集》的配书源码

unit Unit2;

interface
Procedure COSFT(var Y:array of real; N, ISIGN: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 REALFT(var DATA:array of real; N, ISIGN:integer);
var
    I,I1,I2,I3,I4,N2P3:integer;
    WR,WI,C1,C2,THETA,WPR,WPI,WIS,WRS,WRI,H1R,H1I,H2R,H2I,WTEMP:real;
begin
    THETA:=6.28318530717959 / 2 / N;
    C1:=0.5;
    If ISIGN = 1 Then
    begin
        C2:=-0.5;
        FOUR1(DATA, N, 1);
    end
    Else
    begin
        C2:=0.5;
        THETA:=-THETA;
    end; 
    WPR:=-2* Sqr(Sin(0.5 * THETA));
    WPI:=Sin(THETA);
    WR:=1 + WPR;
    WI:=WPI;
    N2P3:= 2 * N + 3;
    For I:=2 To (N div 2) + 1 do
    begin
        I1:=2 * I - 1;
        I2:=I1 + 1;
        I3:=N2P3 - I2;
        I4:=I3 + 1;
        WRS:=WR;
        WIS:=WI;
        H1R:=C1 * (DATA[I1] + DATA[I3]);
        H1I:=C1 * (DATA[I2] - DATA[I4]);
        H2R:=-C2 * (DATA[I2] + DATA[I4]);
        H2I:=C2 * (DATA[I1] - DATA[I3]);
        DATA[I1]:=H1R + WRS * H2R - WIS * H2I;
        DATA[I2]:=H1I + WRS * H2I + WIS * H2R;
        DATA[I3]:=H1R - WRS * H2R + WIS * H2I;
        DATA[I4]:=-H1I + WRS * H2I + WIS * H2R;
        WTEMP:=WR;
        WR:=WR * WPR - WI * WRI + WR;
        WI:=WI * WPR + WTEMP * WPI + WI;
    end;
    If ISIGN = 1 Then
    begin
        H1R:=DATA[1];
        DATA[1]:=H1R + DATA[2];
        DATA[2]:=H1R - DATA[2];
    end
    Else
    begin
        H1R:=DATA[1];
        DATA[1]:=C1 * (H1R + DATA[2]);
        DATA[2]:=C1 * (H1R - DATA[2]);
        FOUR1(DATA, N, -1);
    end; 
end;

Procedure COSFT(var Y:array of real; N, ISIGN:integer);
var
    THETA,WR,WI,WPR,SUM,WTEMP,Y1,Y2,WPI,EVEN,ODD,ENF0,SUM0,SUME:real;
    I,J,M,JJ:integer;
begin
    THETA:=3.14159265358979 / N;
    WR:=1;
    WI:=0;
    WPR:=-2 * Sqr(Sin(0.5 * THETA));
    WPI:=Sin(THETA);
    Sum:=Y[1];
    M:=N div 2;
    For J:=1 To M - 1 do
    begin
        WTEMP:=WR;
        WR:=WR * WPR - WI * WPI + WR;
        WI:=WI * WPR + WTEMP * WPI + WI;
        Y1:=0.5 * (Y[J + 1] + Y[N - J + 1]);
        Y2:=(Y[J + 1] - Y[N - J + 1]);
        Y[J + 1]:=Y1 - WI * Y2;
        Y[N - J + 1]:=Y1 + WI * Y2;
        Sum:=Sum + WR * Y2;
    end; 
    REALFT(Y, M, 1);
    Y[2]:=Sum;
    For JJ:=2 To M  do
    begin
        J:=2*JJ;
        Sum:=Sum + Y[J];
        Y[J]:=Sum;
    end; 
    If ISIGN = -1 Then
    begin
        EVEN:=Y[1];
        ODD:=Y[2];
        For JJ:=1 To M -1 do
        begin
            I:=2*JJ + 1;
            EVEN:=EVEN + Y[I];
            ODD:=ODD + Y[I + 1];
        end;
        ENF0:=2 * (EVEN - ODD);
        SUM0:=Y[1] - ENF0;
        SUME:=(2 * ODD / N) - SUM0;
        Y[1]:=0.5 * ENF0;
        Y[2]:=Y[2] - SUME;
        For JJ:=1 To M -1  do
        begin
            I:=2 * JJ + 1;
            Y[I]:=Y[I] - SUM0;
            Y[I + 1]:=Y[I + 1] - SUME;
        end; 
    end; 
end; 

end.