{Fourier signal processing library
author: Nikolai V. Shokhirev http://www.shokhirev.com/nikolai.html
created: 2002.02.02
last modified: 2004.12.12
ŠNikolai V. Shokhirev, 2002-2004 }
unit uFourier;
interface
uses
uMatTypes, uDynArrays;
type
TTransformType = (ttDirect, ttInverse);
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<br>
from the "standard" interval:<br>
MinFreqIndex <= M <= MaxFreqIndex <br>
1 <= IndexDFT <= N </pre>}
function IndexDFT(N, M: TInt): TInt;
{<pre> MinFreqIndex = -((N - 1) div 2) </pre>}
function MinFreqIndex(N: TInt): TInt;
{<pre> MaxFreqIndex = (N div 2) = Nyquist Frequency</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: IFArr1D);
{<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: IFArr1D): ICArr1D;
{<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.Dim1 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: ICArr1D; N: TInt): IFArr1D;
{<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: ICArr1D; 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: ICArr1D): ICArr1D;
{<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: ICArr1D): ICArr1D;
{<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: ICArr1D): 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: ICArr1D): ICArr1D;
{<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: ICArr1D): ICArr1D;
{<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: ICArr1D): 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.
Input: Cinp[1..N]
Work: Cscr[1..N] Cscr = Cinp - OK
Output: Cout [1..N] </pre>}
procedure CFTG (const Cinp, Cout, Cscr: ICArr1D; NC: TInt);
{<pre>John Herbster's implementation of Glassman FFT algorithm. <br>
TransformType ttDirect - Fourier transform
Input
Cinp - input signal: Arbitrary limits, Dimension = N
Cout - Initially can have arbitrary limits, Dimension = N
Work
Cscr - work array: Arbitrary limits, Dimension = N, Cscr = Cinp is OK.
Output
Cout[0..N-1] - zero-based array of the Fourier transform.
if Cinp[n] = Cinp(tn) where tn = h*n then the transform frequencies are <br>
0, w, 2*w, 3*w, . . (N-1)*w, where w = 2.0*Pi/(h*N)
TransformType = ttInverse - inverse Fourier transform
Input
Cinp[0..N-1] - zero-based array of the Fourier transform.
Work
Cscr - work array: Arbitrary limits, Dimension = N, Cscr = Cinp is OK.
Output
Cout[Lo1..Hi1] - reconstructed signal: Arbitrary limits, Dimension = N.
FFTG keeps the initial limits. </pre>}
procedure FFTG (const Cinp, Cout, Cscr: ICArr1D; TransformType: TTransformType);
{<pre>
Input
Cinp[0..N-1]
Output
Cout[Mmin..Mmax], Mmin = -((N - 1) div 2), Mmax = (N div 2)
Cout[m] = Cinp[N-m], Mmin <= m <= -1
Cout[m] = Cinp[m], 0 <= m <= Mmax </pre>}
procedure ZeroToCentric(const Cinp, Cout: ICArr1D);
{<pre>
Input
Cinp[Mmin..Mmax], Mmin = -((N - 1) div 2), Mmax = (N div 2)
Output
Cout[0..N-1] </pre>}
procedure CentricToZero(const Cinp, Cout: ICArr1D);
implementation
end.
Back
Generated by Lore's Source to HTML Converter(http://www.newty.de/lsc/index.html)