Up \Signal

{@abstract(Fourier signal processing library )
author: Nikolai V. Shokhirev )
created: February 02, 2002)
last modified: July 07, 2003)
ŠNikolai V. Shokhirev, 2002-2003 }
unit uFourier;

interface

uses
  uMatTypes, uArrays, uComplexType, uCompoundArrays;

type
  TFactorArray = array [1..32] of TInt;

{<pre> if the input signal is exp(i*2*Pi*M*j/N) <br>
  sampled at j = 0, 1, .. , N-1 then <br>
  the output array from DFT and CFTG is: <br>
 idx |  1 |  2 |  3 | .. | N-1|  N | <br>
  M  |  0 |  1 |  2 | .. | N-2| N-1| <br>
  M  |+/-N|-N+1|-N+2| .. | -2 | -1 | <br>
</pre>}

  {<pre> returns the corresponding index for a given M from the "standard" interval:<br>
    MinFreqIndex <= M <= MaxFreqIndex </pre>}
  function IndexDFT(N, M: TInt): TInt;
{<pre> MinFreqIndex = -((N - 1) div 2) </pre>}
  function MinFreqIndex(N: TInt): TInt;
{<pre> MaxFreqIndex =  (N div 2) </pre>}
  function MaxFreqIndex(N: TInt): TInt;

{<pre> ain input frequency array - output from GFTG >br>
  returns centered frequency array aout <br>
  aout[1] = MinFreq, . . , aout[N] = MaxFreq </pre>}
  procedure CenterFreq(const ain, aout: IVector);

{<pre> DFT - Not fast Fourier transforms -  for testing pirpose only <br>
  see R.W.Hamming, Numerical methods for scientists and engineers, <br>
  McGraw-Hill, N.Y. 1962.</pre>}

{<pre> Forward Discrete Fourier Transform of real data <br>
  It is used in combinatio with InvDFTR and InterpolDFTR<br>
  x[1..N] of TReal - input signal <br>
  result[1..(N+2) div 2] of Complex - Fourier transform <br>
  if x[k] = x(h*(k-1)) then the transform frequencies are <br>
   0, w, 2*w, 3*w, . .  , where w = 2.0*Pi/(h*N)  </pre>}
  function DFTR(const x: IVector): ICVector;
{<pre> Inverse Discrete Fourier Transform to real data <br>
  It is used in combinatio with DFTR <br>
  If c = DFTR(x) then x = InvDFTR(c, N) <br>
  N - Number of real points ( N = x.Dim in DFTR) <br>
  c[1..(N+2) div 2] of Complex - Fourier transform <br>
  result[1..N] of TReal - real signal  </pre>}
  function InvDFTR(const c: ICVector; N: TInt): IVector;
{<pre> Fourier Interpolation of real data
  It is used in combinatio with DFTR <br>
  N - Number of signal points in DFTR <br>
  h - step = distance between signal points <br>
  c[1..(N+2) div 2] of Complex - Fourier transform from DFTR </pre>}
  function InterpolDFTR(t, h: TFloat; const c: ICVector; N: TInt): TFloat;

{<pre> Forward Discrete Fourier Transform of complex data
  It is used in combinatio with InvDFT and InterpolDFT: <br>
  C[1..N] of Complex - input signal <br>
  result[1..N] of Complex - Fourier transform
  if x[k] = x(h*(k-1)) then the transform frequencies are <br>
  0, w, 2*w, 3*w, . .  , where w = 2.0*Pi/(h*N)  </pre>}
  function DFT(const x: ICVector): ICVector;
{<pre> Inverse Discrete Fourier Transform to complex data
  It is used in combinatio with DFT: <br>
  If c = DFT(z) then z = InvDFT(c) <br>
  c[1..N] of Complex - Fourier transform <br>
  result[1..N] of Complex - complex signal  </pre>}
  function InvDFT(const c: ICVector): ICVector;
{<pre> Fourier Interpolation of complex data  <br>
  It is used in combinatio with DFT: <br>
  h - step = distance between signal points <br>
  c[1..N] of Complex - Fourier transform from DFT </pre>}
  function InterpolDFT(t, h: TFloat; const c: ICVector): Complex;

{<pre> Forward DFT for even number of signal points N = 2*M
  It is used in combinatio with InvDFTSym and InterpolSym: <br>
  C[1..N] of Complex - input signal <br>
  result[1..N] of Complex - Fourier transform <br>
  if x[k] = x(h*(k-1)) then the transform frequencies are <br>
  -Pi, -Pi+w, -Pi+2*w, . .  , where w = Pi/(h*M) <br>
  They map to 1, 2, .. , N indices of the result  </pre>}
  function DFTSym(const x: ICVector): ICVector;
{<pre> Inverse DFT for even number of signal points N = 2*M
  It is used in combinatio with DFTSym: <br>
  If c = DFTSym(z) then z = InvDFTSym(c) <br>
  c[1..N] of Complex - Fourier transform <br>
  result[1..N] of Complex - complex signal  </pre>}
  function InvDFTSym(const c: ICVector): ICVector;
{<pre> Fourier Interpolation of complex data with "Symmetric" frequencies <br>
  It is used in combinatio with DFTSym: <br>
  h - step = distance between signal points <br>
  c[1..N] of Complex - Fourier transform from DFT </pre>}
  function InterpolSym(t, h: TFloat; const c: ICVector): Complex;

{<pre> John Herbster's implementation of Glassman FFT algorithm. <br>
  The number of points does not have to be a power of 2. <br>
  Instead, the transform will take advantage of any factoring <br>
  that can be done on the number of points.  Which means that if the number <br>
  of points is prime, the transform will not be "fast" at all. <br>
  <br>
  Returns a complex Fourier transform of the Cinp array into the Cout array. <br>
  The CScr (scratch) array may be same as the input array. The transform will <br>
  be from the time to the frequency domain if NC is positive of the number <br>
  of complex points to be converted or from the frequency to the time domain <br>
  if NC is the negative of the number of complex points.  The only difference <br>
  between the two is just a divide by the number of points. </pre>}
  procedure CFTG (const Cinp, Cout, CScr: ICVector; NC: TInt);

implementation

end.

Up \Signal