Digital Signal Processing 3D Graphics & Rendering Programming Download Natalia
Komentarze (1)  Wizyty (3 364)  Aktualnie osób (5)  Tutaj (1) 
User: Hasło:

 Fast Fourier Transform

Reklamy Google

Opis klasy FFT Delphi - TFastFourier, wyniki szybkości obliczeń dla różnych wielkości danych oraz przykład zastosowania. Do pobrania gotowy komponent obliczający transformatę.

Krótka charakterystyka klasy TFastFourier

Komponenet FastFourier zawiera w sobie Szybką Transformatę Fouriera. Pamiec jest przydzielana tutaj dynamicznie. Dwie zespolone tablice: wejsciowa i wyjsciowa sa typu zespolonego (InBuffer, OutBuffer). Okna czasowe jakie mozna stosowac okreslaja stałe w definicji modulu z przedrostkiem 'ok'. Najwazniejsze procedury to FFT oraz IFFT - obie bezparametrowe. Dla wypelnionej tablicy wejsciowej InBuffer procedura FFT wykonuje szybką transformate Fouriera a wynik zapisuje do tablicy OutBuffer. IFFT robi transformate odwrotna dla danych wejsciowych. Po dokonaniu transformaty mamy do dyspozycji funkcje: Normalize, Real, Imag, Absol, Power, Energy, Phase, LogEnergy - chyba nie trzeba tlumaczyc ich znaczenia. Ich argumentami sa indeksy tablicy OutBuffer.

Implementacja klasy TFastFourier

Przedstawiony komponent FastFourier dla Delphi zawiera stale dla okien czasowych, typ TComplex (liczba zespolona) oraz klasa calego mechanizmu.
Definicje te wygladaja nastepujaco:

 

type
TComplex = Record
    Re : Extended;
    Im : Extended;
end;

type
TFastFourier = class(TComponent)
public
    InBuffer : array of TComplex;
    OutBuffer : array of TComplex;
    procedure FFT;
    procedure IFFT;
    function Real(idx : LongWord) : Extended;
    function Imag(idx : LongWord) : Extended;
    function Absol(idx : LongWord) : Extended;
    function Power(idx, SamplFreq : LongWord) : Extended;
    function Energy(idx : LongWord) : Extended;
    function Phase(idx : LongWord) : Extended;
published
    property NumSamples : LongWord read FNumSamples write SetNumSamples;
end;

Testy czasowe komponentu TFastFourier na AMD 333 MHz, 128 MB RAM:

  • dla N=8192 t=110 milisek
  • dla N=16384 t=220 milisek
  • dla N=32768 t~380 milisek
  • dla N=65536 t=880 milisek
  • dla N=131072 t=1.8 sek
  • dla N=262144 t=4.1 sek
  • dla N=524288 t=8,6 sek
  • dla N=1048576 t=17,6 sek
  • dla N=2097152 t=41,9 sek


    • Program prezentacyjny komponentu TFastFourier. W pasku statusu pokazywane są czasy obliczeń.

    Pliki do pobrania:

    myFFT.exe (379 kB) - aplikacja testująca komponent TFastFourier (459 pobrań)

    Przykład użycia TFastFourier

    uses TFastFourier;
    var FFT : TFastFourier;

    FFT.NumSamples:=128;
    For i:=0 to NumSamples-1 do
    FFT.InBuffer[i].re:=sin(0.5*i);
    FFT.FFT;

    For i:=0 to NumSamples-1 do
    Wykres.Data[i]:=Energy(i);
    Wykres.Rysuj;

    Algorytm transformaty Fouriera w wersji Fast - algorytm o podstawie 2

    Tutaj postaram się przedstawię algorytm szybkiej transformaty Fouriera. Warunkiem do algorytmu FFT jest aby liczba próbek do przetworzenia byla całkowitą potegą liczby 2 - czyli liczba próbek poddawana transformacie moęe być równa 2, 4, 8 ... 128, 256, 512, 1024 ... itd.
    Nie bedę tutaj przytaczała calego wyprowadzenia. Powiem tylko, ęe algorytm wymaga nastepujących procedur i funkcji, których definicje zamieszczam w dalszej części:

  • IsPowerOfTwo(x : LongWord) : Boolean

    - Funkcja zwraca wartosc logiczna True, jesli jej argument x jest potega liczby 2.
  • NumberOfBitsNeeded(Ile : LongWord) : LongWord

    - Funkcja zwraca liczbe bitow potrzebnych do zakodowania liczby argumentu Ile. Ile musi byc potega liczby 2.
  • ReverseBits(index, NumBits : LongWord) : LongWord

    - Funkcja robi "odwrocenie bitowe" dla 'index' przy znajomosci liczby potrzebnych bitow NumBits.
  • FourierTransform ( Angle: Extended)

    - Procedura szybkiej transformaty. To jest wlasnie najwazniejsza funkcja. W zaleznosci od kata Angle procedura robi FFT lub IFFT.
  • FFT

    - Procedura FFT jest wywolaniem procedury FourierTransform z argumentem 2*PI
  • IFFT

    - Procedura IFFT jest wywolaniem procedury FourierTransform z argumentem -2*PI. Dodatkowo przetransformowane probki nalezy podzielic przez ich ilosc (wynika to z definicji Transformaty Fouriera).


    Unit FastFourier;

          {*********************************************************}
          {***   FAST FOURIER TRANSFORM   Component              ***}
          {***   Autor : Natalia Walkowicz, grudzien 2003        ***}
          {***   Maxymalna liczba probek: N = 2 147 483 648      ***}
          {*********************************************************}

    interface

    uses
      Windows, Messages, SysUtils, Classes, Math, Complex;

    type
         TFastFourier = class(TComponent)
         private
              FNumSamples    : LongWord;
              function  IsPowerOfTwo(x : LongWord) : Boolean;
              function  PowerOf2(Liczba : LongWord) : LongWord;
              function  NumberOfBitsNeeded(Ile : LongWord) : LongWord;
              function  ReverseBits(index, NumBits : LongWord) : LongWord;
              procedure FourierTransform(Angle : Extended);
              procedure SetNumSamples(Value : LongWord);
         public
              InBuffer  : array of TComplex;
              OutBuffer : array of TComplex;
              procedure FFT;
              procedure IFFT;
              function  NumberOfSamples(Count : LongWord) : LongWord;
              function  Real(idx : LongWord) : Extended;
              function  Imag(idx : LongWord) : Extended;
              function  Absol(idx : LongWord) : Extended;
              function  Energy(idx, SamplFreq : LongWord) : Extended;
              function  Power(idx : LongWord) : Extended;
              function  Phase(idx : LongWord) : Extended;
              function  Frequency(idx, SamplFreq : LongWord) : Extended;
              procedure Free;
              constructor Create(AOwner : TComponent); override;
              destructor  Destroy; override;
         published
              property NumSamples : LongWord read FNumSamples write SetNumSamples;
    end;

    procedure Register;

    implementation

    procedure Register;
    begin
         RegisterComponents('Engineering', [TFastFourier]);
    end;

    constructor TFastFourier.Create(AOwner : TComponent);
    begin
         inherited Create(AOwner);
         NumSamples:=0;
    end;

    destructor TFastFourier.Destroy;
    begin
         SetLength(InBuffer,0);
         SetLength(OutBuffer,0);
         inherited;
    end;

    procedure TFastFourier.SetNumSamples(Value : LongWord);
    begin
         FNumSamples:=Value;
         SetLength(InBuffer,Value);
         SetLength(OutBuffer,Value);
    end;

    function TFastFourier.NumberOfSamples(Count: LongWord): LongWord;
    var   i, wynik, potega :  LongWord;
    begin
         if IsPowerOfTwo(Count) then
            wynik:=Count;

         if not IsPowerOfTwo(Count) then
         begin
              potega := 1;
              i:=0;
              while ((i<=31) and (potega           begin
                   potega := potega SHL 1;
                   i:=i+1;
              end;
         end;

         Result:=potega;
    end;

    function TFastFourier.PowerOf2(Liczba : LongWord) : LongWord;
    var i, potega : LongWord;
    begin
         potega:=1;

         for i :=1 to Liczba do
             potega:=potega*2;

         Result:=potega;
    end;

    function TFastFourier.IsPowerOfTwo(x : LongWord) : Boolean;
    var   i, y :  LongWord;
    begin
         y:=2;
         for i:=1 to 31 do
         begin
              if x = y then
              begin
                   Result := True;
                   exit;
              end;
              y:=y SHL 1;
         end;

         Result := False;
    end;

    function TFastFourier.NumberOfBitsNeeded(Ile : LongWord) : LongWord;
    var i : LongWord;
    begin
         for i:=0 to 31 do
         begin
              if (Ile AND (1 SHL i)) <> 0 then
              begin
                   Result:=i;
                   Exit;
              end;
         end;
    end;

    function TFastFourier.ReverseBits(index, NumBits : LongWord) : LongWord;
    var i, rev : LongWord;
    begin
         rev:=0;
         for i:=0 to NumBits-1 do
         begin
              rev:=(rev SHL 1) OR (index AND 1);
              index:=index SHR 1;
         end;

         Result:=rev;
    end;

    procedure TFastFourier.FourierTransform(Angle : Extended);
    var NumBits, i, j, k, n, BlockSize, BlockEnd : LongWord;
        delta_angle, delta_ar                    : Extended;
        alpha, beta                              : Extended;
        tr, ti, ar, ai                           : Extended;
    begin
         NumBits:=NumberOfBitsNeeded(FNumSamples);

         for i:=0 to FNumSamples-1 do
         begin
              j:=ReverseBits(i, NumBits);
              OutBuffer[j]:=InBuffer[i];
         end;

         BlockEnd:=1;
         BlockSize:=2;
         while BlockSize<=FNumSamples do
         begin
              delta_angle:=Angle/BlockSize;
              alpha:=sin(0.5* delta_angle);
              alpha:=2.0*alpha*alpha;
              beta:=sin(delta_angle);

              i:=0;
              while i           begin
                   ar:=1.0;    // cos(0)
                   ai:=0.0;    // sin(0)

                   j:=i;
                   for n:=0 to BlockEnd-1 do
                   begin
                        k:=j+BlockEnd;
                        tr:=ar*OutBuffer[k].Re-ai*OutBuffer[k].Im;
                        ti:=ar*OutBuffer[k].Im+ai*OutBuffer[k].Re;
                        OutBuffer[k].Re:=OutBuffer[j].Re-tr;
                        OutBuffer[k].Im:=OutBuffer[j].Im-ti;
                        OutBuffer[j].Re:=OutBuffer[j].Re+tr;
                        OutBuffer[j].Im:=OutBuffer[j].Im+ti;
                        delta_ar:=alpha*ar+beta*ai;
                        ai:=ai-(alpha*ai-beta*ar);
                        ar:=ar-delta_ar;
                        INC(j);
                   end;
                   i:=i+BlockSize;
              end;
              BlockEnd:=BlockSize;
              BlockSize:=BlockSize SHL 1;
         end;
    end;

    procedure TFastFourier.FFT;
    begin
         FourierTransform(2*PI);
    end;

    procedure TFastFourier.IFFT;
    var i: LongWord;
    begin
         FourierTransform(-2*PI);
         for i:=0 to FNumSamples-1 do
         begin
              OutBuffer[i].Re:=OutBuffer[i].Re/FNumSamples;
              OutBuffer[i].Im:=OutBuffer[i].Im/FNumSamples;
         end;
    end;

    function TFastFourier.Imag(idx: LongWord): Extended;
    var z : TComplex;
    begin
         z:=OutBuffer[idx];
         Result:=z.Im;
    end;

    function TFastFourier.Real(idx: LongWord): Extended;
    var z : TComplex;
    begin
         z:=OutBuffer[idx];
         Result:=z.Re;
    end;

    function TFastFourier.Absol(idx : LongWord) : Extended;
    var z : TComplex;
    begin
         z:=OutBuffer[idx];
         if idx=0 then
            Result:=sqrt(z.Re*z.Re+z.Im*z.Im)/FNumSamples;
         if idx>0 then
            Result:=sqrt(z.Re*z.Re+z.Im*z.Im)/FNumSamples*2;
    end;

    function TFastFourier.Power(idx: LongWord): Extended;
    var z : TComplex;
    begin
         z:=OutBuffer[idx];
         if idx=0 then
            Result:=(z.Re*z.Re+z.Im*z.Im)/FNumSamples;
         if idx>0 then
            Result:=(z.Re*z.Re+z.Im*z.Im)/FNumSamples*2;
    end;

    function TFastFourier.Energy(idx, SamplFreq : LongWord): Extended;
    var x : Extended;
    begin
         x:=Power(idx);
         Result:=x*x/FNumSamples*SamplFreq;
    end;

    function TFastFourier.Phase(idx: LongWord): Extended;
    var re, im : Extended;
    begin
         re:=Real(idx);
         im:=Imag(idx);
         if abs(re)>1e-5 then
            Result:=ArcTan(im/re)
         else
            Result:=0;
    end;

    function TFastFourier.Frequency(idx, SamplFreq: LongWord): Extended;
    begin
         if (idx>=0) or (idx         Result:=SamplFreq/FNumSamples*idx
         else
            Result:=0;
    end;

    procedure TFastFourier.Free;
    begin
         SetLength(InBuffer,0);
         SetLength(OutBuffer,0);
    end;

    end.

    Pliki do pobrania:

    FastFourier.rar (3.04 kB) - komponent FFT dla Delphi (519 pobrań)

    Visitors

    free counters

    Reklamy Google

    Komentarze (1)

    Zapraszam do pisania komentarzy w tym temacie.

    #327 - do malinki
       witam mam takie pytanie i takie do myslemioa malinka jest dobra choc obrazki nie wchodzw szystki po zainstalowaniu malinka sie zacina a co do drugiego pytania  powinniscie wymyslec cos takiego na malinmce ze na przyklad aby szlo pisac pozdrowienia dla danego radia na malince bez  potzrby wejscia na strone radiuaikonsoli pozdrowien

    Dodano przez: mmarek Data: 2009-01-03 16:46:51

    Napisz komentarz

    Pamiętaj, aby wypełnić wszystkie pola.

    Temat: Podpis:

    • .....:: E-mail to Natalia ::.....