Jump to content
CHackbart

Simple FFT

By CHackbart, in Algorithms, Data Structures and Class Design

Recommended Posts

CHackbart    13

CHackbart
Posted

Hi,

 

I currently try to write a fft visualizer (without any third party), but the last time I heard from dft/fft is about 20 years ago. This is what I did from scratch, but the implementation does not seem to be correct. Maybe you have an idea?

  

unit umain;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, Forms, Controls, Graphics, Dialogs, ExtCtrls,
  bgrabitmap, BGRABitmapTypes, uacinerella;

type
  TSpectrum = (tspWaveform, tspFFT, tspLogarithmic, tsp3d);

  TComplex = record
    re, im: double;
  end;
  TComplexArray = array of TComplex;
  TDoubleArray = array of double;

  { TForm1 }

  TForm1 = class(TForm)
    Image1: TImage;
    procedure FormActivate(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure FormDestroy(Sender: TObject);
  private
    FInstance: PAc_instance;
    FInput: TFileStream;
    FAudiodecoder: PAc_decoder;
    FBitmap: TBGRABitmap;
    FSpectrumPalette: array [0 .. 255] of TBGRAPixel;
    FSpectrum: TSpectrum;
    FSpectrumPos: integer;
    FHanningWindow: TComplexArray;
    function ReadProc(Buf: PByte; Size: integer): integer;
    procedure ProcessAudio(Buf: PByte; Size: integer);
    procedure ProcessData();

    procedure CreatePalette;
  public
  end;

var
  Form1: TForm1;

implementation

uses Math;

{$R *.lfm}

function read_proc(Sender: Pointer; buf: PByte; size: integer): integer; cdecl;
begin
  Result := TForm1(Sender).ReadProc(buf, size);
end;

function cInit(const re, im: double): TComplex;
begin
  Result.re := re;
  Result.im := im;
end;

function cadd(const a, b: TComplex): TComplex;
begin
  Result.re := a.re + b.re;
  Result.im := a.im + b.im;
end;

function csub(const a, b: TComplex): TComplex;
begin
  Result.re := a.re - b.re;
  Result.im := a.im - b.im;
end;

function cmul(const a, b: TComplex): TComplex;
begin
  Result.re := a.re * b.re - a.im * b.im;
  Result.im := b.re * a.im + a.re * b.im;
end;

function cabs(const a: TComplex): double;
begin
  Result := sqrt(a.re * a.re + a.im * a.im);
end;

function getBit(const Val: DWord; const BitVal: byte): boolean;
begin
  Result := (Val and (1 shl BitVal)) <> 0;
end;

function enableBit(const Val: DWord; const BitVal: byte; const SetOn: boolean): DWord;
begin
  Result := (Val or Int64(1 shl BitVal)) xor (integer(not SetOn) shl BitVal);
end;

function mirrorTransform(n, m: integer): integer;
var
  i: integer;
begin
  Result := 0;
  for i := 0 to m - 1 do
    Result := enableBit(Result, m - 1 - i, getBit(n, i));
end;

function doPermutation(Source: TComplexArray; m: integer): TComplexArray;
var
  i, n: integer;
begin
  n := length(Source);
  SetLength(Result, n);

  for i := 0 to n - 1 do
  begin
    Result[i] := Source[mirrorTransform(i, m)];
  end;
end;

function doStep(k, M: longint; prev: TComplexArray): TComplexArray;
var
  expTerm, substractTerm: TComplex;
  q, j, offset: longint;
begin
  offset := system.round(intpower(2, M - k));

  SetLength(Result, length(prev));

  for q := 0 to system.round(intpower(2, k - 1) - 1) do
  begin // For each block of the matrix
    for j := 0 to (offset - 1) do
    begin // Fo each line of this block
      // First half
      Result[q * 2 * offset + j] :=
        cadd(prev[q * 2 * offset + j], prev[q * 2 * offset + j + offset]);
      // Second half
      expTerm.re := cos((j * PI) / offset);
      expTerm.im := sin((j * PI) / offset);
      substractTerm := csub(prev[q * 2 * offset + j], prev[q * 2 * offset + j + offset]);
      Result[q * 2 * offset + j + offset] := cmul(expTerm, substractTerm);
    end;
  end;
end;

function doFFT(g: TComplexArray; order: integer): TComplexArray;
var
  previousRank, nextRank: TComplexArray;
  i: integer;
begin

  previousRank := g;
  for i := 1 to order do
  begin
    nextRank := doStep(i, order, previousRank);
    previousRank := nextRank;
  end;
  Result := doPermutation(nextRank, order);
end;

function WindowHann(const ACount: integer): TComplexArray;
var
  i: integer;
begin
  SetLength(Result, ACount);
  for i := 0 to ACount - 1 do
    Result[i] := cInit(0.5 * (1.0 - cos((2.0 * PI * i) / ACount)), 0);
end;

function CalculateSpectrum(const ASignalWindow, AData: TComplexArray): TDoubleArray;
var
  dftResult, tempData: TComplexArray;
  len: integer;
  i: integer;
begin
  len := length(ASignalWindow);
  SetLength(tempData, len);
  for i := 0 to Len - 1 do
  begin
    if i < len then
      tempData[i] := cmul(AData[i], ASignalWindow[i])
    else
      tempData[i] := cInit(0, 0);
  end;
  dftResult := doFFT(tempData, 10);
  Setlength(result, Length(dftResult));
  for i := 0 to length(result) - 1 do
    result[i] := cabs(dftResult[i]);
end;

{ TForm1 }

procedure TForm1.FormCreate(Sender: TObject);
var
  i: integer;
  info: TAc_stream_info;
begin
  FAudiodecoder := nil;
  FInput := TFileStream.Create('Tones_100_20000_incrementing.ac4.trp', fmOpenRead);
  FInstance := ac_init();
  ac_open(FInstance, self, nil, @read_proc, nil, nil, nil);
  for i := 0 to FInstance^.stream_count - 1 do
  begin
    ac_get_stream_info(FInstance, i, @info);
    if info.stream_type = AC_STREAM_TYPE_AUDIO then
    begin
      FAudiodecoder := ac_create_decoder(FInstance, i);
      break;
    end;
  end;
  FBitmap := TBGRABitmap.Create;
  FBitmap.SetSize(640, 480);
  CreatePalette;
  FSpectrum := TSpectrum.tspLogarithmic;
  FHanningWindow := WindowHann(1024);
end;

procedure TForm1.FormActivate(Sender: TObject);
begin
  while not application.terminated do
  begin
    ProcessData();
    application.ProcessMessages;
  end;
end;

procedure TForm1.FormDestroy(Sender: TObject);
begin
  ac_free_decoder(FAudiodecoder);
  ac_close(FInstance);
  FInput.Free;
  FBitmap.Free;
end;


procedure TForm1.CreatePalette;
var
  i: integer;
begin
  for i := 1 to 127 do
    with FSpectrumPalette[i] do
    begin
      Alpha := 255;
      Blue := 0;
      Green := 256 - 2 * i;
      Red := 2 * i;
    end;

  for i := 0 to 31 do
  begin
    with FSpectrumPalette[128 + i] do
    begin
      Alpha := 255;
      Red := 0;
      Green := 0;
      Blue := 8 * i;
    end;
    with FSpectrumPalette[128 + 32 + i] do
    begin
      Alpha := 255;
      Red := 8 * i;
      Green := 0;
      Blue := 0;
    end;
    with FSpectrumPalette[128 + 64 + i] do
    begin
      Alpha := 255;
      Red := 255;
      Green := 8 * i;
      Blue := 8 * (31 - i);
    end;
    with FSpectrumPalette[128 + 96 + i] do
    begin
      Alpha := 255;
      Red := 255;
      Green := 255;
      Blue := 8 * i;
    end;
  end;
end;

function TForm1.ReadProc(Buf: PByte; Size: integer): integer;
begin
  Result := FInput.Read(Buf^, Size);
end;

procedure TForm1.ProcessData;
var
  pack: PAc_package;
begin
  if not assigned(FAudioDecoder) then
    exit;
  pack := ac_read_package(FInstance);
  try
    if (pack = nil) or (pack^.stream_index <> FAudiodecoder^.stream_index) then
      exit;
    if (ac_decode_package(pack, FAudiodecoder) > 0) then
      ProcessAudio(FAudiodecoder^.buffer, FAudiodecoder^.buffer_size);
  finally
    ac_free_package(pack);
  end;
end;

procedure TForm1.ProcessAudio(Buf: PByte; Size: integer);

const
  BANDS = 28;

var
  X, Y, Y1, V, B0, B1, SC: integer;
  Sum: single;
  i, z, w: integer;
  Step: integer;
  ptr: PSmallInt;
  signal: TComplexArray;
  fft: array of double;

  function FixRange(const Y: integer): integer;
  begin
    Result := Y * 127 div FBitmap.Height;
  end;

begin
  if Size = 0 then
    exit;

  if FSpectrum <> tsp3d then
   FBitmap.FillTransparent;

  Y := 0;
  X := 0;


  if FSpectrum <> tspWaveForm then
  begin
    ptr := PSmallInt(Buf);
    Step := (Size div 4) div Length(FHanningWindow);
    if Step < 1 then Step := 1;

    SetLength(signal, length(FHanningWindow));
    for i := 0 to high(signal) do
    begin
      signal[i].re := (32768-ptr^)/65536;
      signal[i].im := 0;

      Inc(ptr, step);
    end;
    fft := CalculateSpectrum(FHanningWindow, signal);
  end;

  Step := (Size div 4) div FBitmap.Width;
  if Step < 1 then
    Step := 1;
  ptr := PSmallInt(Buf);

  case FSpectrum of
    tspWaveform:
    begin
      for x := 0 to FBitmap.Width - 1 do
      begin
        V := (32767 - ptr^) * FBitmap.Height div 65536;
        Inc(ptr, step);
        if X = 0 then
          Y := V;
        repeat
          if Y < V then
            Inc(Y)
          else if Y > V then
            Dec(Y);
          FBitmap.SetPixel(X, Y,
            FSpectrumPalette[FixRange(abs(Y - FBitmap.Height div 2) * 2 + 1)]);
        until Y = V;
      end;
    end;
    TSpectrum.tspFFT: // "normal" FFT
    begin
      Y1 := 0;
      for X := 0 to (FBitmap.Width div 2) - 1 do
      begin
        Y := Trunc(sqrt(fft[X + 1]) * 3 * FBitmap.Height - 4);
        if Y > FBitmap.Height then
          Y := FBitmap.Height; // cap it

        Y1 := (Y + Y1) div 2;
        if (X > 0) and (Y1 > 0) then
          // interpolate from previous to make the display smoother
          while (Y1 >= 0) do
          begin
            FBitmap.SetPixel(X * 2 - 1, FBitmap.Height - Y1 - 1,
              FSpectrumPalette[FixRange(Y1 + 1)]);
            Dec(Y1);
          end;

        Y1 := Y;
        while (Y >= 0) do
        begin
          FBitmap.SetPixel(X * 2, FBitmap.Height - Y - 1,
            FSpectrumPalette[FixRange(Y + 1)]);
          // draw level
          Dec(Y);
        end;
      end;
    end;
    TSpectrum.tspLogarithmic: // logarithmic, acumulate & average bins
    begin
      B0 := 0;
      for X := 0 to BANDS - 1 do
      begin
        Sum := 0;
        B1 := Trunc(Power(2, X * 10.0 / (BANDS - 1)));
        if B1 > 1023 then
          B1 := 1023;
        if B1 <= B0 then
          B1 := B0 + 1; // make sure it uses at least 1 FFT bin
        SC := 10 + B1 - B0;

        while B0 < B1 do
        begin
          Sum := Sum + fft[1 + B0];
          Inc(B0);
        end;

        Y := Trunc((sqrt(Sum / log10(SC)) * 1.7 * FBitmap.Height) - 4);
        // scale it
        if Y > FBitmap.Height then
          Y := FBitmap.Height; // cap it

        while (Y >= 0) do
        begin
          w := Trunc(0.9 * (FBitmap.Width / BANDS));
          for z := 0 to w - 1 do
            FBitmap.SetPixel(X * (FBitmap.Width div BANDS) + z, FBitmap.Height - Y - 1,
              FSpectrumPalette[FixRange(Y + 1)]);
          Dec(Y);
        end;
      end;
    end;
    TSpectrum.tsp3d: // "3D"
    begin
      for X := 0 to FBitmap.Height - 1 do
      begin
        Y := Trunc(sqrt(fft[X + 1]) * 3 * FBitmap.Height);
        // scale it (sqrt to make low values more visible)
        if Y > FBitmap.Height then
          Y := FBitmap.Height; // cap it
        if Y < 0 then
          Y := 0;

        if (FSpectrumPos < FBitmap.Width) and (X < FBitmap.Height) then
          FBitmap.SetPixel(FSpectrumPos, X, FSpectrumPalette[128 + FixRange(Y)]);
        // plot it
      end;
      // move marker onto next position
      FSpectrumPos := (FSpectrumPos + 1) mod FBitmap.Width;
      for X := 0 to FBitmap.Height - 1 do
        if (FSpectrumPos < FBitmap.Width) and (X < FBitmap.Height) then
          FBitmap.SetPixel(FSpectrumPos, X, FSpectrumPalette[255]);
    end;
  end;
  Image1.Picture.Assign(FBitmap);
end;

end.

My guess is that the way I fill the complex array with the signal data is not correct. 

Share this post

Link to post

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
Go To Topic Listing Algorithms, Data Structures and Class Design
×