unit AHoFFT;
(*==========================================================================
Fast Fourier Transform
--------------------------------------------------------------------------
written by : Andreas Hofer
version : 1.1.0, 28.5.2002
============================================================================*)
{$BOOLEVAL OFF }
{$EXTENDEDSYNTAX ON }
{$IOCHECKS ON }
{$LONGSTRINGS ON }
{$OPENSTRINGS ON }
{$TYPEDADDRESS OFF }
{$VARSTRINGCHECKS ON }
{$WRITEABLECONST OFF }
INTERFACE
//=== Basic float type ======================================================
{ The default basic float type for this unit is Single. Use the
conditional define DOUBLEFFTFLOAT below if you want to use
Double as basic float type instead. }
// {$DEFINE DOUBLEFFTFLOAT } // for using double instead of float
{$IFNDEF DOUBLEFFTFLOAT }
type
TFFTFloat = Single;
TFFTDouble = Double; // for intermediate calculations
{$ELSE}
type
TFFTFloat = Double;
TFFTDouble = Extended; // for intermediate calculations
{$ENDIF}
//=== Complex type ==========================================================
// Declaration of the Complex number type and some basic operations on it.
type
TFFTComplex = record
real : TFFTFloat;
imag : TFFTFloat;
end;
// Some basic operations on complex numbers
function Cplx(const re,im : TFFTFloat) : TFFTComplex;
function CplxPolar(const r,phi : TFFTFloat) : TFFTComplex;
function CplxAbs(const z : TFFTComplex) : TFFTFloat;
function CplxArg(const z : TFFTComplex) : TFFTFloat;
function CplxConj(const z : TFFTComplex) : TFFTComplex;
function CplxAdd(const a,b : TFFTComplex) : TFFTComplex;
function CplxSub(const a,b : TFFTComplex) : TFFTComplex;
function CplxMul(const a,b : TFFTComplex) : TFFTComplex;
function CplxQuo(const a,b : TFFTComplex) : TFFTComplex;
function CplxSqrt(const x : TFFTComplex) : TFFTComplex;
function CplxPowerInt(const x : TFFTComplex; exponent : Integer) : TFFTComplex;
//=== FFT Algorithms ========================================================
//--- Basic FFT Algorithms
{
This is the basic FFT algorithm for real and complex signals.
The spectrum array and the signal array must have the same length N.
This algorithm works for any number N but it is fastest if N is a
power of 2. Real signal algorithms also work faster for even N.
}
procedure FFT(const signal : Array of TFFTFloat;
var spectrum : Array of TFFTComplex);
// Real signal FFT
procedure IFFT(const spectrum : Array of TFFTComplex;
var signal : Array of TFFTFloat);
// Real signal inverse FFT
procedure FFTCplx(const signal : Array of TFFTComplex;
var spectrum : Array of TFFTComplex);
// Complex signal FFT
procedure IFFTCplx(const spectrum : Array of TFFTComplex;
var signal : Array of TFFTComplex);
// Complex signal inverse FFT
//--- Fractional FFT
{
The DFT function is based on the definition of the Discrete Fourier
Transform DFT :
DFT(freq) = SUM( k:=0, N-1, x[k]*exp(2*pi*i*k*freq) )
It is recommended only if very few values from the frequency spectrum are
needed, otherwise the FFT or FFFT routines are more efficient.
For freq use teh formula
freq = (desired spectral component) / (sample rate frequency)
}
function DFT(const signal : Array of TFFTFloat; freq : TFFTDouble) : TFFTComplex;
// Real signal DFT
function DFTCplx(const signal : Array of TFFTComplex; freq : TFFTDouble) : TFFTComplex;
// Complex signal DFT
{
FFFT = Fractional Fast Fourier Transformation
This procedure calculates even spaced values of the frequency spectrum of
a signal or a part of the signal. It can be used to calculate only a
small part of the spectrum or to calculate frequency spectrum values
with a different spacing than those from the ordinary FFT (zoom FFT).
'signal' is the signal that is used to calculate the spectrum of.
The number of samples is given by the length of the signal array. Use
the Delphi function Slice if you want only a part of the array to be used.
'tStart' indicates the start time of the signal array in seconds. This
is useful if you want to calculate the spectrum of a part of a signal
which is not at the beginning. This parameter affects the phase of the
resulting spectrum.
'tStep' indicates the sample rate of the signal in seconds.
The 'spectrum' array will receive the calculated spectrum.
The number of spectrum values is given by the length of the spectrum array.
Use the Delphi function Slice if you want only a part of the array to be filled.
'fStart' is the desired starting frequency of the spectrum in Hertz.
'fStep' is the desired frequency spacing of the spectrum in Hertz.
For the inverse FFFT 'spectrum' is the input and 'signal' is output.
'scaleFactor' is applied to the signal array at the end. Usually the
scale factor is calculated as fillows:
scaleFactor = 1 / (total number of spectrum lines) = 1 / Length(signal)
}
procedure FFFT(const signal : Array of TFFTFloat;
tStart, tStep : TFFTDouble;
var spectrum : Array of TFFTComplex;
fStart, fStep : TFFTDouble);
// Real signal FFFT
procedure IFFFT(const spectrum : Array of TFFTComplex;
fStart, fStep : TFFTDouble;
var signal : Array of TFFTFloat;
tStart, tStep : TFFTDouble;
scaleFactor : TFFTFloat);
// Real signal inverse FFFT
procedure FFFTCplx(const signal : Array of TFFTComplex;
tStart, tStep : TFFTDouble;
var spectrum : Array of TFFTComplex;
fStart, fStep : TFFTDouble);
// Complex signal FFFT
procedure IFFFTCplx(const spectrum : Array of TFFTComplex;
fStart, fStep : TFFTDouble;
var signal : Array of TFFTComplex;
tStart, tStep : TFFTDouble;
scaleFactor : TFFTFloat);
// Complex signal inverse FFT
//--- ChirpZ-Transform
{ The Chirp-Z Transform calculates the following, but more efficiently
than using a loop like this:
for f := 0 to High(spectrum) do begin
sum := 0;
for k := 0 to High(signal) do begin
sum := sum + signal[k] * a * b^f * c^k * w^(k*f) // a,b,c,w complex
end;
spectrum[f] := sum;
end;
}
procedure ChirpTransform(const signal : Array of TFFTComplex;
wr,wi, ar,ai, br,bi, cr,ci : Extended;
var spectrum : Array of TFFTComplex);
//=== Weighting Windows =====================================================
{ Weighting windows are used to reduce spurious effects that arise from
sharp transients at the beginning or the end of a signal. The optimal
input for an FFT routine would be either exactly one period of a periodic
signal or a signal that starts and ends with a flat line at zero.
Any other signal will introduce sharp transients (because the FFT
routine is always considering the input as a periodical signal) with
frequency components that may distort the spectrum. To reduce that
problem we can apply a weighting window to the signal first.
}
//--- General weighting windows
type
TFFTWindowFunc = function(x : TFFTFloat) : TFFTFloat;
{ Used to define a window function. x goes from 0 to 1 and the result
of the function should also be in the range 0..1. Used with
'ApplyWindow'. See WTriangle and WWelch for an example of such
a function. }
function WTriangle(x : TFFTFloat) : TFFTFloat;
// Also called Parzen or Bartlett window.
// Goes linearly from 0 to 1 and back to 0.
function WWelch(x : TFFTFloat) : TFFTFloat;
// This window is formed by a parabola which goes from 0 to 1 and back to 0.
type
TFFTWindowInfo = record
intw : TFFTFloat; // integral of window
intw2 : TFFTFloat; // integral of window squared
end;
{ Use this info to undo the amplitude effects of the weighting window in
the spectrum. For example to get the effective values of the sine
waves contained in a signal:
FFT(signal, cplxSpectrum);
RectToAbs(cplxSpectrum, realSpectrum);
ScaleVector(realSpectrum, sqrt(2)/(wi.intw*Length(signal)));
}
function ApplyWindow(var data : array of TFFTFloat;
winFunc : TFFTWindowFunc) : TFFTWindowInfo;
// Apply the window function 'winFunc' to data. Replaces the values in data.
//--- Sum cosine windows
{ Sum cosine windows are weighting windows made up of cosines only }
function ApplySumCosineWindow(var data : array of TFFTFloat;
const coeff : array of Double) : TFFTWindowInfo;
// Apply the sum cosine window function defined by 'coeff' to data.
// Replaces the values in data.
function WSumCosine(const coeff : array of TFFTDouble; x : TFFTFloat) : TFFTFloat;
const // info: first sidelobe level [dB]; sidelobe asymptotic behaviour [db/Oct]
wCoeffHann : array[0..1] of Double = (0.5, 0.5); // 31.5dB; 18dB/Oct
wCoeffHamming : array[0..1] of Double = (0.54, 0.46); // 43dB; 6dB/Oct
wCoeffBlackman : array[0..2] of Double = (0.42, 0.5, 0.08); // 58.11dB; 18dB/Oct
wCoeffNuttall : array[0..3] of Double = (0.355768, 0.487396,
0.144232, 0.012604); // 93.3dB; 18dB/Oct
wCoeffFlatTop2 : array[0..2] of Double = (0.28235,
0.52105, 0.19659); // 43.2dB; 6dB/Oct
wCoeffFlatTop3 : array[0..3] of Double = (0.241096, 0.460841,
0.255381, 0.041872); // 66.8dB; 6dB/Oct
wCoeffFlatTop4 : array[0..4] of Double = (0.209671, 0.407331,
0.281225, 0.092669, 0.009104); // 90.5dB; 18dB/Oct
wCoeffTaylor60 : array[0..11] of Double = (0.46979160648, 0.4897086008,
0.040543436437, 0.00046498453306, -0.00089292984937,
0.00063096192603, -0.00039727029523, 0.00023660292287,
-0.00013197019847, 6.6200331699e-05, -2.7248084479e-05,
7.0249981132e-06); // 60dB; 6dB/Oct
wCoeffTaylor80 : array[0..15] of Double = (0.40987520229, 0.49726668004,
0.091629538934, 0.0011747429958, 2.8285040313e-05,
6.6896112991e-05, -8.0865998033e-05, 7.0758365573e-05,
-5.6300382963e-05, 4.2878214708e-05, -3.1677796476e-05,
2.2686692688e-05, -1.5592228864e-05, 1.005086136e-05,
-5.7617932293e-06, 2.4786542281e-06); // 80dB; 6dB/Oct
procedure TaylorWinCoeff(sideLobedB : TFFTFloat; var coeff : arrayof TFFTDouble);
{ Calculate the taylor window coeffinients for a window with a sideLobe
level of 'sideLobedB'. The number of coefficients is determined by
the length of 'coeff'. }
//=== Auxiliary procedures ==================================================
const// multiply with these to convert from rad to deg and back
radToDeg = 180/pi;
degToRad = pi/180;
function NextPowerOf2(n : Cardinal) : Cardinal;
// if 'n' is a power of 2 return n, otherwise return
// the next greater power of 2.
procedure RealToCplx(const realSignal : Array of TFFTFloat;
var complexSignal : Array of TFFTComplex);
// Convert a real signal into a complex one.
procedure RealToCplx2(const realSignal : Array of TFFTFloat;
const imagSignal : Array of TFFTFloat;
var complexSignal : Array of TFFTComplex);
// Convert two real signals into one complex one.
procedure RealOfCplx(const complexSignal : Array of TFFTComplex;
var realSignal : Array of TFFTFloat);
// Extract the real part of a complex signal.
procedure ImagOfCplx(const complexSignal : Array of TFFTComplex;
var imagSignal : Array of TFFTFloat);
// Extract the imaginary part of a complex signal.
procedure RectToPolar(const spectrum : Array of TFFTComplex;
var amplitude : Array of TFFTFloat;
var phase : Array of TFFTFloat);
// Extract the amplitude and phase of a complex signal.
procedure PolarToRect(const amplitude : Array of TFFTFloat;
const phase : Array of TFFTFloat;
var spectrum : Array of TFFTComplex);
// Convert amplitude and phase into a complex signal.
procedure RectToAbs(const spectrum : Array of TFFTComplex;
var amplitude : Array of TFFTFloat);
// Extract the amplitude of a complex signal.
procedure Reverse(var signal : Array of TFFTFloat);
procedure ReverseCplx(var signal : Array of TFFTComplex);
// Reverses the indices, like 'abcdef' -> 'fedcba'
procedure UnwrapPhase(var phase : Array of TFFTFloat; startIndex : Integer);
// Removes 'Jumps' in the phase spectrum by adding/subtracting 360° when a
// difference > 180° is found between neighbouring samples. The function
// starts at 'startIndex' (this value is unchanged) and works to both sides
// from there.
procedure ScaleVector(var data : Array of TFFTFloat; scalar : TFFTFloat);
procedure ScaleVectorCplx(var data : Array of TFFTComplex; scalar : TFFTComplex);
// Multiplies a vector with a scalar
procedure ZeroPad(const signal : Array of TFFTFloat;
startIndex : Integer;
var result : Array of TFFTFloat);
procedure ZeroPadCplx(const signal : Array of TFFTComplex;
startIndex : Integer;
var result : Array of TFFTComplex);
// 'signal' will be placed into 'result' starting at 'startIndex'. All
// the remaining entries of 'result' before or after the signal are filled
// with zero.
IMPLEMENTATION//************************************************************