Einfache FFT Implementierung Danielson-Lanzcos

[Martin R (Gast)]

Hallo,
Auf meinem TMS320VC5505 läuft bereits die FFT.
Zur Verständnis habe ich in C noch eine weitere Routine angewandt wie 
hier beschrieben:
http://www.codeproject.com/KB/recipes/howtofft.aspx
http://www.mathcs.org/java/programs/FFT/FFTInfo/c12-2.pdf
der bitrevers Teil läuft einwandfrei.
Leider bekomme ich nach der Danielson-Lanzcos Routine ein falsches 
Ergebnis:
Als Eingangsfolge nehme ich:
float data_buf[16] ={

1000,  /* 0 */
0,  /* 1 */
200,  /* 2 */
0,  /* 3 */
1000,  /* 4 */
0,  /* 5 */
200,  /* 6 */
0,  /* 7 */
1000,  /* 8 */
0,  /* 9 */
200,  /* 10 */
0,  /* 11 */
1000,  /* 12 */
0,  /* 13 */
300,  /* 14 */
0,  /* 15 */
};

Ergebnis in Matlab:
4900 + 0i
70.7107 + 70.7107i
0 + 100i
-70.7107 + 70.7107i
3100 + 0i
-70.7107 - 70.7107i
0 - 100i
70.7107 - 70.7107i

Ergebnis vom DSP:
1000 + 0i
3900 - 1012.132i
-2121.322 - 100i
800 + 5878.679i
212.132 + 0i
2500 - 587.868i
212.132 + 100i
800 + 1012.132i

Das kann doch nicht sein, oder? Danielson-Lanczos hat doch da sicherlich 
nichts falsches gemacht ;-)

Gruß

[Xeraniad X. (xeraniad)]

Bei der Ausgabe fällt z. B. die unveränderte Array-Zelle data_buf [0] 
auf. Ist da bei der Indizierung ein Offset +1 reingerutscht? 
Möglicherweise gibt es auch ein Problem bei der Initialisierung des 
Array. Wenn Du den verwendeten C-Code postest, werde ich ihn mal auf nem 
PC (Linux, GNU g++) nachbilden & mit dem gegebenen Input-Array die 
Reproduktion des Problems versuchen. Dann kann evtl. entschieden werden, 
ob ein (Tipp?-)Fehler im Algorithmus oder ob ein Hardware-Problem 
vorliegt.
Es würde mich auch interessieren, wie data_buf in einen Array von 16 
komplexen Zahlen kopiert wird. Es wird ja ein in-place-Algorithmus sein.
Oder war es die Absicht, 8 reelle Abtastwerte zu übergeben, wobei die 
Realteile bei geraden und die Imaginärteile (= 0) bei ungeraden Indices 
gespeichert werden?

[Martin R (Gast)]

Habe das gerade noch mal in CodeBlocks nachgebildet und in den Anhang 
verpackt.
Bei genauerem Hinsehen ist mir auch aufgefallen dass der Code der 
Lanzcos Routine (nicht der bitreverse Teil) sich bei den beiden von mir 
angegebenen Quellen leicht unterscheidet.
Es ist absicht, dass der Imaginärteil der Quelle 0 ist. Es sollen nur 
reelle Quellen verarbeitet werden.

[Xeraniad X. (xeraniad)]

Hallo

Danke fuer den Quell-Code. Es gelang, die Prozedur "four1" aufzurufen. 
Die Erkenntnisse sind wie folgt.
° Ich habe nur den Algorithmus "four1" in 
http://www.mathcs.org/java/programs/FFT/FFTInfo/c12-2.pdf betrachtet.
° Der Algorithmus startet die Array-Indices bei 1, nicht 0 wie bei C 
üblich.
Die Zelle mit Index 0 bleibt unverwendet, dafür muss die Kapazität 17 
für data_buf deklariert werden, weil der Algorithmus auf data_buf [16] 
zugreift. Beträgt die Kapazität nur 16, erfolgen (auch schreibend) 
Zugriffe auf andere oder nicht deklarierte Variablen, mit den 
entsprechenden Konsequenzen.
° isign = -1 heisst FFT, für inverse FFT muss isign = +1 gewählt werden.
° Im Beispiel waren die acht Abtastwerte [10, 2, 10, 2, 10, 2, 10, 3] 
gegeben.
MathLab erwartet die Abtastwerte in dieser Form, nicht mit den 
Imaginär-Nullen dazwischen.
Der Algorithmus "four1" erwartet die Abtastwerte (bzw. seinen input) in 
der Form
 ₀   ₁  ₂   ₃  ₄    ₅  ₆   ₇  ₈    ₉  ₁₀  ₁₁ ₁₂   ₁₃ ₁₄   ₁₅ ₁₆
[?, 10, 0,  2, 0,  10, 0,  2, 0,  10, 0,  2,  0,  10, 0,   3, 0]
übergeben.
° Die acht komplexen Fourier-Koeffizienten für die gegebenen Abtastwerte 
sind
49, exp(j·¼·π), j, exp(j·¾·π), 31, exp(-j·¾·π), -j, exp(-j·¼·π).
° Der Algorithmus "four1" gibt die Resultat-Werte in der Form
 ₀   ₁  ₂           ₃          ₄   ₅  ₆           ₇         ₈
[?, 49, 0,   0.707107, 0.707107,   0, 1,  -0.707107, 0.707107,
 ₉ ₁₀          ₁₁         ₁₂  ₁₃  ₁₄         ₁₅          ₁₆
31, 0,  -0.707107, -0.707107,  0, -1,  0.707107,  -0.707107]
zurück.
° Fazit: Der Algorithmus ist korrekt. Wegen Fortran-Array-Index-Style 
und Speicherung von Real -und Imaginärteilen im gleichen Array muss der 
Aufruf etwas exotisch erfolgen.

[Xeraniad X. (xeraniad)]

P. S.
° Zuvor hatte ich noch nicht bemerkt, dass das von Dir gepostete 
MathLab-Resultat für die Abtastwerte [1000, 200, 1000, 200, 1000, 200, 
1000, 300] selbstverständlich korrekt ist.
° Der von Dir erwähnte Unterschied zwischen den beiden Quellen (
http://www.codeproject.com/KB/recipes/howtofft.aspx und 
http://www.mathcs.org/java/programs/FFT/FFTInfo/c12-2.pdf) bezieht sich 
genau auf den Array-Start-Index.
Die erste Quelle wurde von der zweiten kopiert & modifiziert, sodass 
Array-Indices C-konform bei 0 starten (jedenfalls ab der zweiten 
Bit-Reverse-Version). Dann befinden sich die Realteile bei geraden und 
Imaginärteile bei ungeraden Indices, und es kann data_buf [16] 
deklariert werden. (Nebenbei fragt man sich, weshalb die Prozedur "FFT" 
einen Rückgabewert "float" haben soll.)
Die zweite Quelle (und "four1" im geposteten C-Programm) indizieren 
Arrays von 1 an (z. B. wie FORTRAN oder PASCAL), dann befinden sich 
Realteile bei ungeraden Indices und Imaginärteile bei geraden Indices. 
Ausserdem wird Speicherplatz für eine Realzahl verschwendet.
Jetzt passt alles zusammen.