{*********************************}
{ Generated by c2pas32  v0.9b     }
{ (c) 2001 Oleg Bulychov          }
{ Gladiators Software             }
{ http://www.astonshell.com/      }
{*********************************}

unit fft_c;

interface

(* 
    Fft routines. Originally by J.G.G. Dobbe, modifed to speed up
    by Alexey V. Medvedev 2:5022/78.6
*)
(* PCFFT.C  -- by J.G.G. Dobbe -- Performs an FFT on two arrays (Re, Im) of
               type float (can be changed). This unit is written in C and
               doesn't call assembler routines.
*)
(* --------------------- Local variables -------------------------- *)
(* cos{-2pi/N} for N = 2, 4, 8, ... 16384 *)
(* cos{2pi/N} for N = 2, 4, 8, ... 16384 *)
var {was static}
CosArray: array [0..Pred(28)] of float = (-1.0,0.0,0.70710678118655,0.92387953251129,0.98078528040323,0.99518472667220,0.99879545620517,0.99969881869620,0.99992470183914,0.99998117528260,0.99999529380958,0.99999882345170,0.99999970586288,0.99999992646572,-1.0,0.0,0.70710678118655,0.92387953251129,0.98078528040323,0.99518472667220,0.99879545620517,0.99969881869620,0.99992470183914,0.99998117528260,0.99999529380958,0.99999882345170,0.99999970586288,0.99999992646572); 
(* sin{-2pi/N} for N = 2, 4, 8, ... 16384 *)
(* sin{2pi/N} for N = 2, 4, 8, ... 16384 *)
SinArray: array [0..Pred(28)] of float = (0.0,-1.0,-0.70710678118655,-0.38268343236509,-0.19509032201613,-0.09801714032956,-0.$4906767432742,-0.$A5852852291,-0.$52E6B828572,-0.$C5D88464915,-0.$63795676297,-0.$35B98018628,-0.$1F699031874,-0.$38349518757,0.0,1.0,0.70710678118655,0.38268343236509,0.19509032201613,0.09801714032956,0.$4906767432742,0.$A5852852291,0.$52E6B828572,0.$C5D88464915,0.$63795676297,0.$35B98018628,0.$1F699031874,0.$38349518757); (* --------------------- Function implementations ----------------- *)
(* Another (faster) implementation of Shuffle2Arr                   *)

function SWAP(a: integer; i: integer; j: integer): integer; {<= !!!6 unknown macro}

implementation


procedure shuffle_2_arr(real: pfloat;  imag: pfloat;  n: integer); 
var
temp: float; 
i: integer;
 
j: integer;
 
i0: integer; 
i1: integer; 
i2: integer; 
i3: integer; 
i4: integer; 
i5: integer; 
i6: integer; 
i7: integer; 
i8: integer; 
i9: integer; 
i10: integer; 
i11: integer; 
i12: integer; 
i13: integer; 

function TWICE(k: integer): integer; {<= !!!6 unknown macro}
begin
  
  i:=0; 
  j:=0; 
  
  
  
  
  
  
  
  
  
  
  
  
  
  case n of
    1:
    begin
      TWICE(0); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    2:
    begin
      TWICE(0); 
      TWICE(1); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    3:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    4:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    5:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    6:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    7:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    8:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    9:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    10:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    11:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      TWICE(10); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    12:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      TWICE(10); 
      TWICE(11); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    13:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      TWICE(10); 
      TWICE(11); 
      TWICE(12); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    14:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      TWICE(10); 
      TWICE(11); 
      TWICE(12); 
      TWICE(13); 
      begin 
        if j>i
        then
        begin 
          SWAP(real,i,j); 
          SWAP(imag,i,j); 
        end;
        inc(i); 
      end;
    end;
    else
    begin
    end;
    
  end;{case?}
end;
(* --------------------- Fft -------------------------------------- *)
(* Function     : Actual FFT algorithm. "Re" and "Im" point to start of real 
                 and imaginary arrays of numbers, "Pwr" holds the array sizes
                 as a power of 2 while "Dir" indicates whether an FFT (Dir>=1)
                 or an inverse FFT must be performed (Dir<=0).
  Return value : The transformed information is returned by "Re"
                 and "Im" (real and imaginary part respectively).
  Comment      : modifed to speed up by Alexey V. Medvedev
*)

procedure Fft(Re: pfloat;  Im: pfloat;  Pwr: integer;  Dir: integer); 
var
pr1: pfloat; 
pi1: pfloat; 
pr2: pfloat; 
pi2: pfloat; 
c: float; 
s: float; 
tempr: float; 
tempi: float; 
Qr: float; 
Qi: float; 
Section: integer; 
AngleCounter: integer; 
FlyCount: integer; 
N: integer; (*    Shuffle2Arr(Re, Im, Pwr);    original version: Shuffle before (i)FFT *)
begin
  
  
  
  
  
  
  
  
  
  
  
  
  
  shuffle_2_arr(Re,Im,Pwr); 
  N:= 1 shl Pwr; (* new version, a bit faster *)
  if Dir>=1
  then
  AngleCounter:= 0; 
  else
  (* FFT *)
  AngleCounter:= 14; 
  Section:= 1; (* Inverse FFT *)
  while Section<N
  do
  begin 
    c:= CosArray[AngleCounter]; 
    s:= SinArray[AngleCounter{++} inc(AngleCounter); ]; 
    Qr:= 1; 
    Qi:= 0; 
    for{while} FlyCount:=0 to Pred(Section) { FlyCount++}
    do
    begin 
      pr1:= Re+FlyCount; 
      pi1:= Im+FlyCount; 
      repeat
      begin 
        pr2:= pr1+Section; 
        pi2:= pi1+Section; 
        tempr:= Qr*(*pr2)-Qi*(*pi2); 
        tempi:= Qr*(*pi2)+Qi*(*pr2); {*}pr2^:=*pr1-tempr; 
        {*}pr1^:=*pr1+tempr; 
        {*}pi2^:=*pi1-tempi; 
        {*}pi1^:=*pi1+tempi; 
        pr1:= pr2+Section; 
        pi1:= pi2+Section; 
        (* Perform 2-Point DFT  *)
        (* For Re-part *)
        (* For Im-part *)
      end;
      until not {0=}(pr1-Re<N);
      tempr:= Qr; 
      Qr:= Qr*c-Qi*s; 
      Qi:= Qi*c+tempr*s; 
    end;
    Section:= Section*2; 
  end;
  if Dir<=0
  then
  begin 
    (* Normalize for *)
    (* inverse FFT only *)
    tempr:= 1.0 div N; 
    for{while} Section:=0 to Pred(N) { Section++}
    do
    begin 
      Re[Section]:= Re[Section] * (tempr); 
      Im[Section]:= Im[Section] * (tempr); 
    end;
  end;
end;
(*==============================================================================================*)
(*
Complex version of fft-routines by Alexey V. Medvedev
Uses array of 'struct complex' instead of 2 arrays of float 
-- this is faster
*)
type
complex = record
re: float; 
im: float; 
end;

function comp_SWAP(a: integer; i: integer; j: integer): integer; {<= !!!6 unknown macro}
complex = record
end;


procedure comp_shuffle_2_arr(Cp: pcomplex;  n: integer); 
var
i: integer;
 
j: integer;
 
i0: integer; 
i1: integer; 
i2: integer; 
i3: integer; 
i4: integer; 
i5: integer; 
i6: integer; 
i7: integer; 
i8: integer; 
i9: integer; 
i10: integer; 
i11: integer; 
i12: integer; 
i13: integer; 
begin
  
  i:=0; 
  j:=0; 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  case n of
    1:
    begin
      TWICE(0); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    2:
    begin
      TWICE(0); 
      TWICE(1); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    3:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    4:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    5:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    6:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    7:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    8:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    9:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    10:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    11:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      TWICE(10); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    12:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      TWICE(10); 
      TWICE(11); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    13:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      TWICE(10); 
      TWICE(11); 
      TWICE(12); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    14:
    begin
      TWICE(0); 
      TWICE(1); 
      TWICE(2); 
      TWICE(3); 
      TWICE(4); 
      TWICE(5); 
      TWICE(6); 
      TWICE(7); 
      TWICE(8); 
      TWICE(9); 
      TWICE(10); 
      TWICE(11); 
      TWICE(12); 
      TWICE(13); 
      begin 
        if j>i
        then
        begin 
          comp_SWAP(Cp,i,j); 
        end;
        inc(i); 
      end;
    end;
    else
    begin
    end;
    
  end;{case?}
end;
(* --------------------- Fft -------------------------------------- *)
(* Function     : Actual FFT algorithm. struct complex *Cp point to start of 
                 array of complex numbers, "Pwr" holds the array sizes
                 as a power of 2 while "Dir" indicates whether an FFT (Dir>=1)
                 or an inverse FFT must be performed (Dir<=0).
  Return value : The transformed information is returned by new values
                in Cp.
  Comment      : coded by Alexey V. Medvedev (based on Fft)
*)

procedure comp_Fft(Cp: pcomplex;  Pwr: integer;  Dir: integer); 
var
p1: pcomplex; 
p2: pcomplex; 
c: float; 
s: float; 
tempr: float; 
tempi: float; 
Qr: float; 
Qi: float; 
Section: integer; 
AngleCounter: integer; 
FlyCount: integer; 
N: integer; 
begin
  
  
  
  
  
  
  
  
  
  
  
  
  comp_shuffle_2_arr(Cp,Pwr); 
  N:= 1 shl Pwr; 
  if Dir>=1
  then
  AngleCounter:= 0; 
  else
  (* FFT *)
  AngleCounter:= 14; 
  Section:= 1; (* Inverse FFT *)
  while Section<N
  do
  begin 
    c:= CosArray[AngleCounter]; 
    s:= SinArray[AngleCounter{++} inc(AngleCounter); ]; 
    Qr:= 1; 
    Qi:= 0; 
    for{while} FlyCount:=0 to Pred(Section) { FlyCount++}
    do
    begin 
      p1:= Cp+FlyCount; 
      repeat
      begin 
        p2:= p1+Section; 
        tempr:= Qr*(p2.re)-Qi*(p2.im); 
        tempi:= Qr*(p2.im)+Qi*(p2.re); 
        p2.re:= p1.re-tempr; 
        p1.re:= p1.re+tempr; 
        p2.im:= p1.im-tempi; 
        p1.im:= p1.im+tempi; 
        p1:= p2+Section; 
        (* Perform 2-Point DFT  *)
        (* For Re-part *)
        (* For Im-part *)
      end;
      until not {0=}(p1-Cp<N);
      tempr:= Qr; 
      Qr:= Qr*c-Qi*s; 
      Qi:= Qi*c+tempr*s; 
    end;
    Section:= Section*2; 
  end;
  if Dir<=0
  then
  begin 
    (* Normalize for *)
    (* inverse FFT only *)
    tempr:= 1.0 div N; 
    for{while} Section:=0 to Pred(N) { Section++}
    do
    begin 
      Cp[Section].re:= Cp[Section].re * (tempr); 
      Cp[Section].im:= Cp[Section].im * (tempr); 
    end;
  end;
end;
(* ---------------------------------------------------------------- *)
{$if 1}
(* 
   Very primitive test programm
*)
(* --------------------- Include directives ----------------------- *)
{$include <stdio.h>}
{$include <stdlib.h>}
{$include <unistd.h>}
{$include <time.h>}
(* --------------------- Constant definition ---------------------- *)
const
SIZE = 16384; (* Size of data arrays (re, im) *)
BASE = 14; 
(* --------------------- Variable definitions --------------------- *)
var
re: array [0..Pred(SIZE)] of float; (* Array holding Real part of data *)
im: array [0..Pred(SIZE)] of float; (* Array holding Imaginary part of data *)
type

complex = record
end;(* complex variant *)
(*------------------------------------------------------------------*)


function main(): integer; 
const
NUMC = 20; 
var
i: integer; 
t1: time_t; 
t2: time_t; 
begin
  
  
  
  for{while} i:=0 to Pred(SIZE) { i++}
  do
  begin 
    (* Clear data arrays *)
    sig[i].re:= re[i]:={!!!a type cast? =>} {float(}rand() div {!!!a type cast? =>} {float(}RAND_MAX; 
    sig[i].im:= im[i]:={!!!a type cast? =>} {float(}rand() div {!!!a type cast? =>} {float(}RAND_MAX; 
    (* random values*)
  end;
  sig[0].re:= re[0]:=100.0; (* Fill array with pulse signal *)
  
  printf('Waiting 1 second...'); 
  fflush(stdout); 
  sleep(1); 
  printf(''#13#10'Begin test.'#13#10''); 
  t1:= time(0); 
  for{while} i:=1 to Pred(NUMC) {i++}
  do
  begin 
    Fft(re,im,BASE,1); 
    Fft(re,im,BASE,-1); 
    (* FFT on data points *)
  end;
  (* Inverse FFT on freq data *)
  (*        comp_Fft(sig,BASE,1); *)
  (*        comp_Fft(sig,BASE,-1); *)
  t2:= time(0); 
  printf('End. Time: %g ms'#13#10'',1000. div 2.*(t2-t1) div NUMC); 
  begin
    result:= 0; 
    exit;
  end;
end;
(* ---------------------------------------------------------------- *)
{$endif}

{SWAP(a,i,j) (temp=a[i], a[i]=a[j], a[j]=temp)}
function SWAP(a: integer; i: integer; j: integer): integer; 
begin
  result:= (temp:=a[i],a[i]:=a[j],a[j]:=temp)
end;

{TWICE(k) for(i##k = 2; -- i##k >= 0; j ^= (1<<k) )}
function TWICE(k: integer): integer; 
begin
  result:= for(ik:=2;--ik>=0;j xor =(1 shl k))
end;

{comp_SWAP(a,i,j) (temp=a[i], a[i]=a[j], a[j]=temp)}
function comp_SWAP(a: integer; i: integer; j: integer): integer; 
begin
  result:= (temp:=a[i],a[i]:=a[j],a[j]:=temp)
end;

end.
