FFT in Delphi: Fast Fourier Transform Tutorial with Window Functions and Analysis Code

• Home
• Delphi
• FFT in Delphi: Fast Fourier Transform Tutorial with Window Functions and Analysis Code

FFT in Delphi: Fast Fourier Transform Tutorial with Window Functions and Analysis Code

This guide explains how to implement FFT (Fast Fourier Transform) in Delphi, covering everything from computing real and imaginary components to amplitude and frequency-domain analysis, with complete sample code.
It also includes implementations of various window functions such as the Hanning, Hamming, Blackman–Harris, and Flat Top windows.
This is a practical technical reference for anyone who wants to start FFT analysis in Delphi with minimal effort.

Reference: Sample code for a voice changer application that performs pitch shifting using FFT

FFT Unit (Fast Fourier Transform)

//Fast Fourier Transform
unit UFFT;

interface
uses System.SysUtils,System.Math;
type
  TFFTData=array of Extended;

//Fast Fourier Transform (input signal, output real part, output imaginary part)
procedure fft(InRe:TFFTData;var OutRe,OutIm:TFFTData);
procedure _fft(InRe:TFFTData;var OutRe,OutIm:TFFTData);

// ******************* Window Functions *****************************
// Hanning window
procedure WinHanning(var data: TFFTData);
// Hamming window
procedure WinHamming(var data: TFFTData);
// Gaussian window
procedure WinGauss(var data: TFFTData; m: integer = 1);
// Blackman–Harris window
procedure WinBlackmanHarris(var data: TFFTData);
// Blackman–Nuttall window
procedure WinBlackmanNuttall(var data: TFFTData);
// Flat Top window
procedure WinFlatTop(var data: TFFTData);
// Half-sine window
procedure WinHalfSin(var data: TFFTData);


implementation


procedure fft(InRe:TFFTData;var OutRe,OutIm:TFFTData);
var i  :Integer;
    InN: Integer; // number of input samples
    n  : Integer; // adjusted sample count (power of two)
begin
  InN:=Length(InRe);
  //If the number of samples is not a power of two, pad with zeros
  i:=1;
  while InN > Power(2,i) do inc(i);
  n:=trunc(IntPower(2,i));
  if InN < n then
  begin
    setlength(InRe,n);
    for i := InN to n-1 do InRe[i]:=0;
  end;
  //Perform FFT
  _fft(InRe,OutRe,OutIm);
end;

procedure _fft(InRe:TFFTData;var OutRe,OutIm:TFFTData);
var n:Integer;
    i:integer;
    InIm:TFFTData; //imaginary part of complex numbers
    ct1,ct2,ct3:integer;
    TmpRe,TmpIm:extended;
    nfft:array[0..3] of integer;
    fcos,fsin:TFFTData;
    tmp:extended;
    noblk:integer;
    cntb:array[0..1] of integer;
begin
  n:=Length(InRe);
  setlength(InIm,n);
  for i := 0 to n-1 do InIm[i]:=0;

  ct2:=1;
  for ct1 := 1 to length(InIm)-2 do
  begin
    TmpRe:=0;
    TmpIm:=0;
    if ct1<ct2 then
    begin
      TmpRe:=InRe[ct1-1];
      InRe[ct1-1]:=InRe[ct2-1];
      InRe[ct2-1]:=TmpRe;
      TmpIm:=InIm[ct1-1];
      InIm[ct1-1]:=InIm[ct2-1];
      InIm[ct2-1]:=TmpIm;
    end;
    ct3:=length(InIm) div 2;
    while ct3<ct2 do
    begin
      ct2:=ct2-ct3;
      ct3:=ct3 div 2;
    end;
    ct2:=ct2+ct3;
  end;
                                            //Error adjustment
  nfft[0]:=floor(Log2(length(InIm))/Log2(2)+0.0000001);
  SetLength(fcos,n);
  SetLength(fsin,n);
  fcos[0]:=1;
  fsin[0]:=0;

  for ct1 := 1 to nfft[0] do
  begin
    nfft[2]:=floor(System.math.Power(2,ct1));
    nfft[1]:=n div nfft[2];
    nfft[3]:=nfft[2] div 2;
    for ct2 := 1 to nfft[3] do
    begin
      tmp:=-Pi/nfft[3]*ct2;
      fcos[ct2]:=cos(tmp);
      fsin[ct2]:=sin(tmp);
    end;
    for ct2 := 1 to nfft[1] do
    begin
      noblk:=nfft[2]*(ct2-1);
      for ct3 := 0 to nfft[3]-1 do
      begin
        cntb[0]:=noblk+ct3;
        cntb[1]:=cntb[0]+nfft[3];
        TmpRe:=InRe[cntb[1]]*fcos[ct3]-InIm[cntb[1]]*fsin[ct3];
        TmpIm:=InIm[cntb[1]]*fcos[ct3]+InRe[cntb[1]]*fsin[ct3];
        InRe[cntb[1]]:=InRe[cntb[0]]-TmpRe;
        InIm[cntb[1]]:=InIm[cntb[0]]-TmpIm;
        InRe[cntb[0]]:=InRe[cntb[0]]+TmpRe;
        InIm[cntb[0]]:=InIm[cntb[0]]+TmpIm;
      end;
    end;
  end;
  SetLength(OutRe,Length(InRe));
  SetLength(OutIm,Length(InIm));
  Move(InRe[0],OutRe[0], Length(InRe)*SizeOf(InRe[0]));
  Move(InIm[0],OutIm[0], Length(InIm)*SizeOf(InIm[0]));
end;

//Hanning window
procedure WinHanning(var data:TFFTData);
var i,n:integer;
begin
  n:=length(data);
  for i := 0 to n-1 do
  begin
    data[i]:=(
               0.5 - 0.5*Cos(2*Pi*i/(n-1))
             )*data[i];
  end;
end;
//Hamming window
procedure WinHamming(var data:TFFTData);
var i,n:integer;
begin
  n:=length(data);
  for i := 0 to n-1 do
  begin
    data[i]:=(
               0.54 - 0.46 * Cos(2*Pi*i/(n-1))
             )*data[i];
  end;
end;
//Gaussian window
procedure WinGauss(var data:TFFTData;m:integer=1);
var i,n:integer;
begin
  n:=length(data);
  for i := 0 to n-1 do
  begin
    data[i]:=Exp(
               -2 * power(m,2) / power(n-1,2) * power(i-(n-1)/2,2)
             )*data[i];
  end;
end;
//Blackman–Harris window
procedure WinBlackmanHarris(var data:TFFTData);
var i,n:integer;
begin
  n:=length(data);
  for i := 0 to n-1 do
  begin
    data[i]:=(0.35875-0.48829*cos(2*Pi*i/(n-1))
                     +0.14128*cos(4*Pi*i/(n-1))
                     -0.01168*cos(6*Pi*i/(n-1))
             )*data[i];
  end;
end;
//Blackman–Nuttall window
procedure WinBlackmanNuttall(var data:TFFTData);
var i,n:integer;
begin
  n:=length(data);
  for i := 0 to n-1 do
  begin
    data[i]:=(0.3635819-0.4891775*cos(2*Pi*i/(n-1))
                       +0.1365995*cos(4*Pi*i/(n-1))
                       -0.0106411*cos(6*Pi*i/(n-1))
             )*data[i];
  end;
end;
//Flat Top window
procedure WinFlatTop(var data:TFFTData);
var i,n:integer;
begin
  n:=length(data);
  for i := 0 to n-1 do
  begin
    data[i]:=(1-1.930*Cos(2*Pi*i/(n-1))
               +1.290*Cos(4*Pi*i/(n-1))
               -0.388*Cos(6*Pi*i/(n-1))
               +0.032*Cos(8*Pi*i/(n-1))
             )*data[i];
  end;
end;
//Half-sine window
procedure WinHalfSin(var data:TFFTData);
var i,n:integer;
begin
  n:=length(data);
  for i := 0 to n-1 do
  begin
    data[i]:=Sin(pi*i/(n-1))*data[i];
  end;
end;

end.

Reference URL: https://gijyutsu-keisan.com/excel/numcal/fourier/fourier_2_2.php

Since FFT requires the number of input samples to be a power of two, the fft function adjusts the input length by padding zeros before calling _fft.
In addition, FFT assumes that the beginning and end of the input signal are continuous. When this condition is not met, window functions must be applied.

Using the Fast Fourier Transform

unit Unit1;

interface

uses
  Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
  Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, Vcl.ExtCtrls;

type
  TForm1 = class(TForm)
    Button1: TButton;
    Button2: TButton;
    procedure Button2Click(Sender: TObject);
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

uses UFFT,system.math;

procedure TForm1.Button1Click(Sender: TObject);
var InRe:TFFTData;
    // Real and imaginary parts of complex numbers (i = √-1)
    //     z = Re + Im * i
    OutRe,OutIm:TFFTData;
    n,i:integer;
    bs:extended;// Sampling interval (seconds)
    hz:extended;// Sampling frequency (Hz) = 1 / bs
    bhz:extended;

    OutHz:extended;  //Frequency
    OutAmp:extended; //Amplitude
begin
  n:=200;  //Number of samples
  bs:=0.01;//Sampling interval (10 ms per sample)
  hz:=1/bs;//Sampling frequency (100 Hz)

  //Create test data
  setlength(InRe,n);
  for i := 0 to n-1 do
  begin
    InRe[i]:=4*sin(i*bs*2*pi*10)+ //10Hz
             5*sin(i*bs*2*pi*20)+ //20Hz
             6*sin(i*bs*2*pi*40); //40Hz
  end;
  // If the data is not continuous at the beginning and end,
  // apply one of the window functions
  //WinHanning(InRe);        // Hanning window
  //WinHamming(InRe);        // Hamming window
  //WinGauss(InRe, 1);       // Gaussian window (1–3σ)
  //WinBlackmanHarris(InRe); // Blackman–Harris window
  //WinBlackmanNuttall(InRe);// Blackman–Nuttall window
  //WinFlatTop(InRe);        // Flat Top window
  //WinHalfSin(InRe);        // Half-sine window


  //Perform FFT
  fft(InRe,OutRe,OutIm);

  bhz:=hz/length(OutRe);
  for i := Low(OutRe) to Length(OutRe) div 2 - 1 do
  begin
    //Frequency
    OutHz :=bhz*i;
    //Amplitude = √(Re² + Im²)
    OutAmp:=Sqrt(power(OutRe[i],2)+power(OutIm[i],2));
  end;
end;

end.

For a sampling frequency Fs:
Nyquist frequency = Fs / 2
The FFT output is symmetric around the Nyquist frequency, so the maximum frequency that can be obtained from the FFT is the Nyquist frequency.
The amplitude is calculated as the square root of the sum of the squared real and imaginary components.
The phase can be obtained by taking the arctangent of (imaginary / real).

In the example above, the sampling interval is 10 ms (sampling frequency 100 Hz, Nyquist frequency 50 Hz), and the input signal is a combination of: 10 Hz (amplitude 4), 20 Hz (amplitude 5), and 40 Hz (amplitude 6).

Since FFT assumes that the first and last samples of the input data are continuous, this example applies a Hanning window to the data before calling the FFT function.

If the data is uniformly sampled, FFT can be used for frequency analysis. However, if the sampling is uneven or contains missing values, the Lomb–Scargle Periodogram(Japanese) is recommended instead.