8 Bit FFT übernehmen

[Josef (Gast)]

Hallo zusammen,

ich versuche gerade eine FFT auf einem 8 Bit Controller zu 
implementieren. Da ich das Rad ja nicht neu erfinden muss, versuche ich 
bereits vorhanden Code zu übernehmen: 
http://forum.arduino.cc/index.php/topic,38153.0.html

Als ersten Schritt probiere ich das ganze erstmal auf PC zum laufen zu 
bekommen. Am Code habe ich folgendes angepasst:
-Das Array "normal" als char angelegt statt im Flash
-"Normales" Auslesen des Arrays statt mit pgm_read_word_near()
-Funktionsaufruf hinzugefügt
-Ausgabe der Ergebnisse per printf

Inzwischen lässt sich das Projekt kompilieren aber das Ergebnis scheint 
nicht zu stimmen. Wenn ich es richtig verstehe, müsste das Ergebnis doch 
ins Array geschrieben werden, dann sollte doch bei einem Sinus ein 
bestimmter Wert im Array stark ausschlagen? Stattdessen sind aber alle 
denkbaren Werte im Array vertreten (inclusive negativer Werte). Falls 
sich jemand mit der FFT auskennt und mir einen Tipp geben könnte, wäre 
ich wirklich sehr dankbar, ich komme hier einfach nicht weiter.

Die Kommentare habe ich herausgenommen, weil der Code hier schon so lang 
war. Im Original (siehe Link) sind die ausführlichen Kommentare ja noch 
zu lesen.

#include <stdio.h>
#include <stdlib.h>


int fix_fft(char fr[], char fi[], int m, int inverse);
int fix_fftr(char f[], int m, int inverse);

int main(int argc, char** argv) 
{
  #define N_WAVE      256    /* full length of Sinewave[] */
  #define LOG2_N_WAVE 8      /* log2(N_WAVE) */

  int x;

char Sinewave[N_WAVE-N_WAVE/4] = {
0, 3, 6, 9, 12, 15, 18, 21, 
24, 28, 31, 34, 37, 40, 43, 46, 
48, 51, 54, 57, 60, 63, 65, 68, 
71, 73, 76, 78, 81, 83, 85, 88, 
90, 92, 94, 96, 98, 100, 102, 104, 
106, 108, 109, 111, 112, 114, 115, 117, 
118, 119, 120, 121, 122, 123, 124, 124, 
125, 126, 126, 127, 127, 127, 127, 127, 

127, 127, 127, 127, 127, 127, 126, 126, 
125, 124, 124, 123, 122, 121, 120, 119, 
118, 117, 115, 114, 112, 111, 109, 108, 
106, 104, 102, 100, 98, 96, 94, 92, 
90, 88, 85, 83, 81, 78, 76, 73, 
71, 68, 65, 63, 60, 57, 54, 51, 
48, 46, 43, 40, 37, 34, 31, 28, 
24, 21, 18, 15, 12, 9, 6, 3, 

0, -3, -6, -9, -12, -15, -18, -21, 
-24, -28, -31, -34, -37, -40, -43, -46, 
-48, -51, -54, -57, -60, -63, -65, -68, 
-71, -73, -76, -78, -81, -83, -85, -88, 
-90, -92, -94, -96, -98, -100, -102, -104, 
-106, -108, -109, -111, -112, -114, -115, -117, 
-118, -119, -120, -121, -122, -123, -124, -124, 
-125, -126, -126, -127, -127, -127, -127, -127, 

/*-127, -127, -127, -127, -127, -127, -126, -126, 
-125, -124, -124, -123, -122, -121, -120, -119, 
-118, -117, -115, -114, -112, -111, -109, -108, 
-106, -104, -102, -100, -98, -96, -94, -92, 
-90, -88, -85, -83, -81, -78, -76, -73, 
-71, -68, -65, -63, -60, -57, -54, -51, 
-48, -46, -43, -40, -37, -34, -31, -28, 
-24, -21, -18, -15, -12, -9, -6, -3, */
};


inline char FIX_MPY(char a, char b)
{
    /* shift right one less bit (i.e. 15-1) */
    int c = ((int)a * (int)b) >> 6;
    /* last bit shifted out = rounding-bit */
    b = c & 0x01;
    /* last shift + rounding bit */
    a = (c >> 1) + b;
    return a;
}


int fix_fft(char fr[], char fi[], int m, int inverse)
{
    int mr, nn, i, j, l, k, istep, n, scale, shift;
    char qr, qi, tr, ti, wr, wi;

    n = 1 << m;

    /* max FFT size = N_WAVE */
    if (n > N_WAVE)
        return -1;

    mr = 0;
    nn = n - 1;
    scale = 0;

    /* decimation in time - re-order data */
   for (m=1; m<=nn; ++m) {
        l = n;
        do {
            l >>= 1;
        } while (mr+l > nn);
        mr = (mr & (l-1)) + l;

        if (mr <= m)
            continue;
        tr = fr[m];
        fr[m] = fr[mr];
        fr[mr] = tr;
        ti = fi[m];
        fi[m] = fi[mr];
        fi[mr] = ti;
    }

    l = 1;
    k = LOG2_N_WAVE-1;
    while (l < n) 
  {
        if (inverse) 
    {
          //          variable scaling, depending upon data 
            shift = 0;
            for (i=0; i<n; ++i) 
      {
                j = fr[i];
                if (j < 0)
                    j = -j;
                m = fi[i];
                if (m < 0)
                    m = -m;
                if (j > 16383 || m > 16383) 
        {
                    shift = 1;
                    break;
                }
            }
            if (shift)
                ++scale;
        } 
    else 
    {
            shift = 1;
        }    
        istep = l << 1;
        for (m=0; m<l; ++m) 
    {
        j = m << k;
   //   wr =  pgm_read_word_near(Sinewave + j+N_WAVE/4);
        wr =  Sinewave [j+N_WAVE/4];
        wi = Sinewave[j];
        if (inverse)
            wi = -wi;
        if (shift) 
      {
                wr >>= 1;
                wi >>= 1;
            }
            for (i=m; i<n; i+=istep) 
      {
                j = i + l;
                tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
                ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
                qr = fr[i];
                qi = fi[i];
                if (shift) 
        {
                    qr >>= 1;
                    qi >>= 1;
                }
                fr[j] = qr - tr;
                fi[j] = qi - ti;
                fr[i] = qr + tr;
                fi[i] = qi + ti;
            }
        }
        --k;
        l = istep;
    }
    return scale;
}

int fix_fftr(char f[], int m, int inverse)
{
    int i, N = 1<<(m-1), scale = 0;
    char tt, *fr=f, *fi=&f[N];

    if (inverse)
        scale = fix_fft(fi, fr, m-1, inverse);
    for (i=1; i<N; i+=2) {
        tt = f[N+i-1];
        f[N+i-1] = f[i];
        f[i] = tt;
    }
    if (! inverse)
        scale = fix_fft(fi, fr, m-1, inverse);
    return scale;
}

  x = fix_fftr(&Sinewave[0], 6, 0);
  
    for(x=0;x<(N_WAVE-N_WAVE/4);x++)
    {
      printf(" %d",Sinewave[x]);  
    }
}

[Georg G. (df2au)]

Beim Lesen des Artikels hast du einen der letzten Sätze übersehen:

Unfortunately, there is a flow in the code, no FFT performed
if you copy/paste as it is.

[Thosch (Gast)]

Falls du es auf einem Atmel AVR implementieren möchtest, schau dir mal 
die Routinen von Chan an: http://elm-chan.org/works/akilcd/report_e.html

Der hat 'ne FFT-Routine in AVR-Assembler geschrieben, aber mit Interface 
zur Verwendung aus C-Code heraus.

Die Library kannst du auf der verlinkten Seite unten runterladen, die 
Sourcen für den Audio-Analyser (als Beispiel für die Verwendung der Lib) 
ebenso.

[Josef (Gast)]

Georg G. schrieb:
> Beim Lesen des Artikels hast du einen der letzten Sätze übersehen:
>
> Unfortunately, there is a flow in the code, no FFT performed
> if you copy/paste as it is.

Verdammt, das habe ich wirklich übersehen. Danke für den Hinweis. Aber 
es stellt sich die Frage, was man tun muss, um den Code zu übernehmen. 
Der User focalist hat es ja scheinbar hinbekommen.

Wie auch immer, solange ich damit nicht weiterkomme, sehe ich mir mal 
die Routinen von Chan an. Fragen dazu folgen ;)

[Josef (Gast)]

Achso, der Controller steht noch nicht ganz fest, aber es soll 
wahrscheinlich irgendein PIC werden. Aber das müsste sich ja recht leid 
adaptieren lassen, oder?

[Karl H. (kbuchegg)]

Vor allen Dingen mal aufräumen im Code.
Wenn ich for Schleifen sehe, die bei 1 anfangen, fangen bei mir die 
Alarmglocken zu schrillen an.
Der Datentyp char ist schlecht gewählt. Kein Mensch weiß´, ob 'char' auf 
diesem Compiler ein signed oder ein unsigned Datentyp ist. Immer 
explizit sein! Es gibt 3(!) kleine Datentypen

* signed char
  Für kleine Integer mit Vorzeichen

* unsigned char
  Für kleine Integer ohne Vorzeichen. Das ist der Datentyp der Wahl, 
wenn man es mit Bytes zu tun hat

* char
  Für alles was mit Textverarbeitung zu tun hat.

Hier hat man es nicht mit Textverarbeitung zu tun. Daher ist 
plain/vanilla char nicht angebracht, sondern ein 'ask for trouble'.

[Josef (Gast)]

Josef schrieb:
> Achso, der Controller steht noch nicht ganz fest, aber es soll
> wahrscheinlich irgendein PIC werden. Aber das müsste sich ja recht leid
> adaptieren lassen, oder?
Moment, da habe ich nicht nachgedacht. Die Libarys gehen ja nur für AVR. 
Dann hilft mir das wohl nicht so einfach weiter, oder?

Karl Heinz Buchegger schrieb:
> Immer explizit sein!
Wird gemacht :)

Karl Heinz Buchegger schrieb:
> ist ziemlich offensichtlicher Blödsinn. Das sollte offenbar ein
> Swap-Code sein, so wie das hier steht, ist es eine 0-Operation.
Ich verstehe nicht ganz, was du meinst, fr[m] und fr[mr] werden doch 
getauscht, oder nicht? Wie müsste es denn anders heißen?

[Karl H. (kbuchegg)]

Josef schrieb:

> Karl Heinz Buchegger schrieb:
>> ist ziemlich offensichtlicher Blödsinn. Das sollte offenbar ein
>> Swap-Code sein, so wie das hier steht, ist es eine 0-Operation.
> Ich verstehe nicht ganz, was du meinst, fr[m] und fr[mr] werden doch
> getauscht, oder nicht? Wie müsste es denn anders heißen?


Ich hab diesen Einwand wieder zurück gezogen, weil ich nach dem 5-ten 
mal Code lesen gemerkt habe, dass das eine mal die Indexvariable m 
heißt, während sie das andere mal mr heißt. So einfach sieht man das 
nicht und das fällt für mich genauso unter Code-aufräumen: 
Variablennamen so wählen, dass man aus 2 Meter Entfernung mit 
Sonnenbrille erkennen kann, dass es sich um unterschiedliche Variablen 
handelt. Persönlich mache ich seit langem einen großen Bogen um 
veröffentlichten Code, bei dem 90% aller Variablen einen ein- oder 
zweibuchstabigen Namen tragen.


Aus dem hier werde ich nicht schlau

            shift = 0;
            for (i=0; i<n; ++i) 
      {
                j = fr[i];
                if (j < 0)
                    j = -j;
                m = fi[i];
                if (m < 0)
                    m = -m;
                if (j > 16383 || m > 16383)

mal angenommen, dass ein char tatsächlich ein signed char ist. Dann 
ändert sich der Zahlenwert bei der Zuweisung an den int j nicht. Auch 
nicht bei negativen Zahlen. Aber wie soll ja dann größer als 16383 
werden?
Angenommen char ist in Wirklichkeit in unsigned char. Dann werden 
'negtive Werte' aus dem Original im Bereich 128 bis 255 abgebildet. 
Wieder gibt es keine Chance jemals größer als 16383 zu werden.

So richtig hab ich noch nicht durchschaut, was hier eigentlich passieren 
soll.

[Josef (Gast)]

OK, bei dem Code kann man sich wirklich schnell man vertun. Aufräumen 
würde ich trotzdem nur ungern, weil mir das Risiko hoch erscheint, dass 
ich dabei nur weitere Fehler einprogrammiere. Das ist jetzt schon das 
dritte Code-Beispiel, dass ich erfolglos versuche zum Laufen zu bringen 
- Wirklich frustrierend! Falls jemand noch etwas für PIC Controller 
kennt, wäre das auch super.

Bei dem Codestück muss ich dir Recht geben, die if-Abfrage müsste 
tatsächlich nie erfüllt werden ... Ratlosigkeit

[Karl H. (kbuchegg)]

Ich müsste mich jetzt auch wieder in die FFT und den Butterfly 
Algorithmus einlesen, wie der konkret funktioniert.
Wenn du das ganze sowieso am PC laufen hast, hindert dich ja nichts 
daran, den Code mit printf zu spicken, erst mal mit einer kleineren 
Datenmenge zu operieren und dann anhand der Ausgaben kontrollieren, ob 
die richtigen Werte (sowohl Index als auch Zahlenwerte) miteinander 
verrechnet werden.

[Georg G. (df2au)]

Karl Heinz Buchegger schrieb:
> Angenommen char ist in Wirklichkeit in unsigned char. Dann werden
> 'negtive Werte' aus dem Original im Bereich 128 bis 255 abgebildet.
> Wieder gibt es keine Chance jemals größer als 16383 zu werden.

Im Kopf der fft.cpp schreibt der Autor, dass alle Werte short int sind 
und von -32768 bis +32767 laufen. Float war ihm zu langsam und zu 
unbequem. Vermutlich ist einer der Fehler im Code, dass fr[] und fi[] 
als char definiert sind. In seinem Beispiel Code weist er auch int Werte 
zu und nicht char.

[Josef (Gast)]

Karl Heinz Buchegger schrieb:
> Ich müsste mich jetzt auch wieder in die FFT und den Butterfly
> Algorithmus einlesen, wie der konkret funktioniert.
Geht mir auch so, nur dass ich nicht sicher bin, ob ich überhaupt in der 
Lage bin den Stoff wirklich zu verstehen.

Karl Heinz Buchegger schrieb:
> Wenn du das ganze sowieso am PC laufen hast, hindert dich ja nichts
> daran, den Code mit printf zu spicken, erst mal mit einer kleineren
> Datenmenge zu operieren und dann anhand der Ausgaben kontrollieren, ob
> die richtigen Werte (sowohl Index als auch Zahlenwerte) miteinander
> verrechnet werden.
Dafür mache ich es erstmal am PC. Bisher ist der FFT-Code aber eine 
Black Box für mich. So kann ich natürlich auch nichts debuggen.

Georg G. schrieb:
> Im Kopf der fft.cpp schreibt der Autor, dass alle Werte short int sind
> und von -32768 bis +32767 laufen. Float war ihm zu langsam und zu
> unbequem. Vermutlich ist einer der Fehler im Code, dass fr[] und fi[]
> als char definiert sind. In seinem Beispiel Code weist er auch int Werte
> zu und nicht char.
Stimmt. Gerade habe ich einfach mal blauäugig alles zu unsigned short 
int gemacht, aber natürlich ohne Erfolg. Ich frage mich nur, wie es sein 
kann, dass der Code beim Autor und auch anderen funktioniert hat. Dann 
kann doch eigentlich nichts so grundsätzliches falsch sein.

[Georg G. (df2au)]

Als Zwischenergebnis:
In dem Paket sind zwei FFT Funktionen, für real/imaginär Daten und für 
nur-real-Daten.

Die nur-real-Daten Funktion hat in der Theorie den Charme, das sie 
drastisch schneller ist. Leider funktioniert sie in der angegebenen Form 
nicht.

Die real/imaginär Daten Funktion funktioniert gut, wenn man nur real 
Daten vorgibt und das imaginäre Feld mit Nullen füllt. Mit imaginär 
Daten ungleich Null habe ich es nicht versucht.

Es sind einige verbesserungsfähige Stellen in dem Programm. Als 
Beispiel: Man muss kein int16_t aus einer Tabelle von int8_t Werten 
lesen. Die Kommentare sind teilweise - wohl aus historischen Gründen - 
irreführend.

Vielleicht findet sich ja noch jemand, der die Macken sucht und hier 
publiziert.

[Frank K. (fchk)]

Josef schrieb:
> Achso, der Controller steht noch nicht ganz fest, aber es soll
> wahrscheinlich irgendein PIC werden. Aber das müsste sich ja recht leid
> adaptieren lassen, oder?

Wenn es ein dsPIC wird (und das wäre bei einer FFT vielleicht sinnig), 
hat sich Microchip schon darum gekümmert:

http://www.microchip.com/stellent/idcplg?IdcService=SS_GET_PAGE&nodeId=2680&dDocName=en023598

Die Bibliothek ist beim C30 Compiler dabei. Du darfst sie einfach 
benutzen.

fchk