kramann.info
© Guido Kramann

Login: Passwort:










11 Ü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.
Vorschfift zur Transformation eines Zeitverlaufs in den Frequenzbereich mit Hilfe der Digitalen Fouriertransformation.

Bild 11-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:
Aufteilung in Real- und Imaginärteil

Bild 11-2: Aufteilung in Real- und Imaginärteil

  • Bisher wurde der inverse Vorgang, die Rücktransformation, nicht behandelt.
  • Die Formel hierzu lautet:
Inverse DFT

Bild 11-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.
Grundlage der FFT ist die Möglichkeit aus zwei schon berechneten Hälften einer DFT eine ganze zu erzeugen.

Bild 11-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 11-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:
DFT für N=1

Bild 11-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 11-2: Nicht rekursive Implementierung eines FFT-Algorithmus