#ifndef frb_fourier_H
#define frb_fourier_H

/* Fourier transform operations */
/* Last edited on 2004-12-31 12:23:11 by stolfi */

#include <frb_types.h>

void frb_fourier_transform ( double_vec_t *V, double_vec_t *F );
  /*
    Computes the Fourier transform {F} of a sample vector {V},
    destroying {V} in the process.
    
    This routine actually computes the Hartley transform:

      {V[i] == Sum { F[j] * basis(N,j)(i) : j == 0..N-1 }}

    where

      {basis(N,j)(i) == sqrt(2/N) * cos((2*Pi*j*i / N) + Pi/4)}

    and

      {N == (V.nel) == (F.nel)}.
      
    The {basis} functions are orthonormal, meaning that
    
      {Sum { basis(N,j)(i)*basis(N,k)(i) : i == 0..N-1 } == dirac(j,k)}
      
    for all {j} and {k} in {0..N-1}. Therefore, the Fourier
    transform {F} can be computed by scalar products:
    
      {F[j] == Sum { V[i] * basis(N,j)(i) : i == 0..N-1 }}
      
    The actual frequency of component {F[j]} is 
    {min(j MOD N, (N-j) MOD N)}.  */
  

void frb_power_spectrum ( double_vec_t *F, double_vec_t *P );
  /* 
    Given the Hartley-Fourier transform {F} of a signal,
    stores in {P} its power spectrum. If {F} has {n} elements,
    {P} must have {floor(n/2)+1} elements. {P[0]} will be
    the squared norm of {F[0]}, and {P[n/2]} will be the 
    power at the maximum frequency {n/2}, which comes from
    one or two components of {F}. */
   
unsigned frb_fourier_compute_order ( double band, double step, double length );
  /*
    Selects a suitable number of resampling points for computing the
    FFT of a filtered signal, given the filter cutoff wavelength
    {band}, the minimum spacing {step}, and the total curve length
    {length}. */

void frb_fourier_diff ( double_vec_t *F, unsigned order );
  /*
    Replaces the FFT of a curve by the FFT of its {order}th
    derivative, after discarding the maximum frequency term. */

#endif