Problem bei FFT

[3_of_8]

Das ist meine momentan (funktionierende) Implementation einer DFT. Wenn möglich wird eine FFT angewandt. Ich empfehle, das ganze nicht mit mehr als 2000 Werten aufzurufen, es sei denn, die Anzahl der Werte ist eine Zweierpotenz oder lässt sich sehr oft durch 2 teilen.

Aufgerufen wird das einfach mit DFT(a), wobei a ein Array of Integer, Single oder ein TComplexArray ist.

Nach dem Aufruf enthält das Array die Amplitudenwerte der einzelnen Frequenzen, in Index 1 steht der Wert der Frequenz (Samplingrate/Abtastwerte) Hz, in Index 2 2*(Samplingrate/Abtastwerte) Hz und in Index I steht der Wert für I*(Samplingrate/Abtastwerte) Hz. Man sollte nur Werte bis zum Index (Abtastwerte/2), also bis (Samplingrate/2) Hz, der sogenannten Nyquist-Frequenz abfragen, da einem ansonsten der Alias-Effekt auftritt.

(Vielleicht wäre das was für die CL?)

Delphi-Quellcode:

unit complex;

interface

uses Math;

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

  TComplexArray=array of TComplex;

function AddC(a, b: TComplex): TComplex; overload;
function AddC(a: TComplex; b: Single): TComplex; overload;
function SubC(a, b: TComplex): TComplex; overload;
function SubC(a: TComplex; b: Single): TComplex; overload;
function MulC(a, b: TComplex): TComplex; overload;
function MulC(a: TComplex; b: Single): TComplex; overload;
function MakeC(re, im: Single): TComplex;

procedure DFT(var a: TComplexArray); overload;
procedure DFT(var a: array of Single); overload;
procedure DFT(var a: array of Integer); overload;

implementation

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

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

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

function SubC(a: TComplex; b: Single): TComplex;
begin
  Result.re:=a.re-b;
  Result.im:=a.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 MulC(a: TComplex; b: Single): TComplex;
begin
  Result.re:=a.re*b;
  Result.im:=a.im*b;
end;

function MakeC(re, im: Single): 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 DoDFT(var a: TComplexArray; n, lo: Integer; w: TComplex);
var b, f: array of TComplex;
    I, J, index: Integer;
begin
  setlength(b, n);
  for I:=0 to n-1 do
  begin
    b[I]:=a[I+lo];
    a[I+lo]:=MakeC(0, 0);
  end;
  setlength(f, n);
  f[0]:=MakeC(1, 0);
  f[1]:=w;
  for I:=2 to n-1 do
    f[I]:=MulC(f[I-1], w);
  for I:=0 to n-1 do
  begin
    index:=0;
    for J:=0 to n-1 do
    begin
      a[I+lo]:=AddC(a[I+lo], MulC(b[J], f[index]));
      inc(index, I);
      if index>n then dec(index, n);
    end;
  end;
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 else if n>1 then
  begin
    DoDFT(a, n, lo, w);
  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;

procedure DFT(var a: array of Single);
var I, index: Integer;
    b: TComplexArray;
begin
  setlength(b, length(a));
  for I:=0 to high(a) do b[I]:=MakeC(a[I], 0);
    DoFFT(b, length(b), 0, MakeC(cos(2*Pi/length(b)),
        sin(2*Pi/length(b))));
  for I:=0 to high(b) do a[I]:=sqrt(sqr(b[I].re)+sqr(b[I].im));
end;

procedure DFT(var a: array of Integer);
var I: Integer;
    b: TComplexArray;
begin
  setlength(b, length(a));
  for I:=0 to high(a) do b[I]:=MakeC(a[I], 0);
    DoFFT(b, length(b), 0, MakeC(cos(2*Pi/length(b)),
        sin(2*Pi/length(b))));
  for I:=0 to high(b) do a[I]:=trunc(sqrt(sqr(b[I].re)+sqr(b[I].im)));
end;

end.

EDIT: Hier noch eine Funktion zum Analysieren der vorherrschenden Frequenz in einem Tonabschnitt.

Delphi-Quellcode:

function GetFrequency(Data: array of Single; SampleRate, SampleIndex, 
    Samples: Integer): Integer;
var I, n, max: Integer;
    a: array of Single;
begin
  max:=0;
  if length(Data)-SampleIndex<Samples then
    raise EAccessViolation.Create('Sound index or length out of bounds.');
  setlength(a, Samples);
  for I:=0 to Samples-1 do
    a[I]:=Data[SampleIndex+I];
  dft(a);
  for I:=0 to (n shr 1)-1 do if a[I]>a[max] then
  begin
    max:=I;
  end;
  Result:=trunc(SampleRate/Samples*max);
end;