Atmel FFT

5.2.3 FFT (Fast Fourier Transformation) (in Arbeit)

Delphi Routine

FFT mit Realzahlen
Const
  Log2n = 5;
  n = 32; // 2Log2n

Type
  RComplex = Record re,im: Real End;
  RArray = Array [0..n-1] of Real;
  RCArray = Array [0..n-1] of RComplex;

Var
  xr: RCArray;
  yr: RArray;

Procedure FFT(Var fr:RCArray; Var pw:RArray);
Var k1,k2,k3,k4,i1,i2,i3,i4,i5:SmallInt;
    t1,t2:Double; ft:RCArray;
Begin
  // Init Sinus/Cosinus
  For i1 := 0 to n-1 do Begin
    t1 := i1 * 2 * pi / n;
    ft[i1].re := Cos(t1);
    ft[i1].im := Sin(t1);
  End;
  // Bitreversing
  i2 := 0;
  For i1 := 0 to n-2 do Begin
    If i2 >= i1 Then Begin
      t1 := fr[i2].re;
      fr[i2].re := fr[i1].re;
      fr[i1].re := t1;
    End;
    k1 := n shr 1;
    While i2 >= k1 do Begin
      i2 := i2 - k1;
      k1 := k1 shr 1;
    End;
    i2 := i2 + k1;
  End;
  // FFT
  k1 := 1;
  k2 := n shr 1;
  k3 := log2n - 1;
  k4 := 1;
  While k1 <> n do Begin
    For i4 := 0 to k2-1 do Begin
      For i5 := 0 to k1-1 do Begin
        i1 := (i4 shl k4) + i5;
        i2 := i1 + k1;
        i3 := i5 shl k3;
        t1 := fr[i2].re * ft[i3].re - fr[i2].im * ft[i3].im;
        t2 := fr[i2].re * ft[i3].im + fr[i2].im * ft[i3].re;
        fr[i2].re := fr[i1].re - t1;
        fr[i2].im := fr[i1].im - t2;
        fr[i1].re := fr[i1].re + t1;
        fr[i1].im := fr[i1].im + t2;
      End;
    End;
    k1 := k1 shl 1;
    k2 := k2 shr 1;
    k3 := k3 - 1;
    k4 := k4 + 1;
  End;
  // Normierung
  For i1 := 0 to n-1 do Begin
    fr[i1].re := fr[i1].re / n;
    fr[i1].im := fr[i1].im / n;
  End;
  // Powerspektrum
  For i1 := 0 to n-1 do Begin
    t1 := fr[i1].re * fr[i1].re + fr[i1].im * fr[i1].im;
    pw[i1] := Sqrt(t1);
  End;
End;

Daraus abgeleitet die Integer FFT. Der Trick liegt in der Variablen fb (rote Schrift)

Delphi Routine

FFT mit Integerzahlen
Const
  Log2n = 5;
  n = 32; // 2Log2n
  Log2fb = 14;
  fb = 16384; // 2Log2fb

Type
  IComplex = Record re,im: Integer End;
  IArray = Array [0..n-1] of Integer;
  ICArray = Array [0..n-1] of IComplex;

Var
  xr: ICArray;
  yr: IArray;

Procedure IntFFT(Var fr:ICArray; Var pw:IArray);
Var k1,k2,k3,k4,i1,i2,i3,i4,i5:SmallInt;
    t1,t2:Integer; ft:ICArray;
Begin
  // Init Sinus/Cosinus
  For i1 := 0 to n-1 do Begin
    ft[i1].re := Round(Cos(i1*2*pi/n) * fb);
    ft[i1].im := Round(Sin(i1*2*pi/n) * fb);
  End;
  // Bitreversing
  i2 := 0;
  For i1 := 0 to n-2 do Begin
    If i2 >= i1 Then Begin
      i3 := fr[i2].re;
      fr[i2].re := fr[i1].re;
      fr[i1].re := i3;
    End;
    k1 := n shr 1;
    While i2 >= k1 do Begin
      i2 := i2 - k1;
      k1 := k1 shr 1;
    End;
    i2 := i2 + k1;
  End;
  // FFT
  k1 := 1;
  k2 := n shr 1;
  k3 := log2n - 1;
  k4 := 1;
  While k1 <> n do Begin
    For i4 := 0 to k2-1 do Begin
      For i5 := 0 to k1-1 do Begin
        i1 := (i4 shl k4) + i5;
        i2 := i1 + k1;
        i3 := i5 shl k3;
        t1 := fr[i2].re * ft[i3].re - fr[i2].im * ft[i3].im;
        t2 := fr[i2].re * ft[i3].im + fr[i2].im * ft[i3].re;
        t1 := t1 div fb;
        t2 := t2 div fb;
        fr[i2].re := fr[i1].re - t1;
        fr[i2].im := fr[i1].im - t2;
        fr[i1].re := fr[i1].re + t1;
        fr[i1].im := fr[i1].im + t2;
      End;
    End;
    k1 := k1 shl 1;
    k2 := k2 shr 1;
    k3 := k3 - 1;
    k4 := k4 + 1;
  End;
  // Normierung
  For i1 := 0 to n-1 do Begin
    fr[i1].re := fr[i1].re div n;
    fr[i1].im := fr[i1].im div n;
  End;
  // Powerspektrum
  For i1 := 0 to n-1 do Begin
    t1 := fr[i1].re * fr[i1].re + fr[i1].im * fr[i1].im;
    pw[i1] := IntSqrt(t1);
  End;
End;

   

Floatingpoint FFT

Round(Float FFT)

Integer FFT

  

n

 Input

Re-Teil

Im-Teil

Power

Re-Teil

Im-Teil

Power

Re

Im

Power

 Abs Fehler
0 1023 127,875 0,000 127,875 128 0 128 127 0 127 -0,875
1 1023 119,440 36,232 124,814 119 36 125 119 36 124 -0,814
2 1023 96,343 64,375 115,871 96 64 116 96 64 115 -0,871
3 1023 64,547 78,651 101,746 65 79 102 64 78 100 -1,746
4 0 31,969 77,179 83,538 32 77 84 31 77 83 -0,538
5 0 6,141 62,353 62,655 6 62 63 6 62 62 -0,655
6 0 -7,938 39,907 40,689 -8 40 41 -7 39 39 -1,689
7 0 -9,091 17,007 19,284 -9 17 19 -9 16 18 -1,284
8 0 0,000 0,000 0,000 0 0 0 0 0 0 0,000
9 0 13,958 -7,460 15,826 14 -7 16 13 -7 14 -1,826
10 0 26,665 -5,304 27,187 27 -5 27 26 -5 26 -1,187
11 0 33,328 3,283 33,490 33 3 33 33 3 33 -0,490
12 0 31,969 13,242 34,603 32 13 35 31 13 33 -1,603
13 0 23,858 19,580 30,864 24 20 31 23 19 29 -1,864
14 0 12,805 19,164 23,048 13 19 23 12 19 22 -1,048
15 0 3,569 11,764 12,293 4 12 12 3 11 11 -1,293
16 0 0,000 0,000 0,000 0 0 0 0 0 0 0,000
17 0 3,569 -11,764 12,293 4 -12 12 3 -11 11 -1,293
18 0 12,805 -19,164 23,048 13 -19 23 12 -19 22 -1,048
19 0 23,858 -19,580 30,864 24 -20 31 23 -19 29 -1,864
20 0 31,969 -13,242 34,603 32 -13 35 31 -13 33 -1,603
21 0 33,328 -3,283 33,490 33 -3 33 33 -3 33 -0,490
22 0 26,665 5,304 27,187 27 5 27 26 5 26 -1,187
23 0 13,958 7,460 15,826 14 7 16 13 7 14 -1,826
24 0 0,000 0,000 0,000 0 0 0 0 0 0 0,000
25 0 -9,091 -17,007 19,284 -9 -17 19 -9 -16 18 -1,284
26 0 -7,938 -39,907 40,689 -8 -40 41 -7 -39 39 -1,689
27 0 6,141 -62,353 62,655 6 -62 63 6 -62 62 -0,665
28 0 31,969 -77,179 83,538 32 -77 84 31 -77 83 -0,538
29 0 64,547 -78,651 101,746 65 -79 102 64 -78 100 -1,746
30 0 96,343 -64,375 115,871 96 -64 116 96 -64 115 -0,871
31 0 119,440 -36,232 124,814 119 -36 125 119 -36 124 -0,814

Hieraus läßt sich jetzt die Assembler Version entwickeln.