FFT Algorithmus in C für ARM7

[Steffen Bauert (Gast)]

Hallo,

Hintergrund: Ich habe mir mal ein kleines Bord gemacht mit nem 
AT91SAM7S256, Display, AD-Wandler, SD-Card .. bla bla.
Damit habe ich schon einige Sachen experimentiert.

Jetzt wollte ich mal was mit einer FFT machen. Also eine Frequenz an den 
AD-Wandler und dann auf dem Display so ein schönes Spektrum anzeigen 
lassen ;)

Hat jemand einen Tipp oder gar eine FFT Funktion in C die ich mir dann 
mal anschauen kann die auf einem ARM läuft? Ich weis was die FFT macht, 
aber wie blicke ich noch nicht richtig durch. Daher wäre vielleicht eine 
FFT in C (in Assembler hab ich viele gefunden, aber auch nicht für ARM) 
für die Implementierung gut.

Vielen Dank!

[Düsentrieb (Gast)]

http://www.dspguru.com/sw/opendsp/lib.htm
http://www.gidf.de/

[Andreas S. (andreas) (Admin)]

http://www.mikrocontroller.net/articles/Digitale_Signalverarbeitung#FFT
-> http://sourceforge.net/projects/kissfft/

[Steffen Bauert (Gast)]

Hallo,
danke erstmal.. aber so richtig blicke ich noch nicht durch (eigentlich 
garnicht ;)


1.
Überall wird das Eingabe-Array als komplexe Zahl dargestellt (Realteil + 
Imaginärteil). Jetzt kommen meine Daten vom AD-Wandler (16Bit, 2er 
Komplement also 32768 bis -32768. Den Wert vom ADC mache ich ja (so wie 
ich es verstehe) in den Realteil des Arrays (also die geraden Indizes 
0,2,4...). Was mache ich in den Imaginärteil? Kann ich den einfach auf 0 
setzten?

2.
Bei wikipedia ist pseudocode angegeben..
http://de.wikipedia.org/wiki/Schnelle_Fourier-Transformation#Implementierung_.28Pseudocode.29
So wie ich das sehe kann man doch den Faktor e hoch -2(Pi)ik/n als 
konstante machen. k ist dabei der aktuelle Index (?), n die Anzahl der 
Samples. Aber was ist i? Bzw. wie berechne ich das? Wenn ich oben im 
Eingangsarray i = 0 setzte immer wäre doch der gesamte faktor immer 0.. 
überhauptnichtversteh


3.
Wenn ich jetzt die meine Funktion fft(in, out) habe. Sehe ich das 
richtig, das ich die Amplitude dann mit wurzel(re² + im²) bekomme?
Also z.B. wurzel(out[0]² + out[1]²) ist die Amplitude, die rechne ich 
dann in eine Spannung um analog dazu wie ich den ADC Wert in eine 
Spannung umrechnen würde und habe die Amplitude in Volt.. ?

Ich würde das ganze gerne selber Implementieren und auch irgendwie 
verstehen ;)

Danke für eure Hilfe!!!

[Christian (Gast)]

Hallo,

also das i (in der elektrotechnik auch als j bekannt) gehöhrt zur 
Polardarstellung einer komplexen Zahl. (z.b.: 35*e^i90) 35 ist der 
Radius des Zeiers und das 90 ist der Winkel(kann grad oder radiand sein, 
bei Software radiand, im Matheheft Grad). Das i kann man nicht berechnen 
aber i^2 gibt -1 wenns dir weiter hilft.

[Martin Wantke (Gast)]

Das e hoch i*p ist die Eulersche Formel, auch Eulersche Identität 
genannt. Das heißt: schreibt man z = e^(i*p) und z ist eine komplexe 
Zahlen (was ja im Grunde zwei reelle Zahlen sind->reelle und imaginäre), 
dann ist der reelle Teil = cos(p) und der imaginäre = sin(p).
Du übergibst der Fourier-Transformation ein Array aus komplexen Zahlen, 
wobei im Realteil alle Abtastwerte vom deinem AD-Wandler gespeichert 
sind und alle imaginären Teile des Arrays auf 0 gesetzt werden. Wenn du 
das Spektrum suchst benötigst von den transformierten Werten nur die 
Hälfte + 1. Hast du z.B. der FFT 32 komplexe Zahlen übergeben, brauchst 
du zum Schluss die ersten 17, die anderen Werte sind Spiegelungen mit 
negativer imaginären Teil. Jetzt brauchst du nur noch die Auslenkung der 
ersten 17 komplexen Zahlen ermitteln indem du jeweils r = Wurzel((real + 
real) + (imag * imag)) rechnet und dann hast du dein Spektrum.

Aus dem Pseudocode von Wikipedia hatte ich folgende C++ Routinen 
geschrieben, wurden aber auf wiki wieder entfernt:

vector<complex<double> > DFT(vector<complex<double> >f)
{
    int n = f.size();
    vector<complex<double> > r(n);
    complex<double> temp;
    for(int k = 0; k != n; k++)
    {
        for(int j = 0; j != n; j++)
        {
            double x = (-2.0  M_PI  k * j) / n;
            temp = complex<double>(cos(x), sin(x));
            r[k] = r[k] + (f[j] * temp);
        }
    }
    return r;
}

vector<complex<double> > FFT(vector<complex<double> >f)
{
    int n = f.size();
    if((n & n - 1) != 0)
    {
        cout << "n ist keine Zweierpotenz." << endl;
        vector<complex<double> > r(0); return r;
    }
    else if(n == 1)
    {
        return f;
    }
    vector<complex<double> > gPar(n / 2);
    vector<complex<double> > uPar(n / 2);
    for(int k = 0; k != n / 2; k++)
    {
        gPar[k] = f[k * 2];
        uPar[k] = f[k * 2 + 1];
    }
    vector<complex<double> > g = FFT(gPar);
    vector<complex<double> > u = FFT(uPar);
    complex<double> temp;

    vector<complex<double> > r(n);
    for(k = 0; k != n / 2; k++)
    {
        double x = (-2.0  PI  k) / n;
        temp = complex<double>(cos(x), sin(x));
        r[k] = g[k] + (u[k] * temp);
        r[k + (n / 2)] = g[k]- (u[k] * temp);
    }
    return r;
}

Ich hab die dft in Excel nachgebildet mit 8 Abtastwerten, dadurch hab 
ich die Transformation besser verstanden.
Auf
http://manfred.informatik.hu-berlin.de/php/fftw.php
kannst du die Werte der fft überprüfen.