|
Übungen zu Mikrocontroller Anwendungen
- Bei der Analsyse von Zeitsignalen spielt die Fouriertransformation eine wichtige Rolle.
- Ihre schnelle Variante, die FFT, läßt sich auch mit den begrenzten Resourcen eines Mikrocontrollers umsetzen.
- Zudem wird in der Vorlesung Regelungssysteme bald eine weitere Transformationsvorschrift, die Laplace-Transformation, eingeführt werden.
- Deshalb bietet es sich an in dem für dieses Semester angesetzten Übungsteil zu Mikrocontroller-Anwendungen, die bereits behandelte DFT wieder aufzufrischen und zu ergänzen, um sich im Anschluß dann mit der FFT zu befassen.
|
- Im letzten Semester haben Sie die Digitale Fouriertransformation kennengelernt.
- Die Transformationsvorschrift für eine Meßreihe mit N Messungen und einer Samplingrate von fs, bzw. T=1/fs, lautete:
- n: Laufindex über sämtliche Frequenzanteile.
- k: Laufindex des Zeitsignals über sämtliche Zeitpunkte.
|
Bild 0-1: Vorschfift zur Transformation eines Zeitverlaufs in den Frequenzbereich mit Hilfe der Digitalen Fouriertransformation.
- Trennt man Real- und Imaginärteil der Formel mit Hilfe der Euler-Beziehung, ergibt sich:
|
Bild 0-2: Aufteilung in Real- und Imaginärteil
- Bisher wurde der inverse Vorgang, die Rücktransformation, nicht behandelt.
- Die Formel hierzu lautet:
|
Bild 0-3: Inverse DFT
Übung 1
- Verwenden Sie die Euler-Beziehung, um anzugeben, wie die Rücktransformation u(kT) aus dem Real- und Imaginärteil der Vorwärtstransformation erfolgen kann.
- Erzeugen Sie mit C++ eine Sinusschwingung der Amplitude +/-1 mit einer Frequenz von 100Hz und einer Samplingrate von 44100Hz mit der Dauer von 0,1 Sekunde.
- Transformieren Sie mit C++ die entstandenen Samples mit Hilfe der DFT in den Frequenzbereich.
- Speichern Sie die entstandenen Daten in eine Datei und visualisieren Sie diese mit Hilfe von Scilab.
- Führen Sie nun ebenfalls in C++ die inverse DFT auf das Transformationsergebnis aus. Vergleichen Sie das Ergenis mit der Ausgangskurve, um die Richtigkeit Ihres Algorithmus zu überprüfen.
- Hinweis - Kompilieren eines C++-Quelltextes in der Konsole: g++ -o test test.cpp
- Hinweis - Daten in einer Datei speichern in C++: siehe vorgegebenes Material in Übung 1, Regelungssysteme, Kapitel 1.3
- Hinweis - Daten in Scilab von einer Datei laden und visualisieren: siehe vorgegebenes Material in Übung 1, Regelungssysteme, Kapitel 1.3
|
FFT - Fast-Fourier-Transformation
- Bei der DFT treten redundante Berechnungen auf.
- Diese Redundanzen werden besonders groß, wenn die Anzahl N der zu transformierenden Elemente der Zeitreihe einer Potenz von 2 entspricht, also N={2,4,8,...,4096,8192,usw.)
- Nutzt man dies aus, so gelingt es den Rechenaufwand von der Ordnung N2, d.h. O(N2) auf die Ordnung NlogN, d.h. O(NlogN), reduzieren.
- Anstatt also der Rechenaufwand quadratisch mit der Anzahl N der Elemente der zu transformierenden Zeitreihe u(ti) ansteigt, steigt er nur noch mit NlogN an.
- Im ersten Moment mag diese nicht gerade ungeheuerlich groß erscheinen, doch ein konkretes Beispiel soll diesen Eindruck wiederlegen:
- Angenommen eine Audioaufnahme von 10 Sekunden Dauer, die in CD-Qualität aufgezeichnet wurde soll transformiert werden.
- Die entsprechende Samplingrate wäre 44100 Samples pro Sekunde und damit 441000.
- Die nächste Zweierpotenz läge bei 219 = 524288.
- Berechnet man als Grundeinheit 1ms, so bräuchte die DFT ca. 75 Stunden, die FFT jedoch nur ca. 90 Minuten.
- Grundlage der FFT bildet die Tatsache, dass es möglich ist, wenn man getrennt voneinander die DFT aller N/2 geraden Samples und die aller N/2 ungeraden bereits vorliegen hat, die aller N Samples auf folgende Weise bekommt:
- Hinweis: n=0..N-1, die Transformierten mit geraden Indices (0,2,4,6,..) sind beschrieben 2n, die mit Ungeraden durch 2n+1.
|
Bild 0-4: Grundlage der FFT ist die Möglichkeit aus zwei schon berechneten Hälften einer DFT eine ganze zu erzeugen.
Vorüberlegungen zu einer Implementierung des FFT-Algorithmus
- Ein Algorithmus für die FFT wird bei Transformierten der Länge 1 beginnen und sie zu solchen der Länge 2 zusammenfügen.
- Dann kämen die Längen 4,8,16, usw. dran.
- Um sich Klarheit darüber zu verschaffen, was zu was zusammengeführt wird, ist dies für die Länge N=8 illustriert:
- Die Darstellung erfolgt rückwärts beginnend mit dem letzten Schritt.
- Davor werden auch alle Elemente mit geraden Indices wiederum in Elemente mit geraden und ungeraden Indices aufgeteilt gedacht und das gleiche für die ungeraden.
- Es zeichnet sich ab, dass die Schrittweite von einem Element mit jeder früheren Rechenebene verdoppelt, während sich die Längen halbieren.
- Auf jeder Ebene muß die Berechnung für jede mögliche Verschiebung innerhalb der Schrittweite durchgeführt werden, also sozusagen für jede Phase der Schrittperiode.
- Für die Schrittweite 4 in der Ebene0 z.B. für die Phasen 0,1,2,3.
|
01234567
gugugugu Ebene2 - letzter Schritt: gerade und ungerade Indices der Länge N/2
g.u.g.u. Ebene1 - vorletzter Schritt: Länge N/4
.g.u.g.u
g...u... Ebene0 - Länge N/8
.g...u..
..g...u.
...g...u
Code 0-1: Berechnungsphasen einer FFT für N=8.
Aufgabe
- Um den FFT Algorithmus implementieren zu können, muß bestimmt werden, wie die DFT bei einer Länge von N=1 aussieht, da dies die Startlänge auf Ebene 0 ist.
- Weisen Sie nach, dass folgendes gilt:
|
Bild 0-5: DFT für N=1
Übung2
- Als weitere Vorbereitung zur Implementierung eines FFT-Algorithmus soll zunächst nur der letzte Schritt auf der letzten Ebene funktionierend für das Beispiel aus Übung 1 implementiert werden und mit dem Ergebnis aus Übung 1 vergleichen werden.
- Schreiben Sie also ein Programm, in dem zunächste die Transformierte von G für alle geraden Indices von u bestimmt wird (Ggerade_re, Ggerade_im)
- ..und das gleiche für die ungeraden Indices durchgeführt wird (Gungerade_re, Gungerade_re).
- Schließlich sollen beide Teile nach obiger Regel für FFT zu einer Transformierten von 0..N-1 zusammengeführt werden und dessen Spektrum gespeichert werden.
|
Algorithmus für die Fast-Fourier-Transformation (FFT)
- Das nachfolgende Programm setzt die FFT für das Beispiel aus Übung 1 so um, wie das zuvor beschrieben wurde.
- Beginnend mit DFTs der Länge 1, verdoppelt sich mit jeder Ebene die Länge der gebildeten Transformierten.
- Um für N eine Zweierpotenz zu haben, wurde die Länge auf 4096 Samples gesetzt.
- Um die Ergebnisse von FFT und DFT vergleichen zu können, indem man die Verläufe der Spektren superpositioniert, muß auch die DFT-Berechnung für 4096 Samples durchgeführt werden.
- Häufig wird in Pseudocode eine rekursive Implementierung für die FFT angegeben.
- Dies ist jedoch für praktische Anwendungen mit größeren N nicht sinnvoll, da für jede Ebene noch einmal der komplette Speicher für N Samples bis zum Ende des Algorithmus behalten wird.
|
#include<iostream>
#include<fstream>
#include<math.h>
#include<stdlib.h>
using namespace std;
void fft(double* Ghilf_re, double* Ghilf_im, double* G_re, double* G_im, double* u, int N)
{
//1. Packe in G_re=u und G_im=0 hinein, was einer DFT der Länge 1 entspricht:
for(int i=0;i<N;i++)
{
G_re[i]=u[i];
G_im[i]=0.0;
}
//2. Gehe alle Ebenen durch.
//Die letzte Hat eine Verschränkung und verschränkt Gerade und Ungerade der Länge N/2
//Die vorlestzte hat 2 der Länge N/4
//4 N/8
//8 N/16
//...
//N/2 1
//Berechnung der Anzahl der Ebenen:
//Bei 8=2hoch3 gibt es 3 Ebenen, delta beginnt mit N/2 und phase von 0..N/2-1
//Einfacher: man geht alle deltas und periodenlängen in einer while-Schleife durch:
int laenge = 1; //Länge einer der beiden Zeitreihen, die verbunden werden sollen
int delta = N/2; //Schrittweite von der geraden zur ungeraden Reihe, der doppelte
//Wert wäre der Abstand von einem Wert zum nächsten der gleichen Reihe.
while(laenge<N) //Schleife, die alle Ebenen durchgeht
{
for(int phase=0;phase<delta;phase++) //Es existiert jede Verrückung um 1 als weitere zu
//berechnende Verschrränkung auf der gleichen Ebene
{
//FFT auf der aktuellen Ebene
//erste Hälfte 0..laenge
for(int n=0;n<laenge;n++)
{
//Berechnung der Indices für gerade und ungerade Elemente
//der gerade miteinander zu verknüpfenden Elemente:
Ghilf_re[n*delta+phase] = G_re[n*delta*2+phase] //gerader Anteil
+cos(M_PI*(double)n/(double)laenge)*G_re[n*delta*2+delta+phase]
+sin(M_PI*(double)n/(double)laenge)*G_im[n*delta*2+delta+phase];
Ghilf_im[n*delta+phase] = G_im[n*delta*2+phase+0] //gerader Anteil
+cos(M_PI*(double)n/(double)laenge)*G_im[n*delta*2+delta+phase]
-sin(M_PI*(double)n/(double)laenge)*G_re[n*delta*2+delta+phase];
}
//zweite Hälfte laenge..2*laenge
for(int n=0;n<laenge;n++)
{
Ghilf_re[(n+laenge)*delta+phase] = G_re[n*delta*2+phase]
-cos(M_PI*(double)n/(double)laenge)*G_re[n*delta*2+delta+phase]
-sin(M_PI*(double)n/(double)laenge)*G_im[n*delta*2+delta+phase];
Ghilf_im[(n+laenge)*delta+phase] = G_im[n*delta*2+phase]
-cos(M_PI*(double)n/(double)laenge)*G_im[n*delta*2+delta+phase]
+sin(M_PI*(double)n/(double)laenge)*G_re[n*delta*2+delta+phase];
}
}
//Zurückkopieren von Ghilf nach G vorbereitend für die nächste Ebene
for(int i=0;i<N;i++)
{
G_re[i]=Ghilf_re[i];
G_im[i]=Ghilf_im[i];
}
laenge*=2;
delta/=2;
}
}
int main()
{
int sample_rate = 44100;
double f=100.0;
// int N=4410;
int N=4096; //2hoch12
int NH=N/2;
double* u = new double[N];
double* u_rueck = new double[N];
double* zeit = new double[N];
double* trans_re = new double[N];
double* trans_im = new double[N];
double* G_re = new double[N];
double* G_im = new double[N];
double* Ghilf_re = new double[N];
double* Ghilf_im = new double[N];
double t = 0.0;
double dt = 1.0/(double)sample_rate;
for(int i=0;i<N;i++)
{
u[i] = sin(2.0*M_PI*f*t);
zeit[i]=t;
t+=dt;
}
ofstream uin("uin.txt");
for(int i=0;i<N;i++)
uin<<t<<" "<<u[i]<<endl;
uin.close();
fft(Ghilf_re,Ghilf_im,G_re,G_im,u,N);
ofstream uspek("uspek_fft003.txt");
for(int i=0;i<N;i++)
{
double fi = (double)i*(double)sample_rate/(double)N;
uspek<<fi<<" "<<sqrt(G_re[i]*G_re[i]+G_im[i]*G_im[i])<<endl;
}
uspek.close();
}
Code 0-2: Nicht rekursive Implementierung eines FFT-Algorithmus
|