realft.pas

来自「Delphi Pascal 数据挖掘领域算法包 数值算法大全」· PAS 代码 · 共 56 行

PROCEDURE realft(VAR data: gldarray; n,isign: integer);
(* Programs using routine REALFT must define the type
TYPE
   gldarray = ARRAY [1..2*n] OF real;
where 2*n is the dimension of the input data array. When
routine FOUR1 is used with REALFT, its data type 'gldarray'
should be set as in this program. *)
VAR
   wr,wi,wpr,wpi,wtemp,theta: double;
   i,i1,i2,i3,i4: integer;
   c1,c2,h1r,h1i,h2r,h2i,wrs,wis: real;
BEGIN
   theta := 6.28318530717959/(2.0*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.0*sqr(sin(0.5*theta));
   wpi := sin(theta);
   wr := 1.0+wpr;
   wi := wpi;
   FOR i := 2 TO (n DIV 2)+1 DO BEGIN
      i1 := i+i-1;
      i2 := i1+1;
      i3 := n+n+3-i2;
      i4 := i3+1;
      wrs := sngl(wr);
      wis := sngl(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*wpi+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;