with Ada.Text_IO;            use Ada.Text_IO;
with Ada.Text_IO.C_Streams;
with Ada.Long_Float_Text_IO; use Ada.Long_Float_Text_IO;
with Ada.Integer_Text_IO;    use Ada.Integer_Text_IO;
with Interfaces.C;           use Interfaces.C;
with Interfaces.C.Strings;   use Interfaces.C.Strings;
with Ada.Strings.Unbounded;  use Ada.Strings.Unbounded;
with Ada.Command_Line;       use Ada.Command_Line;
with Ada.Numerics.Elementary_Functions ; use Ada.Numerics.Elementary_Functions ;

with GNAT.Source_Info; use GNAT.Source_Info;
with gsl ;
with gsl.complex ;
with gsl.complex.math; 
with gsl.fft.complex ;
with gsl.fft.real ;
with gsl.fft.halfcomplex ;

with ehandler ;

procedure Fft is
   use gsl.double_text_io;
   
   verbose : boolean := true;
   Status : Int ;
   procedure Sep is
   begin
      Put(" ; ");
   end Sep ;
   procedure Complex_Print( c : gsl.complex.gsl_complex ) is
   begin
      Put( gsl.complex.Re(c) ); Sep ;
      Put( gsl.complex.Im(c) ); Sep ;  
      Put( gsl.complex.math.cabs(c)); Sep ;
      Put( gsl.complex.math.arg(c)); Sep;    
   end Complex_Print ;

   procedure T1 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 128 ;
      NSqrt : constant double := double(Sqrt( float(N) )) ;
      data : gsl.complex.gsl_complex_complex_array(0..N-1) := (others => (0.0,0.0) );
      temp : gsl.complex.gsl_complex ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      data(0) := (1.0,0.0) ;
      for i in 1..10
      loop
         data(i) := (1.0,0.0) ;
         data(N-i) := (1.0,0.0) ;
      end loop ;

      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      for i in data'range
      loop
         Put(i) ; Sep ;
         Complex_Print(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      Status := gsl.fft.complex.radix2_forward( data , 1 , data'Length );

      Create(logfile,Out_File,myname & ".fft.csv");
      Set_Output(logfile);
      for i in data'range
      loop
         temp := gsl.complex.math.div_real(data(i) , NSqrt );
         data(i) := temp ;
         Put(i) ; Sep ;
         Complex_Print(temp);
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      Status := gsl.fft.complex.radix2_backward( data , 1 , data'Length );

      Create(logfile,Out_File,myname & ".rfft.csv");
      Set_Output(logfile);      
      for i in data'range
      loop
         Put(i) ; Sep ;
         temp := gsl.complex.math.div_real(data(i) , NSqrt );
         data(i) := temp ;
         Complex_Print(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

   end T1 ;
   procedure T2 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 630 ;
      NSqrt : constant double := double(Sqrt( float(N) )) ;
      data : gsl.complex.gsl_complex_complex_array(0..N-1) := (others => (0.0,0.0) );
      temp : gsl.complex.gsl_complex ;
      wt : access gsl.fft.complex.gsl_fft_complex_wavetable;
      ws : access gsl.fft.complex.gsl_fft_complex_workspace;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      data(0) := (1.0,0.0) ;
      for i in 1..10
      loop
         data(i) := (1.0,0.0) ;
         data(N-i) := (1.0,0.0) ;
      end loop ;
      wt := gsl.fft.complex.wavetable_alloc(N);
      ws := gsl.fft.complex.workspace_alloc(N);
      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      for i in data'range
      loop
         Put(i) ; Sep ;
         Complex_Print(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      Status := gsl.fft.complex.forward( data , 1 , data'Length , wt , ws );

      Create(logfile,Out_File,myname & ".fft.csv");
      Set_Output(logfile);
      for i in data'range
      loop
         temp := gsl.complex.math.div_real(data(i) , NSqrt );
         data(i) := temp ;
         Put(i) ; Sep ;
         Complex_Print(temp);
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      Status := gsl.fft.complex.backward( data , 1 , data'Length , wt , ws );

      Create(logfile,Out_File,myname & ".rfft.csv");
      Set_Output(logfile);      
      for i in data'range
      loop
         Put(i) ; Sep ;
         temp := gsl.complex.math.div_real(data(i) , NSqrt );
         data(i) := temp ;
         Complex_Print(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      gsl.fft.complex.workspace_free(Ws);
      gsl.fft.complex.wavetable_free(wt);
   end T2 ;

   procedure T3 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 100 ;
      NSqrt : constant double := double(Sqrt( float(N) )) ;
      data : gsl.complex.gsl_complex_complex_array(0..N-1) := (others => (0.0,0.0) );
      temp : gsl.complex.gsl_complex ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      data(0) := (1.0,0.0) ;
      for i in 1..10
      loop
         data(i) := (1.0,0.0) ;
         data(N-i) := (1.0,0.0) ;
      end loop ;
      ehandler.Setup ;
      
      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      for i in data'range
      loop
         Put(i) ; Sep ;
         Complex_Print(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      begin
         Status := gsl.fft.complex.radix2_forward( data , 1 , data'Length );
      exception
         when others =>
               Put("Exception in "); Put(myname); Put(" exiting..."); New_Line ;
               ehandler.Reset ;
               return ;
      end ;
      Create(logfile,Out_File,myname & ".fft.csv");
      Set_Output(logfile);
      for i in data'range
      loop
         temp := gsl.complex.math.div_real(data(i) , NSqrt );
         data(i) := temp ;
         Put(i) ; Sep ;
         Complex_Print(temp);
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      Status := gsl.fft.complex.radix2_backward( data , 1 , data'Length );

      Create(logfile,Out_File,myname & ".rfft.csv");
      Set_Output(logfile);      
      for i in data'range
      loop
         Put(i) ; Sep ;
         temp := gsl.complex.math.div_real(data(i) , NSqrt );
         data(i) := temp ;
         Complex_Print(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

   end T3 ;

   procedure T4 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 128 ;
      --NSqrt : constant double := double(Sqrt( float(N) )) ;
      data : gsl.double_array(0..N-1) := (others => 0.0);

      ws : access gsl.fft.real.gsl_fft_real_workspace;
      wt : access gsl.fft.real.gsl_fft_real_wavetable;
      wth : access gsl.fft.halfcomplex.gsl_fft_halfcomplex_wavetable;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;

      for i in N/3 .. 2 * N/3
      loop
         data(i) := 1.0 ;
      end loop ;

      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      for i in data'range
      loop
         Put(i) ; Sep ;
         Put(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      ws := gsl.fft.real.workspace_alloc(N);
      wt := gsl.fft.real.wavetable_alloc(N);
      wth := gsl.fft.halfcomplex.wavetable_alloc(N);

      Status := gsl.fft.real.transform(data,1,N,wt,ws);
      -- Low pass filter
      for i in 11..data'Last
      loop
         data(i) := 0.0 ;
      end loop ;
      Status := gsl.fft.halfcomplex.inverse(data,1,N,wth,ws);

      Create(logfile,Out_File,myname & ".filtered.csv");
      Set_Output(logfile);
      for i in data'range
      loop
         Put(i) ; Sep ;
         Put(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      gsl.fft.halfcomplex.wavetable_free(wth);
      gsl.fft.real.wavetable_free(wt);
      gsl.fft.real.workspace_free(ws);

   end T4 ;
   procedure T5 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 256 ;
      --NSqrt : constant double := double(Sqrt( float(N) )) ;
      data : gsl.double_array(0..N-1) := (others => 0.0);
      freqs : gsl.complex.gsl_complex_complex_array(0..N-1) := (others => (0.0,0.0)) ;
      t : double := 0.0 ;
      tdelta : double := 1.0/double(N) ;
      ws : access gsl.fft.real.gsl_fft_real_workspace;
      wt : access gsl.fft.real.gsl_fft_real_wavetable;
      wth : access gsl.fft.halfcomplex.gsl_fft_halfcomplex_wavetable;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;

      for i in data'Range 
      loop
         data(i) := double(Sin(Float(Ada.Numerics.Pi * 2.0 * t))) + 
                    2.0 * double(Sin(Float(Ada.Numerics.Pi * 8.0 * t))) ;
         t := t + tdelta ;
      end loop ;

      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      for i in data'range
      loop
         Put(i) ; Sep ;
         Put(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      ws := gsl.fft.real.workspace_alloc(N);
      wt := gsl.fft.real.wavetable_alloc(N);
      wth := gsl.fft.halfcomplex.wavetable_alloc(N);

      Status := gsl.fft.real.transform(data,1,N,wt,ws);
      Status := gsl.fft.halfcomplex.unpack(data,freqs,1,N);
      Create(logfile,Out_File,myname & ".freq.csv");
      Set_Output(logfile);
      for i in 1..16
      loop
         Put(i) ; Sep ;
         Put(data(i)); Sep ;
         Complex_Print(freqs(i));
         Put(gsl.complex.math.cabs(freqs(i))*2.0/double(N)) ;
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      Status := gsl.fft.halfcomplex.inverse(data,1,N,wth,ws);

      Create(logfile,Out_File,myname & ".filtered.csv");
      Set_Output(logfile);
      for i in data'range
      loop
         Put(i) ; Sep ;
         Put(data(i));
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      gsl.fft.halfcomplex.wavetable_free(wth);
      gsl.fft.real.wavetable_free(wt);
      gsl.fft.real.workspace_free(ws);

   end T5 ;
begin

   T1 ;
   T2 ;
   T3 ;
   T4 ;
   T5 ;

end Fft;