====== std.fft — schnelle Fourier-Transformation ====== Radix-2 für Zweierpotenzen, Bluestein für jede andere Länge; dazu Faltung und Korrelation über die Frequenzebene, Fensterfunktionen und Spektren. → [[lyx_-_programmiersprache:units:mathematik|Mathematik-Units]] · [[lyx_-_programmiersprache:units:convolve|std.convolve]] · [[lyx_-_programmiersprache:units:wavelet|std.wavelet]] · [[lyx_-_programmiersprache:units:complex|std.complex]]


import std.fft;

Schicht 0, hängt an ''std.math'', ''std.math.constants'' und ''std.alloc''. 21 Funktionen. Alle Beispiele mit ''lyxc 1.0.21A'' übersetzt und ausgeführt. ---- ===== Darstellung ===== Ein komplexes Signal sind **zwei** ''f64''-Felder gleicher Länge — Realteil und Imaginärteil, getrennt. Die Transformation arbeitet **an Ort und Stelle**: nach ''FftForward'' stehen in denselben Feldern die Frequenzanteile, das Eingangssignal ist überschrieben.


var re: int64 := alloc(n * 8);
var im: int64 := alloc(n * 8);
// ... füllen ...
FftForward(re, im, n);      // re/im enthalten jetzt das Spektrum

Wer das Signal danach noch braucht, kopiert es vorher. ---- ===== Ein vollständiger Durchlauf =====


unit main;
import std.io;
import std.string;
import std.alloc;
import std.math;
import std.math.constants;
import std.fft;

fn ZF(t: pchar, v: f64): void { Print(t); PrintF64(v); }
fn ZI(t: pchar, v: int64): void { Print(t); PrintLn(IntToStr(v)); }

fn main(): int64 {
    var n: int64 := 16;
    var re: int64 := alloc(n * 8);
    var im: int64 := alloc(n * 8);

    // Sinus mit 3 Perioden über 16 Abtastwerte
    var i: int64 := 0;
    while (i < n) {
        var t: f64 := (i as f64) / (n as f64);
        var arg: f64 := 2.0 * PI * 3.0 * t;
        pokef64(re + i * 8, SinF64(arg));
        pokef64(im + i * 8, 0.0);
        i := i + 1;
    }

    var e1: f64 := 0.0;
    i := 0;
    while (i < n) { var v: f64 := peekf64(re + i*8); e1 := e1 + v * v; i := i + 1; }

    ZI("IsPowerOfTwo(16)   : ", FftIsPowerOfTwo(16) as int64);
    ZI("NextPowerOfTwo(17) : ", FftNextPowerOfTwo(17));

    FftForward(re, im, n);
    ZI("PeakBin            : ", FftPeakBin(re, im, n));
    ZF("BinFrequency(3, 16000 Hz): ", FftBinFrequency(3, n, 16000.0));

    var e2: f64 := 0.0;
    i := 0;
    while (i < n) {
        var r: f64 := peekf64(re + i*8);
        var m: f64 := peekf64(im + i*8);
        e2 := e2 + r * r + m * m;
        i := i + 1;
    }
    ZF("Energie Zeitbereich: ", e1);
    ZF("Energie Freq / n   : ", e2 / (n as f64));

    FftInverse(re, im, n);
    ZF("nach Hin+Rueck x0  : ", peekf64(re));
    ZF("nach Hin+Rueck x4  : ", peekf64(re + 32));
    return 0;
}


IsPowerOfTwo(16)   : 1
NextPowerOfTwo(17) : 32
PeakBin            : 3
BinFrequency(3, 16000 Hz): 3000.000000
Energie Zeitbereich: 8.000000
Energie Freq / n   : 8.000000
nach Hin+Rueck x0  : 0.000000
nach Hin+Rueck x4  : -0.999999

Drei Aussagen stecken darin: * **''PeakBin'' ist 3** — genau die Zahl der Perioden im Fenster. Bei 16 000 Hz Abtastrate entspricht Bin 3 einer Frequenz von 3000 Hz, und ''FftBinFrequency'' rechnet das aus. * **Die Energie stimmt überein** (8 = 8). Das ist der **Satz von Parseval**: die Energie im Zeitbereich gleicht der im Frequenzbereich, geteilt durch ''n''. Eine Prüfung, die **keine fremde Referenz braucht** — sie gilt aus sich heraus und deckt einen falschen Normierungsfaktor sofort auf. * **Hin und zurück stellt wieder her.** ''x4'' war ''−1'' und ist es nach beiden Schritten wieder (bis auf die letzte Stelle). ---- ===== Welches Verfahren ===== ^ Funktion ^ Länge ^ Aufwand ^ | ''FftRadix2'' | nur Zweierpotenzen | O(n log n) | | ''FftBluestein'' | **jede** Länge | O(n log n), größere Konstante | | ''FftDft'' | jede | O(n²) — zum Gegenprüfen | | ''FftForward'' / ''FftInverse'' | jede | wählt selbst | **''FftForward'' und ''FftInverse'' sind der normale Weg.** Sie nehmen Radix-2, wo es geht, und Bluestein sonst. Die Einzelverfahren braucht man nur, wenn man das Verhalten festlegen will — oder ''FftDft'' als langsame, aber unabhängige Gegenrechnung beim Prüfen. ''FftIsPowerOfTwo'' und ''FftNextPowerOfTwo'' helfen beim Auffüllen, wenn man Bluestein vermeiden will. ''FftForwardReal(data, n, outRe, outIm)'' nimmt ein rein reelles Signal und spart das Nullfeld für den Imaginärteil. ---- ===== Spektren auswerten ===== ^ Funktion ^ Liefert ^ | ''FftMagnitude(re, im, n, out)'' | Betrag je Bin | | ''FftPower(re, im, n, out)'' | Betragsquadrat | | ''FftPhase(re, im, n, out)'' | Phasenwinkel | | ''FftPeakBin(re, im, n)'' | Bin mit dem größten Betrag | | ''FftBinFrequency(bin, n, sampleRate)'' | Bin → Frequenz in Hz | ---- ===== Fenster =====


var win: int64 := alloc(n * 8);
FftWindowHann(win, n);
FftApplyWindow(data, win, n);

^ Fenster ^ Eigenschaft ^ | ''FftWindowHann'' | guter Allrounder | | ''FftWindowHamming'' | niedrigerer erster Nebenzipfel | | ''FftWindowBlackman'' | stärkste Dämpfung, breiterer Hauptzipfel | > **Ein Fenster ändert die Amplituden.** Wer nach dem Fenstern absolute Pegel ablesen will, muss durch die Fensterverstärkung teilen — ''FftWindowGain(window, n)'' liefert sie. Ohne diesen Schritt sind die Beträge systematisch zu klein. > > Gefenstert wird, weil die FFT das Fenster als **periodisch fortgesetzt** ansieht. Passt die Signalfrequenz nicht genau in das Fenster, entsteht am Rand ein Sprung, und dessen Energie verteilt sich über das ganze Spektrum („Leckeffekt"). Im Beispiel oben passt der Sinus mit genau 3 Perioden exakt hinein — deshalb ist dort kein Fenster nötig und der Peak ist scharf. ---- ===== Faltung und Korrelation ===== ^ Funktion ^ Zweck ^ | ''FftConvolveReal(a, na, b, nb, out)'' | lineare Faltung zweier reeller Folgen | | ''FftConvolveCircular(are, aim, bre, bim, n)'' | zyklische Faltung, an Ort und Stelle | | ''FftCorrelateReal(a, b, n, out)'' | Kreuzkorrelation | ''FftConvolveReal'' liefert ''na + nb − 1'' Werte — der Zielpuffer muss so groß sein. > **Der Unterschied zwischen linear und zyklisch ist kein Detail.** Die FFT faltet von Natur aus **zyklisch**: was rechts hinausläuft, kommt links wieder herein. Für eine lineare Faltung müssen beide Folgen vorher auf mindestens ''na + nb − 1'' mit Nullen aufgefüllt werden — ''FftConvolveReal'' erledigt das, ''FftConvolveCircular'' ausdrücklich nicht. Wer die zyklische Fassung für Signalfilterung nimmt, bekommt am Anfang des Ergebnisses das Ende des Signals eingemischt. Für Filter mit Randbehandlung, 2D und Overlap-Add ist [[lyx_-_programmiersprache:units:convolve|std.convolve]] die passendere Unit — sie setzt auf ''std.fft'' auf. ---- ===== Offener Punkt ===== > **Die inneren Schleifen summieren noch naiv.** ''std.fft'' ist an einer Stelle nicht auf kompensiertes Summieren umgestellt (''std.convolve'' an acht). Bei sehr langen Signalen mit stark verschiedenen Größenordnungen lohnt eine Gegenprobe über [[lyx_-_programmiersprache:units:sum|std.sum]]. Für übliche Längen und Wertebereiche fällt es nicht ins Gewicht — die Parseval-Prüfung oben stimmt exakt. ---- Letzte Aktualisierung: 2026-08-16 · alle Beispiele mit ''lyxc 1.0.21A'' übersetzt und ausgeführt; Peak-Bin, Frequenzzuordnung und Rücktransformation gegen eine unabhängige Referenzrechnung geprüft, die Energiegleichheit als referenzfreie Invariante nachgerechnet.