MODULE BITSpectrum EXPORTS Main;
(* Last edited on 1999-12-22 21:10:14 by hcgl *)

(* Writes the Fourier spectrum (coefs and power) of a uniformly sampled signal. *)

(*
  This program reads a sequence of samples of a real-valued signal f(t).

  It is assumed that the sampling interval "dt" is constant, and the
  signal has period "N*dt" where "N" is a power of two.  The number of
  samples must be "N" or "N+1"; if the latter, the first and last
  samples must be equal.

  The program computes the Fourier transform of the sampled signal
  in the ``smart'' basis

      "basis(N,k)[i] = sqrt(2/N) * cos((2*Pi*k*i / N) + Pi/4)"

  and writes it to disk in two formats: basic coefficients
  "F[k]" (extension ".fft"), and power "P[k]" for each frequency
  "k" (extension ".fpw").  The power is defined as 

      "P[0]   = F[0]^2"
      "P[k]   = F[k]^2 + F[N-k]^2"  for  "k" in "1..N/2-1",
      "P[N/2] = F[N/2]^2"

*)

IMPORT  ParseParams, Process, FileRd, Wr, FileWr;
FROM PZShape IMPORT ComputeSpectrum, ComputePowerSpectrum;
IMPORT  OSError, Thread, Stdio;
(* FROM Math IMPORT sqrt; *)
IMPORT  PZLRChain;

FROM Stdio IMPORT stderr;
FROM PZTypes IMPORT LONG, LONGS;
 
<* FATAL Wr.Failure, Thread.Alerted, OSError.E *>

TYPE
  Options = RECORD
      inFile: TEXT;    (* Input signal file name (WITH extension). *)
      outName: TEXT;   (* Output plot file name (WITHOUT extension). *)
    END;
    
PROCEDURE Main() =
  BEGIN
    WITH
      o = GetOptions(),
      rp = PZLRChain.Read(FileRd.Open(o.inFile), headerOnly := FALSE),
      cmt = "/// original signal\n" & rp.cmt & "\\\\\\ end original signal",
      p = rp.c^,
      a = ComputeSpectrum(p)^,
      unitF = SelectUnit(a),
      pwr = ComputePowerSpectrum(a)^,
      unitP = SelectUnit(pwr)
    DO  
     (* CheckPowerPreservation(p,a); *)
     (* Wr.PutText(stderr, "unitF = " & Fmt.LongReal(unitF) & "\n"); *)
     PZLRChain.Write(
        FileWr.Open(o.outName & ".fft"), "fourier spectrum (smart basis)\n" & cmt, a,
        stride := 0.0d0, unit := unitF
      );
     (* Wr.PutText(stderr, "unitP = " & Fmt.LongReal(unitP) & "\n"); *)
     PZLRChain.Write(
        FileWr.Open(o.outName & ".fpw"), "power spectrum\n" & cmt, pwr,
        stride := 0.0d0, unit := unitP
      )
    END
  END Main;

PROCEDURE SelectUnit(READONLY s: LONGS): LONG =
  CONST MaxScaled = 256.0d0*256.0d0*256.0d0;  (* 2^24 *)
  VAR unit: LONG := 1.0d0/(MaxScaled*MaxScaled*MaxScaled);
  BEGIN
    FOR i := 0 TO LAST(s) DO
      WITH eps = ABS(s[i])/MaxScaled DO
        WHILE unit < eps DO unit := 2.0d0*unit END
      END
    END;
    RETURN unit
  END SelectUnit;


PROCEDURE GetOptions(): Options =
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO
      TRY
        pp.getKeyword("-inFile");
        o.inFile := pp.getNext();
        
        pp.getKeyword("-outName");
        o.outName := pp.getNext();

        pp.finish();      
        
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: BITSpectrum \\\n");
          Wr.PutText(stderr, "  -inFile NAME.EXT \\\n");
          Wr.PutText(stderr, "  -outName NAME\n");
        Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;

(*
PROCEDURE CheckPowerPreservation(READONLY p, a: LONGS) =
  VAR sp, sa: LONG := 0.0d0;
  BEGIN
    FOR i := 0 TO LAST(p) DO sp := sp + p[i]*p[i] END;
    Wr.PutText(stderr, "norm(p) = " & Fmt.LongReal(sqrt(sp)) & "\n");
    FOR i := 0 TO LAST(a) DO sa := sa + a[i]*a[i] END;
    Wr.PutText(stderr, "norm(a) = " & Fmt.LongReal(sqrt(sa)) & "\n");
    Wr.PutText(stderr, "ratio p/a = " & Fmt.LongReal(sqrt(sp/sa)) & "\n");
  END CheckPowerPreservation;
*)
  
BEGIN
  Main()
END BITSpectrum.