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?
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.
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.
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.
