Eine einfache FFT – Fast Fourier Transformation – mit dem Arduino Uno

Ein Rechtecksignal kann aus verschiedenen Sinussignalen erstellt werden indem man die einzelnen Sinusschwingungen addiert. Um ein Rechtecksignal zu erhalten nimmt man die Grundfrequenz f mit der Amplitude 1 und addiert die dreifache Frequenz 3*f mit der Amplitude 1/3 und 5*f mit 1/5 Amplitude usw. Das ergibt dann eine Rechteckschwingung.

 f(t) = sin(\omega t) + 1/3 sin(3 \omega t) + 1/5 sin(5 \omega t) + 1/7 sin(7 \omega t) ...

Diese Fourierreihe beschreibt ganz allgemein die Frequenzbestandteile einer Schwingung. Siehe Wikipedia.

Hat man bereits eine Schwingung vorliegen und möchte sie zerlegen in die einzelnen Frequenzen kann man ein Fourierzerlegung machen, auch Fourier Transformation genannt. Für Mikrokontroller und andere Programme wurde eine schnelle Frequenzzerlegung geschrieben, die FFT oder Fast Fourier Transformation um diese Zerlegung fast im Echtzeit machen zu können.

Eine FFT kann A/D Wandler Messwerte aus dem Zeitbereich (Wellenform) in den Frequenzbereich (Frequenzspektrum) übertragen. Dazu werden die Messwerte der A/D Wandlung aufgenommen, mit der fix_fft(..) Funktion berechnet und als Spektrum ausgegeben.

Eine gute Anleitung dazu ist auf dieser Seite. Ich habe ein paar Veränderungen vorgenommen, sodass das Programm auf dem Arduino Uno und der Arduino Umgebung 1.6.1 laufen kann. Die 8-Bit char Werte habe ich durch 16-Bit int Werte ersetzt, damit mein Arduino später mit 10 Bit A/D Wandlung arbeiten kann.

Als ersten Schritt berechne ich die Welle mit einer vorgegebenen Frequenz 40Hz, schicke diese durch die fix_fft() Funktion und gebe sie über den Serial Monitor aus mit der Funktion print_fft(). Dann warte ich eine Sekunde, addiere ich 40 Hz und wiederhole das ganze.

Das Programm

// FFT Berechnungen von dieser Seite
// http://members.shaw.ca/el.supremo/Arduino/fft/ffttest.htm
// Änderungen:
// 16 Bit integer Werte, dann kann später der 10-Bit Wandler des Arduino Uno verwendet werden
//
// Matthias Busse 11.4.2015 Version 0.2

#define N_WAVE	256    // volle Länge der Sinewave[] möglich 16, 64, 256, 1024 ...
#define LOG2_N_WAVE 8  // log2(N_WAVE) 
#define F_SAMPLE 5120  // 5120/128 = 40Hz per FFT "bin"
#define AMP_CARRIER 68 // amplitude of each "carrier" wave

int im[N_WAVE/2];
int re[N_WAVE/2];
int Sinewave[N_WAVE/2] = { // 128 Werte, die halbe positive Sinuswelle
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,
};

int i; 
double freq, phase;    // Frequenz, Phase

void setup() {
  Serial.begin(57600);
  freq = 40;
}

void loop() {
  phase = 0.0;
  for(i=0;i < 128;i++) {
    im[i] = 0; // Imaginärteil ist immer 0
    re[i] = (int)(AMP_CARRIER*sin(phase*2.0*PI)); // Welle 1 berechnen
    phase += freq/F_SAMPLE; // Phase von Welle 1 updaten
    if(phase >= 1)phase -= 1;
  }
  fix_fft(re,im,7,0); // die FFT der Messwerte berechnen
  print_fft((char *)"Eine Frequenz aufgenommen mit 5120Hz");
  delay(1000);
  freq += 40;
}

int fix_fft(int fr[], int fi[], int m, int inverse){
// fix_fft() - Forward / Inverse Fast Fourier Transform.
// fr[n],fi[n] Real und Imaginär Teil Felder für Eingabe und Ergebnis
// mit 0 <= n < 2**m; 
// inverse=0 heisst forward transform (FFT), oder =1 für iFFT.
int mr, nn, i, j, l, k, istep, n, scale, shift;
int qr, qi, tr, ti, wr, wi;
int idx;
  n = 1 << m;
  if (n > N_WAVE) return -1;  /* max FFT size = N_WAVE */
  mr = 0;
  nn = n - 1;
  scale = 0;
  for (m=1; m <= nn; ++m) { /* decimation in time - re-order data */
    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) {
      shift = 0; /* variable scaling, depending upon data */
      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 { // nicht invers
     shift = 1;
    }
    istep = l << 1;
    for (m=0; m < l; ++m) {
      j = m << k;
      if((idx = j+N_WAVE/4) >= 128) wr = -Sinewave[idx - 128];
      else wr = Sinewave[idx];
      if(j >= 128) wi = Sinewave[j];
      else 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;
}

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

void print_fft(char *title) {
  int i,j,largest;
  char str[65];
  char linfo[6];

  str[64] = 0;  
  largest = 0;
  for (i=0; i < 64;i++){ // Die größe Amplitude gint den maximalen Y-Achsenwert
    re[i] = sqrt(re[i] * re[i] + im[i] * im[i]);
    if(re[i] > largest) largest = re[i];
  }
  Serial.println("");
  Serial.println(title);
  Serial.println("");
  for(j=largest;j >= 0;j--) { // die Grafik ausgeben, mit der höchsten Amplitude anfangan
    for(i=0;i < 64;i++) {
      if(re[i] >= j)str[i] = '*'; // Wenn die Amplitude dieses Wertes mindesten so gross ist wie die aktuelle, einen * ausgeben
      else str[i] = ' '; // sonst ein Leerzeichen ausgeben
    }
    sprintf(linfo,"%3d ",j);
    Serial.print(linfo);
    Serial.println(str);
  }
  Serial.print("    ");
  for(i=0;i < 64;i++) { // Die X-Achsen Skala unten
    Serial.print(i%10);
  }
  Serial.println("");
  Serial.print("    ");
  for(i=0;i < 64;i++) {
    if(i < 10)Serial.print(" ");
    else Serial.print((i/10)%10);
  }
  Serial.println("");
}

Und hier die Ausgabe der Sinusschwingung auf dem Seriellen Monitor

fft-serial

Die X-Achse hat 40 Hz Schritte, sodass 8 = 8*40 = 320 Hz darstellt, und das ist hier abgebildet. Die anderen Frequenzen haben eine deutlich kleinere Amplitude und bewegen sich häufig bei 0-1.

Das oben beschriebene Rechtecksignal kann folgendermassen mit 80 Hz gebildet werden.

  for(i=0;i < 32;i++) {
    im[i] = 0; // Imaginärteil ist immer 0
    re[i] = (int)(AMP_CARRIER); // Welle 1 berechnen
    phase += freq/F_SAMPLE; // Phase von Welle 1 updaten
    if(phase >= 1)phase -= 1;
  }
  for(i=32;i < 64;i++) {
    im[i] = 0; // Imaginärteil ist immer 0
    re[i] = (int)(-AMP_CARRIER); // Welle 1 berechnen
    phase += freq/F_SAMPLE; // Phase von Welle 1 updaten
    if(phase >= 1)phase -= 1;
  }
  for(i=64;i < 96;i++) {
    im[i] = 0; // Imaginärteil ist immer 0
    re[i] = (int)(AMP_CARRIER); // Welle 1 berechnen
    phase += freq/F_SAMPLE; // Phase von Welle 1 updaten
    if(phase >= 1)phase -= 1;
  }
  for(i=96;i < 128;i++) {
    im[i] = 0; // Imaginärteil ist immer 0
    re[i] = (int)(-AMP_CARRIER); // Welle 1 berechnen
    phase += freq/F_SAMPLE; // Phase von Welle 1 updaten
    if(phase >= 1)phase -= 1;
  }

Und hier das zugehörige Spektrum

fft-rechteck

Es wurden verwendet:
Hardware: Arduino Uno
Software: Arduino 1.6.1

von Matthias Busse

Ein Gedanke zu „Eine einfache FFT – Fast Fourier Transformation – mit dem Arduino Uno

  1. Pingback: Arduino FFT auf dem 5110 Display ausgeben | Shelvin – Elektronik ausprobiert und erläutert

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert

Diese Website verwendet Akismet, um Spam zu reduzieren. Erfahre mehr darüber, wie deine Kommentardaten verarbeitet werden.