INTERFACE PZShape;
(* Last edited on 1999-08-07 00:40:56 by hcgl *)

(* Converts a curve segment into a "shape function" for Fourier analysis. *)

IMPORT PZLR3Chain, PZLRChain;
FROM PZTypes IMPORT LONGS;

PROCEDURE Signal(READONLY c: PZLR3Chain.T): REF PZLRChain.T;
  (*
    Converts an open 2D curve into a sequence of samples that can be
    used as a periodic "signal" in the Hartley-Shannon theorem.

    The curve must be given as "N+1" samples "c[0..N]" uniformly spaced
    along the curve, where "N" is a power of two. The output signal is
    computed by the following algorithm:

      1. If "n=0" return the sequence "(0,0)".

      2. Divide the curve into two curves "a[0..N/2]", "b[0..N/2]",
         each with with "N/2" steps and "N/2+1" points,
         at a the middle sample point "r = c[N/2]". 

      3. Convert the two curves recursively into signals "A[0..N/2]", 
         "B[0..N/2]" each with "N/2+1" samples.

      4. Let "p=c[0]" and "q=c[N]" be the endpoints of the curve, and
         "theta" be the angle between the vectors "r-p" and "q-r",
         reasured CLOCKWISE.

      5. Let "C[0..N]" be a tent sequence with peak value 
         "h = L*theta/4" in the center sample and 0 at the extremities,
         where "L" is the total curve length.  Add the "A" sequence to 
         "C[0..N/2]", and the "B" seqence to the "C[N/2..N]", both scaled
         by 1/2.

    In step 4 the angle "theta" should include the correct number
    of full turns.  The current implementation doesn't do that...

  *)

TYPE
  PowerSpectrum = LONGS;  (* Indexed by frequency *)
  Spectrum = LONGS;       (* Indexed by coefs *)
  
PROCEDURE ComputeSpectrum(READONLY p: PZLRChain.T): REF Spectrum;
  (*
    Computes the Fourier spectrum of the shape function "p",
    in the ``smart'' Fourier basis.
    
    Assumes the sampling interval "dt" is constant, and "p" has either
    "N=2^n" or "N+1=2^n+1" samples; if the latter, the first and last
    sample must be equal.  In either case the period is assumed to be
    "N*dt"

    The Fourier transform uses the ``smart'' basis

        "basis(N,k)[i] = sqrt(2/N) * cos((2*Pi*k*i / N) + Pi/4)"
  *)
  
PROCEDURE ComputePowerSpectrum(READONLY a: Spectrum): REF PowerSpectrum;
  (*
    Returns a vector "pwr[f][k]" with the power (square of amplitude)
    for each coordinate "k" and each frequency "f", from 0 to "fMax = N DIV 2".
    Elements "pwr[0]" and "pwr[fMax]" include only one term each
    of the series; all other output elements include two conjugate
    terms.  More precisely, the power spectrum 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"

  *)

END PZShape.
(*
  Copyright © 2001 Universidade Estadual de Campinas (UNICAMP).
  Authors: Helena C. G. Leitão and Jorge Stolfi.
  
  This file can be freely distributed, used, and modified, provided
  that this copyright and authorship notice is preserved, and that any
  modified versions are clearly marked as such.
  
  This software has NO WARRANTY of correctness or applicability for
  any purpose. Neither the authors nor their employers chall be held
  responsible for any losses or damages that may result from its use.
*)