#ifndef pz_fourier_H
#define pz_fourier_H

/* Fourier transform operations */
/* Last edited on 2015-01-20 16:49:10 by stolfilocal */

#include <pz_types.h>

void pz_fourier_transform(double_vec_t *V, double_vec_t *F);
  /* Computes the Fourier (Hartley) 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.ne) == (F.ne)}.
      
    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)}.  */

unsigned pz_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 pz_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