#ifndef pz_shape_H
#define pz_shape_H

/* Converts a curve segment into a "shape function" for Fourier analysis. */
/* Last edited on 2008-02-08 12:11:55 by stolfi */

#include <pz_r3_chain.h>
#include <pz_double_chain.h>
#include <pz_types.h>

pz_double_chain_t *pz_shape_get(pz_r3_chain_t *c);
  /* !!! Replace by the doubly-integrated curvature (see {FRB}) !!!
    
    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... */

typedef double_vec_t pz_power_spectrum_t; 
  /* Power spectrum of a signal, indexed by (absolute) frequency.
    Specifically, if the signal has {n} samples, the spectrum
    has {floor(n/2)} elements; element {i} is the
    power in frequencies {i} and {-i}. */

typedef double_vec_t pz_fourier_t; 
  /* Discrete Fourier (Hartley) transform of a signal, indexed by
    (signed) frequency modulo the number of samples {n}. Specifically,
    if the signal has {n} samples, its transform has {n} elements, and
    element {i} is the amplitude of the Hartley basis function with
    frequency {i}. */

pz_fourier_t *pz_fourier_compute(pz_double_chain_t *p);
  /* Computes the Fourier transform of the shape function "p",
    in the Hartley 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, sample {[0]} and 
    sample {[N]} must be equal.  In either case the period is assumed to be
    "N*dt"

    The Fourier transform uses the Hartley basis

        "basis(N,k)[i] == sqrt(2/N) * cos((2*Pi*k*i / N) + Pi/4)"
        
    which is orthonormal under the scalar product 
    {<f|g> = SUM{ f[i]*g[i] : i IN 0..N-1}. 
  */
  

pz_power_spectrum_t *pz_power_spectrum_compute(pz_fourier_t *F)
  /* Returns a vector "pwr[f][k]" with the power (square of amplitude)
    for each coordinate "k" and each frequency "f", from 0 to "fMax == floor(N/2)".
    Element "pwr[0]" includes only one term of the series (the constant term).
    If {N} is even, element "pwr[fMax]" also includes only one element
    (the Nyquist limit frequency {N/2}, which had better be absent from the signal). 
    All other output elements include two components with 
    the same absolute frequency.  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"

  */

#endif