Anyone got a Fourier (or suchlike) routine for getting frequencies & amplitudes of music? I have a grandkids’ birthday coming up!
1 Like
{ Unit FFTs
This unit provides a forward and inverse FFT pascal implementation
for complex number series.
The formal definition of the complex DFT is:
y[k] = sum(x[m]*exp(-i*2*pi*k*m/n), m = 0..n-1), k = 0..n-1
Copyright: Nils Haeck M.Sc. (email: n.haeck@simdesign.nl)
For more information visit http://www.simdesign.nl
Original date of publication: 10 Mar 2003
This unit requires these other units:
- Complexs: Complex number unit
- Types: Additional mathematical variable types
- SysUtils: Delphi system utilities
****************************************************************
The contents of this file are subject to the Mozilla Public
License Version 1.1 (the "License"); you may not use this file
except in compliance with the License. You may obtain a copy of
the License at:
http://www.mozilla.org/MPL/
Software distributed under the License is distributed on an
"AS IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or
implied. See the License for the specific language governing
rights and limitations under the License.
}
unit FFTs;
interface
uses
Complexs, Types, SysUtils, Dialogs;
const
cMaxFactorCount = 20;
cMaxPrimeFactor = 1021;
cMaxPrimeFactorDiv2 = (cMaxPrimeFactor + 1) div 2;
resourcestring
sErrPrimeTooLarge = 'Prime factor for FFT length too large. Change value for cMaxPrimeFactor in FFTs unit';
{ ForwardFFT:
Perform a complex FFT on the data in Source, put result in Dest. This routine
works best for Count as a power of 2, but also works usually faster than DFT
by factoring the series. Only in cases where Count is a prime number will this
method be identical to regular complex DFT.
The largest prime factor in Count should be less or equal to cMaxPrimeFactor.
The remaining factors are handled by optimised partial FFT code, that can be
found in the FFT_X procedures
Inputs:
Source: this can be any zero-based array type of TComplex
Count: The number of elements in the array.
Outputs:
Dest: this can be any zero-based array type of TComplex, and will contain
the FFT transformed data (frequency spectrum). Source may be equal to
Dest. In this case, the original series will be overwritten with the new
fourier-transformed series.
}
procedure ForwardFFT(const Source: array of TComplex; var Dest: array of TComplex; Count: integer);
{ Perform the inverse FFT on the Source data, and put result in Dest. This is based
on the forward FFT with some additional customisation. The result of a forward
FFT followed by an inverse FFT should yield the same data, except for rounding
errors.
}
procedure InverseFFT(const Source: array of TComplex; var Dest: array of TComplex; Count: integer);
implementation
const
// Some helper constants for the FFT optimisations
c31: TFloat = -1.5000000000000E+00; // cos(2*pi / 3) - 1;
c32: TFloat = 8.6602540378444E-01; // sin(2*pi / 3);
u5: TFloat = 1.2566370614359E+00; // 2*pi / 5;
c51: TFloat = -1.2500000000000E+00; // (cos(u5) + cos(2*u5))/2 - 1;
c52: TFloat = 5.5901699437495E-01; // (cos(u5) - cos(2*u5))/2;
c53: TFloat = -9.5105651629515E-01; //- sin(u5);
c54: TFloat = -1.5388417685876E+00; //-(sin(u5) + sin(2*u5));
c55: TFloat = 3.6327126400268E-01; // (sin(u5) - sin(2*u5));
c8: TFloat = 7.0710678118655E-01; // 1 / sqrt(2);
type
// Base 1 and Base 0 arrays
TIdx0FactorArray = array[0..cMaxFactorCount] of integer;
TIdx1FactorArray = array[1..cMaxFactorCount] of integer;
// Factorise the series with length Count into FactorCount factors, stored in Fact
procedure Factorize(Count: integer; var FactorCount: integer; var Fact: TIdx1FactorArray);
var
i: integer;
Factors: TIdx1FactorArray;
const
// Define specific FFT lengths (radices) that we can process with optimised routines
cRadixCount = 6;
cRadices: array[1..6] of integer =
(2, 3, 4, 5, 8, 10);
begin
if Count = 1 then begin
FactorCount := 1;
Factors[1] := 1;
end else begin
FactorCount := 0;
end;
// Factorise the original series length Count into known factors and rest value
i := cRadixCount;
while (Count > 1) AND (i > 0) do begin
if Count mod cRadices[i] = 0 then begin
Count := Count div cRadices[i];
inc(FactorCount);
Factors[FactorCount] := cRadices[i];
end else
dec(i);
end;
// substitute factors 2*8 with more optimal 4*4
if Factors[FactorCount] = 2 then begin
i := FactorCount - 1;
while (i > 0) AND (Factors[i] <> 8) do
dec(i);
if i > 0 then begin
Factors[FactorCount] := 4;
Factors[i] := 4;
end;
end;
// Analyse the rest value and see if it can be factored in primes
if Count > 1 then begin
for i := 2 to trunc(sqrt(Count)) do begin
while Count mod i = 0 do begin
Count := Count div i;
inc(FactorCount);
Factors[FactorCount] := i;
end;
end;
if (Count > 1) then begin
inc(FactorCount);
Factors[FactorCount] := Count;
end;
end;
// Reverse factors so that primes are first
for i := 1 to FactorCount do
Fact[i] := Factors[FactorCount - i + 1];
end;
{ Reorder the series in X to a permuted sequence in Y so that the later step can
be done in place, and the final FFT result is in correct order.
The series X and Y must be different series!
}
procedure ReorderSeries(Count: integer; var Factors: TIdx1FactorArray; var Remain: TIdx0FactorArray;
const X: array of TComplex; var Y: array of TComplex);
var
i, j, k: integer;
Counts: TIdx1FactorArray;
begin
FillChar(Counts, SizeOf(Counts), 0);
k := 0;
for i := 0 to Count - 2 do begin
Y[i] := X[k];
j := 1;
k := k + Remain[j];
Counts[1] := Counts[1] + 1;
while Counts[j] >= Factors[j] do begin
Counts[j] := 0;
k := k - Remain[j - 1] + Remain[j + 1];
inc(j);
inc(Counts[j]);
end;
end;
Y[Count - 1] := X[Count - 1];
end;
procedure FFT_2(var Z: array of TComplex);
var
T1: TComplex;
begin
T1 := ComplexAdd(Z[0], Z[1]);
Z[1] := ComplexSub(Z[0], Z[1]);
Z[0] := T1;
end;
procedure FFT_3(var Z: array of TComplex);
var
T1, M1, M2, S1: TComplex;
begin
T1 := ComplexAdd(Z[1], Z[2]);
Z[0] := ComplexAdd(Z[0], T1);
M1 := ComplexScl(c31, T1);
M2.Re := c32 * (Z[1].Im - Z[2].Im);
M2.Im := c32 * (Z[2].Re - Z[1].Re);
S1 := ComplexAdd(Z[0], M1);
Z[1] := ComplexAdd(S1, M2);
Z[2] := ComplexSub(S1, M2);
end;
procedure FFT_4(var Z: array of TComplex);
var
T1, T2, M2, M3: TComplex;
begin
T1 := ComplexAdd(Z[0], Z[2]);
T2 := ComplexAdd(Z[1], Z[3]);
M2 := ComplexSub(Z[0], Z[2]);
M3.Re := Z[1].Im - Z[3].Im;
M3.Im := Z[3].Re - Z[1].Re;
Z[0] := ComplexAdd(T1, T2);
Z[2] := ComplexSub(T1, T2);
Z[1] := ComplexAdd(M2, M3);
Z[3] := ComplexSub(M2, M3);
end;
procedure FFT_5(var Z: array of TComplex);
var
T1, T2, T3, T4, T5: TComplex;
M1, M2, M3, M4, M5: TComplex;
S1, S2, S3, S4, S5: TComplex;
begin
T1 := ComplexAdd(Z[1], Z[4]);
T2 := ComplexAdd(Z[2], Z[3]);
T3 := ComplexSub(Z[1], Z[4]);
T4 := ComplexSub(Z[3], Z[2]);
T5 := ComplexAdd(T1, T2);
Z[0] := ComplexAdd(Z[0], T5);
M1 := ComplexScl(c51, T5);
M2 := ComplexScl(c52, ComplexSub(T1, T2));
M3.Re := -c53 * (T3.Im + T4.Im);
M3.Im := c53 * (T3.Re + T4.Re);
M4.Re := -c54 * T4.Im;
M4.Im := c54 * T4.Re;
M5.Re := -c55 * T3.Im;
M5.Im := c55 * T3.Re;
S3 := ComplexSub(M3, M4);
S5 := ComplexAdd(M3, M5);;
S1 := ComplexAdd(Z[0], M1);
S2 := ComplexAdd(S1, M2);
S4 := ComplexSub(S1, M2);
Z[1] := ComplexAdd(S2, S3);
Z[2] := ComplexAdd(S4, S5);
Z[3] := ComplexSub(S4, S5);
Z[4] := ComplexSub(S2, S3);
end;
procedure FFT_8(var Z: array of TComplex);
var
A, B: array[0..3] of TComplex;
Gem: TFloat;
begin
A[0] := Z[0]; B[0] := Z[1];
A[1] := Z[2]; B[1] := Z[3];
A[2] := Z[4]; B[2] := Z[5];
A[3] := Z[6]; B[3] := Z[7];
FFT_4(A);
FFT_4(B);
Gem := c8 * (B[1].Re + B[1].Im);
B[1].Im := c8 * (B[1].Im - B[1].Re);
B[1].Re := Gem;
Gem := B[2].Im;
B[2].Im :=-B[2].Re;
B[2].Re := Gem;
Gem := c8 * (B[3].Im - B[3].Re);
B[3].Im :=-c8 * (B[3].Re + B[3].Im);
B[3].Re := Gem;
Z[0] := ComplexAdd(A[0], B[0]); Z[4] := ComplexSub(A[0], B[0]);
Z[1] := ComplexAdd(A[1], B[1]); Z[5] := ComplexSub(A[1], B[1]);
Z[2] := ComplexAdd(A[2], B[2]); Z[6] := ComplexSub(A[2], B[2]);
Z[3] := ComplexAdd(A[3], B[3]); Z[7] := ComplexSub(A[3], B[3]);
end;
procedure FFT_10(var Z: array of TComplex);
var
A, B: array[0..4] of TComplex;
begin
A[0] := Z[0]; B[0] := Z[5];
A[1] := Z[2]; B[1] := Z[7];
A[2] := Z[4]; B[2] := Z[9];
A[3] := Z[6]; B[3] := Z[1];
A[4] := Z[8]; B[4] := Z[3];
FFT_5(A);
FFT_5(B);
Z[0] := ComplexAdd(A[0], B[0]); Z[5] := ComplexSub(A[0], B[0]);
Z[6] := ComplexAdd(A[1], B[1]); Z[1] := ComplexSub(A[1], B[1]);
Z[2] := ComplexAdd(A[2], B[2]); Z[7] := ComplexSub(A[2], B[2]);
Z[8] := ComplexAdd(A[3], B[3]); Z[3] := ComplexSub(A[3], B[3]);
Z[4] := ComplexAdd(A[4], B[4]); Z[9] := ComplexSub(A[4], B[4]);
end;
{
Synthesize the FFT by taking the even factors and the odd factors multiplied by
complex sinusoid
}
procedure SynthesizeFFT(Sofar, Radix, Remain: integer; var Y: array of TComplex);
var
GroupOffset, DataOffset, Position: integer;
GroupNo, DataNo, BlockNo, SynthNo: integer;
Omega: double;
S, CosSin: TComplex;
Synth, Trig, Z: array[0..cMaxPrimeFactor - 1] of TComplex;
// Local function
procedure InitializeTrigonomials(Radix: integer);
// Initialize trigonomial coefficients
var
i: integer;
W: double;
X: TComplex;
begin
W := 2 * pi / Radix;
Trig[0] := Complex(1.0, 0.0);
X := Complex(cos(W), -sin(W));
Trig[1] := X;
for i := 2 to Radix - 1 do
Trig[i] := ComplexMul(X, Trig[i - 1]);
end;
// Local Function
procedure FFT_Prime(Radix: integer);
// This is the general DFT, which can't be made any faster by factoring because
// Radix is a prime number
var
i, j, k, N, AMax: integer;
Re, Im: TComplex;
V, W: array[0..cMaxPrimeFactorDiv2 - 1] of TComplex;
begin
N := Radix;
AMax := (N + 1) div 2;
for j := 1 to AMax - 1 do begin
V[j].Re := Z[j].Re + Z[n-j].Re;
V[j].Im := Z[j].Im - Z[n-j].Im;
W[j].Re := Z[j].Re - Z[n-j].Re;
W[j].Im := Z[j].Im + Z[n-j].Im;
end;
for j := 1 to AMax - 1 do begin
Z[j] := Z[0];
Z[N-j] := Z[0];
k := j;
for i := 1 to AMax - 1 do begin
Re.Re := Trig[k].Re * V[i].Re;
Im.Im := Trig[k].Im * V[i].Im;
Re.im := Trig[k].Re * W[i].Im;
Im.Re := Trig[k].Im * W[i].Re;
Z[N-j].Re := Z[N-j].Re + Re.Re + Im.Im;
Z[N-j].Im := Z[N-j].Im + Re.Im - Im.Re;
Z[j].Re := Z[j].Re + Re.Re - Im.Im;
Z[j].Im := Z[j].Im + Re.Im + Im.Re;
k := k + j;
if k >= N then
k := k - N;
end;
end;
for j := 1 to AMax - 1 do begin
Z[0].Re := Z[0].Re + V[j].Re;
Z[0].Im := Z[0].Im + W[j].Im;
end;
end;
// main
begin
// Initialize trigonomial coefficients
InitializeTrigonomials(Radix);
Omega := 2 * pi / (Sofar * Radix);
CosSin := Complex(cos(Omega), -sin(Omega));
S := Complex(1.0, 0.0);
DataOffset := 0;
GroupOffset := 0;
Position := 0;
for DataNo := 0 to Sofar - 1 do begin
if Sofar > 1 then begin
Synth[0] := Complex(1.0, 0.0);
Synth[1] := S;
for SynthNo := 2 to Radix - 1 do
Synth[SynthNo] := ComplexMul(S, Synth[SynthNo - 1]);
S := ComplexMul(CosSin, S);
end;
for GroupNo := 0 to Remain - 1 do begin
if (Sofar > 1) AND (DataNo > 0) then begin
Z[0] := Y[Position];
BlockNo := 1;
repeat
inc(Position, Sofar);
Z[BlockNo] := ComplexMul(Synth[BlockNo], Y[Position]);
inc(BlockNo);
until BlockNo >= Radix;
end else begin
for BlockNo := 0 to Radix - 1 do begin
Z[BlockNo] := Y[Position];
inc(Position, Sofar);
end;
end;
case Radix of
2: FFT_2(Z);
3: FFT_3(Z);
4: FFT_4(Z);
5: FFT_5(Z);
8: FFT_8(Z);
10: FFT_10(Z);
else
// Any larger prime number than 5 (so 7, 11, 13, etc, up to cMaxPrimeFactor)
FFT_Prime(Radix);
end; //case
Position := GroupOffset;
for BlockNo := 0 to Radix - 1 do begin
Y[Position] := Z[blockNo];
Inc(Position, Sofar);
end;
GroupOffset := GroupOffset + Sofar * Radix;
Position := GroupOffset;
end;
inc(DataOffset);
GroupOffset := DataOffset;
Position := DataOffset;
end;
end;
procedure ForwardFFT(const Source: array of TComplex; var Dest: array of TComplex; Count: integer);
// Perform a FFT on the data in Source, put result in Dest. This routine works best
// for Count as a power of 2, but also works usually faster than DFT by factoring
// the series. Only in cases where Count is a prime number will this method be
// identical to regular DFT.
type
PComplexArray = ^TComplexArray;
TComplexArray = array[0..0] of TComplex;
var
i: integer;
FactorCount: integer;
SofarRadix: TIdx1FactorArray;
ActualRadix: TIdx1FactorArray;
RemainRadix: TIdx0FactorArray;
begin
if Count = 0 then exit;
// Decompose the series with length Count into FactorCount factors in ActualRadix
Factorize(Count, FactorCount, ActualRadix);
// Check if our biggest prime factor is not too large
if (ActualRadix[1] > cMaxPrimeFactor) then raise EMathError.Create(sErrPrimeTooLarge);
// Setup Sofar and Remain tables
RemainRadix[0] := Count;
SofarRadix[1] := 1;
RemainRadix[1] := Count div ActualRadix[1];
for i := 2 to FactorCount do begin
SofarRadix[i] := SofarRadix[i-1] * ActualRadix[i-1];
RemainRadix[i] := RemainRadix[i-1] div ActualRadix[i];
end;
// Reorder the series so that the elements are already in the right place for synthesis
ReorderSeries(Count{, FactorCount}, ActualRadix, RemainRadix, Source, Dest);
// Synthesize each of the FFT factored series
for i := 1 to FactorCount do
SynthesizeFFT(SofarRadix[i], ActualRadix[i], RemainRadix[i], Dest);
end;
procedure InverseFFT(const Source: array of TComplex; var Dest: array of TComplex; Count: integer);
// Perform the inverse FFT on the Source data, and put result in Dest. It performs
// the forward FFT and then divides elements by N
var
i: integer;
S: TFloat;
TmpSource: array of TComplex;
begin
if Count = 0 then exit;
// Since TmpSource is local, we do not have to free it again,
// it will be freed automatically when out of scope
SetLength(TmpSource, Count);
// Create a copy with inverted imaginary part.
for i := 0 to Count - 1 do
with Source[i] do
TmpSource[i] := Complex(Re, -Im);
ForwardFFT(TmpSource, Dest, Count);
// Scale by 1/Count, and inverse the imaginary part
S := 1.0 / Count;
for i := 0 to Count - 1 do begin
Dest[i].Re := S * Dest[i].Re;
Dest[i].Im := - S * Dest[i].Im;
end;
end;
end.
1 Like
Extra need
{ Unit Complexs
This unit implements complex number arithmic, including the basic
operations addition, substraction, multiplication, division,
magnitude and phase.
Copyright: Nils Haeck M.Sc. (email: n.haeck@simdesign.nl)
For more information visit http://ww.simdesign.nl
Original date of publication: 10 Mar 2003
This unit requires these other units:
- Math: Delphi mathematics unit
- Types: Additional mathematical variable types
****************************************************************
The contents of this file are subject to the Mozilla Public
License Version 1.1 (the "License"); you may not use this file
except in compliance with the License. You may obtain a copy of
the License at:
http://www.mozilla.org/MPL/
Software distributed under the License is distributed on an
"AS IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or
implied. See the License for the specific language governing
rights and limitations under the License.
}
unit Complexs;
interface
uses
Types, Math;
type
TFloat = single;
// Complex numbers, with precision specified in TFloat (Types unit)
TComplex = packed record
Re: TFloat; // Real part
Im: TFloat; // Imaginary part
end;
const
// Zero value
ComplexZero: TComplex = (Re: 0.0; Im: 0.0);
// Set a complex number
function Complex(Re: TFloat; Im: TFloat): TComplex;
// Add complex numbers (Result = C1 + C2)
function ComplexAdd(const C1, C2: TComplex): TComplex;
// Substract complex numbers (Result = C1 - C2)
function ComplexSub(const C1, C2: TComplex): TComplex;
// Multiply complex numbers (Result = C1 * C2)
function ComplexMul(const C1, C2: TComplex): TComplex;
// Scale complex numbers (Result = Scale * C)
function ComplexScl(Scale: TFloat; const C: TComplex): TComplex;
// Get the magnitude of the complex number C
function ComplexMag(const C: TComplex): TFloat;
// Get the phase of the complex number C (in radians, between -pi and pi)
function ComplexPhase(const C: TComplex): TFloat;
implementation
function Complex(Re: TFloat; Im: TFloat): TComplex;
// Set a complex number
begin
Result.Re := Re;
Result.Im := Im;
end;
function ComplexAdd(const C1, C2: TComplex): TComplex;
// Add complex numbers (Result = C1 + C2)
begin
Result.Re := C1.Re + C2.Re;
Result.Im := C1.Im + C2.Im;
end;
function ComplexSub(const C1, C2: TComplex): TComplex;
// Substract complex numbers (Result = C1 - C2)
begin
Result.Re := C1.Re - C2.Re;
Result.Im := C1.Im - C2.Im;
end;
function ComplexMul(const C1, C2: TComplex): TComplex;
// Multiply complex numbers (Result = C1 * C2)
begin
Result.Re := C1.Re * C2.Re - C1.Im * C2.Im;
Result.Im := C1.Im * C2.Re + C1.Re * C2.Im;
end;
function ComplexScl(Scale: TFloat; const C: TComplex): TComplex;
// Scale complex numbers (Result = Scale * C)
begin
Result.Re := Scale * C.Re;
Result.Im := Scale * C.Im;
end;
function ComplexMag(const C: TComplex): TFloat;
// Get the magnitude of the complex number C
begin
Result := sqrt(sqr(C.Re) + sqr(C.Im));
end;
function ComplexPhase(const C: TComplex): TFloat;
// Get the phase of the complex number C (in radians, between -pi and pi)
const
c2Pi = 2 * pi;
cPid2 = 0.5 * pi;
begin
// Both zero
if (C.Re = 0) and (C.Im = 0) then begin
Result := 0;
exit;
end;
// Non-zero case
if abs(C.Re) > abs(C.Im) then begin
Result := ArcTan(C.Im / C.Re); {-45 to 45 deg, 135 to -135 deg}
if C.Re < 0 then Result := Result + pi;
end else begin
Result := cPid2 - ArcTan(C.Re / C.Im); {45 to 135, -45 to -135}
if C.Im < 0 then Result := Result + pi;
end;
if Result > pi then Result := Result - c2pi;
end;
end.
Example use
private
{ Private declarations }
RArry,F1Arry: ARRAY of TComplex;
PROCEDURE TFFTView.Execute(VAR tmpA:ARRAY of DOUBLE; Sze:INTEGER);
VAR cyc:INTEGER;
fval:REAL;
BEGIN
Screen.Cursor:= crHourGlass;
ArrSze:= Sze DIV 2;
SetLength(RArry,Sze);
// Set Ch 1
FOR cyc:= 1 TO Sze DO
BEGIN
RArry[cyc-1].Im:= 0; // Set imaginary part Zero
RArry[cyc-1].Re:= tmpA[cyc-1]; // The amptlidue data of the waveform against time
END;
SetLength(F1Arry,Sze);
ForwardFFT(RArry,F1Arry,Sze); // Create the FFT data Real and Imaginary
Ch1Max:= 0;
FOR cyc:= 1 TO ArrSze DO
BEGIN
fval:= ComplexMag(F1Arry[cyc-1]); // Convert Real and Imaginary to Magnitude
IF (fval>Ch1Max) THEN Ch1Max:= fval;
F1Arry[cyc-1].Re:= fval;
END;
Screen.Cursor:= crDefault;
PaintBoxPaint(Self); // Which draws the F1Arry.Re values to get the frequency spectrum
END;
2 Likes
Many thanks, Charles. I am looking forward to the challenge!
1 Like