Circular spectrum visualizer

[Michael II]

Hoi EWeiss

die von dir genutzte uSpectrum.FFT Funktion rechnet falsch.

Nimm doch eine hier aus dem Forum, zum Beispiel diese hier


Ich habe die Unit aus dem Forum etwas gekürzt (Code hier unten). Diese Unit fügst du zu deinem Projekt hinzu:

Delphi-Quellcode:

unit uDFT;

interface

uses Math;

//
// Autor Matze - siehe: https://www.delphipraxis.net/597828-post1.html
//


type
  TComplex = record
    re, im: Extended;
  end;

  TComplexArray = array of TComplex;

procedure DFT(var a: TComplexArray);

implementation

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


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


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

function MakeC(re, im: extended): TComplex;
begin
  Result.re := re;
  Result.im := im;
end;


procedure shuffle(var a: TComplexArray; n, lo: Integer);
var I, m: Integer;
    b: TComplexArray;
begin
  m := n shr 1;
  setlength(b, m);
  for I := 0 to m - 1 do
    b[i] := a[lo + i];
  for I := 0 to m - 1 do
    a[lo + i + i + 1] := a[lo + i + m];
  for I := 0 to m - 1 do
    a[lo + i + i] := b[i];
end;


procedure DoFFT(var a: TComplexArray; n, lo: Integer; w: TComplex);
var I, m: Integer;
    z, v, h: TComplex;
begin
  if n and (n - 1) = 0 then
  begin
    if n > 1 then
    begin
        m := n shr 1;
        z := MakeC(1, 0);
        for I := lo to lo + m - 1 do
        begin
            h := SubC(a[i], a[i + m]);
            a[i] := AddC(a[i], a[i + m]);
            a[i + m] := MulC(h, z);
            z:=MulC(z,w);
        end;
        v := MulC(w, w);
        DoFFT(a, m, lo, v);
        DoFFT(a, m, lo + m, v);
        shuffle(a, n, lo);
    end;
  end;
end;

procedure DFT(var a: TComplexArray);
begin
    DoFFT(a, length(a), 0, MakeC(cos(2 * Pi / length(a)),
        sin(2 * Pi / length(a))));
end;
end.

Wenn du die verwendeten Dateitypen (dein TComplex und Matzes TComplex) nicht anpassen magst, dann ersetze in uSpectrum.pas die FFT Funktion durch diese hier (Code unten). (uDFT unter uses hinzuzufügen.)

Delphi-Quellcode:

uses ….uDFT; 

…
…


procedure TSpectrum.FFT( var Dat : array of TComplex );
var a : uDFT.TComplexArray;
    i, n : integer;
begin
  n := length( Dat );
  setlength( a, n );

  for i := 0 to n-1 do
  begin
    a[i].re := Dat[i].r;
    a[i].im := Dat[i].i;
  end;

  DFT( a );

  for i := 0 to n-1 do
  begin
    Dat[i].r := a[i].re/n;
    Dat[i].i := a[i].im/n;
  end;
end;

Ich verwende für meine Programme eine iterative Version von Cooley und Tukey. Ich speichere dabei sämtliche Einheitswurzeln einmal in einer Tabelle ab und greife dann auf diese zu. Eine solche iterative Lösung ist bei der von dir gewählten Problemgrösse FFFTSize=2048 aber nur ca. 6 Mal schneller.

Selbst auf meinem langsamen Notebook benötigt obige Lösung nur ca. 2/1000 Sekunden.


Viel Spass beim Codieren.