Forum: Digitale Signalverarbeitung / DSP / Machine Learning FFT mit reelwertigen Einganswerten mit halb so großer FFT


von Andreas W. (andy_w)


Lesenswert?

Hallo,
folgendes Problem: ich will eine FFT mit 4096 Werten in einem LPC1769 
machen. Der hat 2*32K RAM für die FFT braucht man aber bereits 
4096*2*sizeof(float) = 32k Speicher. Da in dem Projekt auch noch etwas 
sonstiger Speicher gebraucht wird, reicht das RAM erwartungsgemäß nicht 
aus.

Für die FFT muß man ja 4096 komplexe Zahlen im Speicher vorsehen, also 2 
float zu je 4 Bit. Da bei mir aber die Einganswerte alle reelwertig sind 
und damit alle imaginären Anteile Null sind, gibt es Tricks, das Ganze 
mit einer FFT mit 2048 Werten durchzuführen. Es gibt z.B. ein Script 
hier:

www.zhaw.ch/~rumc/dsv1/unterlagen/dsv1kap3dftfft.pdf

Auf den letzten beiden Seiten wird es beschrieben. Die reelen 4096 
Eingangswerte werden abwechselnd auf die Real- und Imaginärteile der 
2048-FFT verteilt, dann die FFT durchgeführt und anschließend daraus die 
4096 komplexen Ergebnisse berechnet. Das geht, weil im Ergebnis Werte 
mit den Indizes i und 4096-i konjugiert komplex zueinander sind. Leider 
ist die Beschreibung ausgerechnet in dem Punkt in der Doku fehlerhaft, 
die Indizes gehen dort durcheinander. Die Eingangswerte gehen von 0 bis 
M-1, die FFT-werte von 0 bis M/2-1, soweit ist das korrekt. Aber später 
im Ergebnis werden Indizes von 0 bis M-1 für die Ergebniswerte (das ist 
richtig, es gibt ja 4096 Stück) auch in die Ausgangswerte der FFT 
eingesetzt und das geht nicht, denn die FFT liefert ja nur 2048 Werte, 
also Indizes von 0 bis M/2-1.

Die einzige andere Doku, die diese Thema anschneidet, fand ich über 
Google bei der FH München und TH Wien (die Dokumente sind identisch), 
aber genau bei dem Punkt geraten da wieder die Indizes durcheinander, 
man muß dann wieder FFT-Ergebnisse mit Indizes bis 4095 adressieren...

Ich brauche leider die 4096-FFT, reserviere ich Speicher für eine 
2048-FFT, reicht es aus, jetzt brauche ich "nur" noch eine fehlerfreie 
Doku, die das Verfahren der Schachtelung korrekt beschreibt. Noch besser 
wäre natürlich ein Codebeispiel ;-) Eine Doku würde mir aber reichen. 
Die Rechengeschwindigkeit ist nicht sooo kritisch, ich brauche nur eine 
Lösung, die mit halb soviel Speicher auskommt als die Standardlösung mit 
der normalen FFT.

Das ganze soll für einen digitalen Audioverstärker verwendet werden, 
dort kommen FIR-Filter mit 4095 Koeffizienten (symmetrisch, daher 
lineare Phase) zum Einsatz. Die Filter selber laufen in FPGAs, die 
Koeffizienten dafür werden aber vom LPC1769 berechnet und einmal in die 
FPGAs geladen. Ich berechne einen Sollfrequenzgang und mithilfe einer 
IFFT daraus die Filterkoeffizienten (die müssen noch mit einer 
Fensterfunktion versehen werden). Die IFFT ist ja identisch zur FFT, sie 
unterscheidet sich nur in einem konstanten Faktor für alle 
Ergebniswerte. Bei mir ist sogar der reele Sollfrequenzgang symmetrisch 
bzgl. i und 4096-i, so daß die IFFT nicht nur konjugiert komplexe Werte 
liefert, sondern sogar rein reele, sonst würde das FIR-Filer auch nicht 
funktionieren. Eine Simulation ergibt auch schon gute Ergebnisse. Da die 
FIR-Filter bis in den Baßbereich wirksam sein sollen, müssen die schon 
4095 Koeffzienten haben, sonst kann man den unteren Baßbereich nicht 
mehr individuell beeinflussen.

Gruß
Andy_W

von MK (Gast)


Lesenswert?

Hi,
guck dir das mal an (v.A. ab Seite 47):

http://www.dss.tf.uni-kiel.de/teaching/lecture_digital_signal_processing/lecture_slides/dig_sig_03_spektren.pdf


Scheint sehr detailliert beschrieben zu sein.

lg,
MK

von MK (Gast)


Lesenswert?

... äh Seite 44 ntl.

von frame (Gast)


Lesenswert?

Schau dir doch mal den CMSIS DSP code an.

Ich nutze ihn in Form der STM32F4xx_DSP_StdPeriph_Lib. Der DSP-Teil 
kommt von ARM, sollte in der NXP-Variante also identisch sein.

Dort sind die FFT-Routinen übrigens in 3 Varianten enthalten, für float, 
16 Bit integer und 32 Bit integer. float würde ich für den lpc1769 ohne 
FPU nicht empfehlen - und dann noch 4096 Punkte...
Mein STM32F4 braucht bei 168MHz etwa 3ms für eine 1024 Punkte - FFT, bei 
maximaler Optimierung und mit FPU.

von Andreas W. (andy_w)


Angehängte Dateien:

Lesenswert?

Hallo,
das andere Dokument hatte ich auch schon gefunden, bei der z.B. hat mich 
Google zur FH München geführt. Das ist das zweite, das ich auch schon 
erwähnt hatte. Auch da geraten die Indizes etwas durcheinander, vor 
allem wird dazu kein einziges Wort gesagt, wie das zu interpretieren 
ist.

Aber ich hatte mir genau diese Vorlage vorgenommen und danach die Sacher 
selber hergeleitet, vielleicht zur Info ein Scan vom Schmierblatt mit 
meinen Notizen...

Daraus habe ich dann auch entsprechenden Code gemacht - und er 
funktioniert.
1
//******************************************************************************
2
//
3
//  real FFT/IFFT
4
//
5
//  nur reele, symmetrische Eingabewerte
6
//******************************************************************************
7
void fft_r(sint32 dir, sint32 m, float *xr, float *xi)
8
{
9
  sint32  n, i;
10
  float   f1r, f1i, f2r, f2i;
11
12
  n = 1 << m;
13
  for (i = 0; i < n/2; i++)
14
  {
15
    xr[i] = xr[i+i];
16
    xi[i] = xr[i+i+1];
17
  }
18
  fft(dir, m-1, xr, xi);
19
  xr[n/2] = xr[0];
20
  xi[n/2] = xi[0];
21
  for (i = 0; i < n/2; i++)
22
  {
23
    f1r = 0.5*(xr[i] + xr[n/2-i]);
24
    f1i = 0.5*(xi[i] - xi[n/2-i]);
25
    f2r = 0.5*(xi[i] + xi[n/2-i]);
26
    f2i = 0.5*(xr[n/2-i] - xr[i]);
27
    xr[n/2+i] = f1r + f2r*cos(2*PI/n*i) - f2i*sin(2*PI/n*i);
28
  }
29
  for (i = 0; i < n/2; i++)
30
  {
31
    xr[i] = xr[n/2+i];
32
  }
33
}

Der Code ist noch vorläufig, Hauptsache, der funktioniert erst einmal 
vom Prinzip. xr[] und xi[] sind die Arrays für die Real- und 
Imaginärteile der FFT. Beim Aufruf stehen die 4096 reellen Werte in 
xr[]. Da es sich um gespiegelte Werte handelt, also xr[i] = xr[4096-i], 
kann man später auch nur die Hälfte als Eingangswerte nehmen und den 
Rest entsprechend umkopieren. Diese 4096 Werte werden auf die ersten 
2048 Zellen von xr[] und xi[] verteilt. 2^m ist die Größe der FFT, m ist 
bei mir also 12 und n daher 4096. dir ist 1 für eine FFT und -1 für eine 
IFFT, der Unterschied ist lediglich ein konstanter Faktor nach der FFT. 
Die Parameter für die aufgerufene fft() sind die gleichen wie für diese 
fft_r(), es wird aber fft mit m-1 aufgerufen, also eine 2048-FFT.

Danach müssen die Ergebnisdaten aufbereitet werden. Zuerst trennt man 
das Ergebnis in zwei Teilergebnisse f1 und f2 (jeweils mit Real- und 
Imaginärteil, also f1r, f1i, f2r und f2i) auf, genau so, wie ein paar 
Seiten zuvor beschrieben für eine FFT mit zwei getrennten Datensätzen 
gemeinsam (das Verfahren war verständlich, da es immer gleich viel 
indizes blieben). Diese beiden Teilergebnisse werden dann zu einem 
großen zusammengefaßt (die Zeile nach Berechnung von r1r, r1i, r2r und 
r2i mit cos und sin). Da mein Ergebnis auch nur reelle Komponenten 
enthält, berechne ich auch nur den Realteil, eine weitere Zeile wäre für 
den Imaginärteil nötig, der entsteht, wenn die reellen Einganswerte 
nicht gespiegelt sind.

Meine Herleitung hat ergeben, daß beim Berechnen des Ergebnisses nur die 
untere Hälfte mit Index 0-2047 herauskommt, die obere Hälfte ist 
konjugiert komplexe gespiegelt, also ff[4096-i] = konjugiert(ff[i]). 
Davon steht in den Dokumenten z.B. kein Wort, das wird wohl nur mündlich 
in den Vorlesungen gesagt... xr[2048] und xi[2048] fehlt bei den 
Eingangsdaten, da 4096-2048 = 2048 ist, da sich die Eingangsdaten der 
2048-FFT wie sonst auch aber unendlich wiederholen (es wird ja 
angenommen, daß die periodisch sind), sind die identisch mit xr[0] und 
xi[0].

Das Ergebnis von ff[2048] fehlt aus dem gleichen Grund. Ich brauche das 
Element bei meiner Anwendung nicht, es ist für ein FIR-Filter mit 
Koeffizienten -2047...+2047 (gefenstert). Ich nehme mal an, daß es auch 
identisch mit ff[0] sein könnte, evtl. aber auch komplett Null. Auf 
jeden Fall ist der Imaginärteil Null, denn ff[2048] muß ja konjugiert 
komplex zu sich selbst sein.

Der Code ist schon ausprobiert und er funktioniert, noch auf dem PC 
unter Visual Studio. Nun kann ich den Code noch anpassen. Wahrscheinlich 
werde ich bei float bleiben (so bin ich sicher vor Overflow), die 
Rechenzeit ist nicht so kritisch, die Filterkoeffizienten werden ja 
sozusagen "offline" berechnet und dann in FPGAs geladen, die schuften 
sich dann mit dem FIR-Filter ab... Die FIR-Filter arbeiten dann mit 
Festkomma, intern mit 27 Bit (Spartan 3 Memoryblöcke mit je 2048*9 Bit).

Ich lege sehr viel Wert auf linaerphasige Filter bei Audio, besonders 
bei einem Equalizer, um Raumresonanzen abzuschwächen. Meine Aktivboxen 
mit noch analogen Besselfiltern, die auf möglichst gutmütige Phase 
optimiert sind, haben schon so manchen Besucher vom Sessel gerissen, 
selbst professionelle Abhörmonitore tun sich da schwer. Die 
Abhörmonitore V8 von KRK sind bei mir jedenfall nur gerade eben noch für 
die hinteren Surroundkanäle gut genug...

Gruß
Andy_W

von Anja (Gast)


Lesenswert?

Andreas W. schrieb:
> Die einzige andere Doku, die diese Thema anschneidet,

suche mal nach "real valued in place fft" z.b. nach rfft und rsfft von 
Sorensen. Ursprünglich in Fortran.

z.B. hier:
http://my.fit.edu/~vkepuska/ece5526/SPHINX/SphinxTrain/src/programs/cdcn_norm/rsfft.c

Gruß Anja

von Andreas W. (andy_w)


Lesenswert?

Hallo,
da ich inzwischen eine funktionierende Lösung habe (ich brauche auch für 
die Berechnungen vor und nach der FFT nur noch den halben Speicher, das 
Programm ist entsprechend umgeschrieben)m habe ich mir die letzte Lösung 
nur noch mal angesehen. Die ist vom Code her deutlich umfangreicher als 
z.B. meine FFT plus die Umrechnung. Welche Lösung schneller ist, habe 
ich aber nicht untersucht. Meine ist jedenfalls noch schnell genug 
(einige wenige 1/10s), solange kann ich noch warten, bis eine neue 
Klangeinstellung in die FPGAs geladen ist.

Zuerst hat es auf meinem LPC1769 fast 8s gedauert, bis ich ein Ergebnis 
sah, das Programm hat einen Wunschfrequenzgang berechnet, daraus die 
Filterkoeffizienten und daraus wiederum den Frequenzgang des FIR-Filters 
und auf dem Minigraphikdisplay des Evaluationboards dargestellt. Das war 
mir dann doch zu lange, bis ich den Übeltäter ausgemacht habe: die 
Berechnung des Frequenzgangs des FIR-Filters für 96 Frequenzen. Jedesmal 
muß dafür für alle 2048 Koeffizienten (FIR-Filter mit 4095 Stützstellen, 
aber symmetrische Koeffizienten) einen Cosinus berechnen, also insgesamt 
rund 200000-mal... Diese Berechnung fällt später in der Anwendung aber 
weg, die Filterkoeffizientenberechnung mit IFFT geht deutlich schneller.

Wie man den Frequenzgang eines FFT-Filters berechnet, habe ich hier 
gezeigt:
Beitrag "FIR Filter Betragsfunktion" letzte Formel
Da fällt mir aber kein schnellerer Algorithmus ein, auch die Berechnung 
des Cosinus aus vorhergehenden bringt leider nichts. Die Formel für 
Cosinus eines Vielfachen eines Winkels läuft auf eine binomische Formel 
hinaus, für größere Vielfache noch aufwendiger als der Cosinus selber 
und mit Additionsformel ist Sinus und Cosinus und zusätzlich sqrt nötig, 
bringt auch kaum Vorteile, dafür evtl. aufsummierende Rundungsehler. 
Aber die Funktion brauche ich bestenfalls zur Illustration, wenn der 
User den aktuellen Frequenzgang mal sehen will (der Verstärker bekommt 
ein LCD Graphikdisplay 800*480 mit Touchscreen), das kann auch im 
Hintergrund berechnet werden.

Gruß
Andy_W

von Detlef _. (detlef_a)


Lesenswert?

>>>>
Da fällt mir aber kein schnellerer Algorithmus ein, auch die Berechnung
des Cosinus aus vorhergehenden bringt leider nichts.
<<<<<<
Cosinus iterieren geht super:
Du willst z.B. Schrittweite von 0.123rad:
Die komplexe Zahl
z=e^j*0.123=0.9924+j*0.1227
dreht bei der Multiplikation eine andere komplexe Zahl um genau den 
Winkel.

(1+j*0)*z=0.9924+j*0.1227            -> cos(1*0.123)=0.9924
(0.9924+j*0.1227)*z=0.9699+j*0.2435  -> cos(2*0.123)=0.9699

Kostet 4 Multiplikationen und zwei Additionen.
Die Nummer läuft unter 'Cordic', geht super.

math rulez!
Cheers
Detlef

von Andreas W. (andy_w)


Lesenswert?

Hallo,
stimmt, inzwischen wurde ich auch fündig:

cos(n*a) = Tn(cos(a)) mit Tn(x) Tschebyscheff Polynom n-ter Ordnung 
(steht in Wikipedia). Hört sich erst einmal richtig fies an, aber es 
geht auch einfach, es gilt:
T(n+1)(x) = 2*x*Tn(x)-T(n-1)(x)
setzt man x = cos(a) kann man ebenfalls der Reihe nach iterieren:
x(1) = cos(a) muß einmal berechnet werden. x(0) = 1.
x(n+1) = 2*x(n)*cos(a)-x(n-1)
das kann man mit einer Multiplikation und 2 Additionen hinkriegen. x(m) 
ist dann cos(m*a), so daß man die Koeffizienten der Reihe nach aus den 
jeweils 2 vorhergehenden berechnen kann. Das Tschebyscheff Polynom 
selber muß man nie ermitteln und den Cosinus da durch jagen (auch mit 
Hornerschema dauert das für 2047-te Ordnung sehr lange...), man kann die 
Cosinuswerte selber iterieren.

Heute abend werde ich das mal im LPC1769 ausprobieren, sollte dann doch 
eine ganze Ecke schneller werden.

Gruß
Andy_W

von Andreas W. (andy_w)


Lesenswert?

Hallo,
inzwischen habe ich meinen Algorithmus mit den Tschebyscheffpolynomen 
implementiert. Er funktioniert, ich bekommen die gleichen Ergebnisse wie 
vorher. Und er läuft mehr als 8-mal so schnell für große FIR-Filter, bei 
mir eben mit 2047 Koeffizienten. In der Schleife eine Addition, eine 
Multiplikation, eine Subtraktion und zwei Moves (2 Variable umkopieren), 
mehr ist nicht nötig. Den Cosinus muß man nur einmal vor dem 
Schleifenstart berechnen.

Ja, Mathe rulez! Einige Mathelehrer in der Schule haben mich gefürchtet, 
sonst ist es doch eher umgekehrt ;-)

Gruß
Andy_W

von Detlef _. (detlef_a)


Lesenswert?

Yo, so geht das optimal.

Deine Gleichung:
x(n+1) = 2*x(n)*cos(a)-x(n-1)

ist die Differenzengleichung für ein rekursives, ungedämpftes System, 
das auf der Frequenz a schwingt. Diese Gleichung kann man auch 
heranziehen, um Frequenzen in Signalen zu bestimmen, siehe:
Beitrag "Frequenzlinie ausfiltern und deren Amplitude auswerten"

Dass man diese Gleichung aus Tschebyscheffpolynomen herleiten kann war 
mir neu, wieder was dabeigelernt, thx.

math rulez
Cheers
Detlef

von Estimator (Gast)


Lesenswert?


von Markus F. (Gast)


Lesenswert?

Nettes Doc! Sowas Ähnliches kann ich vielleicht demnächst auch 
gebrauchen!

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.