Home / Tutorials / DSP / DFT per Hand implementieren

DFT per Hand implementieren

Die Fouriertransformation (FT) ist ein grundlegendes und sehr nützliches Werkzeug. Sie bildet die Grundlage für die digitale Signalverarbeitung. Vereinfacht kann man sagen, dass man mit ihrere Hilfe ein Signal aus dem Zeitbereich in den Frequenzbereich transformieren kann.

Ein Zeitbereichssignal sagt uns wie sich etwas im Laufe der Zeit verändert. Das kann eine Aufzeichnung der Temperatur über einen bestimmten Zeitraum sein oder aber auch ein Tonsignal, wie z.B. Musik. Bei diesen Signalen kann man direkt erkennen zu welchem Zeitpunkt etwas passiert, wie sich das Signal mit der Zeit verändert.
Transformiert man dieses Signal nun in den Frequenzbereich, dann erkennt man welche Frequenzen in dem Signal „stecken“. Sowas wird auch Spektralanalyse genannt. Bei Musik kennt man das vielleicht vom Visualizer von Stereoanlagen.

Die Standardformel für die DFT lautet:

{X}_k = \sum_{n=0}^{N-1} x_n e^{\frac{-j2 \pi nk}{N}} 

hierbei gilt:

  • x ist die zu transformierende Zahlenfolge im Zeitbereich
  • X ist die Zahlenfolge im Frequenzbereich
  • N ist die Anzahl an Elemente der Beiden Zahlenfolgen
  • j ist die imaginäre Einheit

Man muss sich den obigen Ausdruck in etwa so vorstellen: Meist hat man ein Signal (Sensorwerte, Tonaufnahme oder ähnliches) von dem man gerne die FT bilden möchte, sagen wir um die Frequenzanteile zu erhalten. Das besagte Signal muss natürlich diskret vorliegen, als Zahlenfolge. Das ist unser x.
Das Ergebnis dieser Transformation ergibt eine komplexe Zahlenfolge {X}, das heißt jedes Element von {X} ist eine komplexe Zahl. Der Betrag dieser Zahl spiegelt die Signalstärke eben dieser Frequenz wieder, die das jeweilige Element repräsentiert. Es steckt auch eine Phaseninformation in dieser komplexen Zahl, die wir erst einmal ignorieren. Das wird erst später interessant.

Zur Veranschaulichung schreiben wir den Ausdruck etwas ausführlicher hin, für eine Samplelänge von 4:

{X}_k = x_0 + x_1 e^{-j \frac{2 \pi}{8}k} + x_2 e^{-j \frac{2 \pi}{8}2k} + x_3 e^{-j \frac{2 \pi}{8}3k}

nehmen wir an x ist hat 4 Samples und ist ein Einheitssprung, dessen Sprung genau in der Mitte befindet. Dann sieht unser Signal so aus

x=[0,0,1,1]

Dann sehen alle Terme so aus:

X[0] = 0 + 0 \cdot 1 + 1 \cdot 1 + 1 \cdot 1
X[1] = 0 + 0 \cdot e^{-j \frac{2 \pi}{8}} + 1 \cdot e^{-j \frac{2 \pi}{8}2} + 1 \cdot e^{-j \frac{2 \pi}{8}3}
X[2] = 0 + 0 \cdot e^{-j \frac{2 \pi}{8}2} + 1 \cdot e^{-j \frac{2 \pi}{8}4} + 1 \cdot e^{-j \frac{2 \pi}{8}6}
X[3] = 0 + 0 \cdot e^{-j \frac{2 \pi}{8}3} + 1 \cdot e^{-j \frac{2 \pi}{8}6} + 1 \cdot e^{-j \frac{2 \pi}{8}9}

Gedanken zur Implementierung

Der Ausdruck e^{-j \frac{2 pi}{N}nk} sieht noch etwas „fies“ aus. Je nachdem für welches System man die Implementierung macht, kann man ihn sich direkt berechnen lassen. Aber gerade bei Mikrocontrollersystemen braucht man dafür meistens eine zusätzliche Library, deren Einbindung man sich wohl überlegen sollte. Wenn man jedoch genau hinsieht ist dieser Term nur von n,k und N abhängig. N legen wir selber fest und für alle Produkte von n cdot k lässt sich der Wert recht schmerzlos vor raus berechnen, da er nicht von unserem Signal abhängt.

Denn Ausdruck können wir noch etwas umschreiben: e^{\frac{-j2 \pi }{N}}=\omega

latex$ omega$ ist die sog. Einheitswurzel. Diese Ist für ein Gegebenes N eine feste komplexe Zahl und man kann jetzt die obige Formel etwas verkürzt schreiben:

{X}_k = \sum_{n=0}^{N-1} x_n \omega^{nk}

Die Summe sieht ausgeschrieben für ein gegebenes k so aus:
{X}_k =x_{0} \cdot \omega^{0}+ x_{1} \cdot \omega^{k} + x_{2} \cdot \omega^{2k} +...+ x_{N-1} \cdot \omega^{(N-1)k}

Die Zahlenfolge x kann man als Vektor auffassen genauso wie die Werte für die \omega Potenzen. Die Summe selber lässt sich dann einfach, für ein gegebenes k, als Vektor-Vektor Multiplikation schreiben:

{X}_k=\vec \omega_{k} \cdot \vec x

Wenn man X nun wieder als Vektor ansieht, kann man alle Vektor-Vektor Multiplikationen als Matrix-Vektor Multiplikation interpretieren.

W = \begin{bmatrix}  1 & 1 & 1 & 1 & \cdots &1 \\  1&\omega&\omega^2&\omega^3&\cdots&\omega^{N-1} \\  1&\omega^2&\omega^4&\omega^6&\cdots&\omega^{2(N-1)}\\ 1&\omega^3&\omega^6&\omega^9&\cdots&\omega^{3(N-1)}\\  \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\  1&\omega^{N-1}&\omega^{2(N-1)}&\omega^{3(N-1)}&\cdots&\omega^{(N-1)(N-1)}\\  \end{bmatrix}

\large{X} = \large{W} \cdot \vec{x}

Die Matrix W wird in der englischen Literatur meist sinnigerweise als DFT-Matrix bezeichnet. Prozessoren die Vektor-Matrix Operation direkt durchführen können (z.B. GPUs) profitieren extrem von dieser Darstellungsweise.

Man erkennt hier besonders gut, wie der Aufwand für die Berechnung einer DFT skaliert. Es bedarf hierbei N komplexe Multiplikationen pro Zeile. Da wir N Zeilen haben, liegt die Anzahl der Multiplikationen bei N^2. Zusätzlich kommen noch (N-1)N Additionen vor. Der Aufwand wächst also recht schnell mit höheren Frequenzen.

In der unteren Abbildung sieht man einen Plot der ersten 6 Zeilen der DFT-Matrix. In rot sieht man den Realanteil, in blau den Imaginärteil.

bb232ed5e0107ba4a7f065c8ef3e0be4

Beispiele
Ich verwende zuerst einmal Numpy, ein Modul für Python, und erstelle mit dem FFT-Befehl eine Standartreferenz mit der wir unsere anderen Implementierungen vergleichen können. Als Testsignal nehmen wir ein Sinussignal mit einer Frequenz von 2Hz und tasten es 100 mal pro Sekunde ab. Diese Signal nennen wir x.

e1239cf55ba8f51621da207c9325272d

import numpy as np
import matplotlib.pyplot as plt

f_0 = 2
Fs = 100

t = np.linspace(0.1, 2*np.pi*f_0, Fs)
x = np.sin(t)
X_large = np.fft.fft(x)

plt.stem(abs(X_large))
plt.title('FFT', fontweight='bold', fontsize=20)
plt.ylabel('Amplitude', fontweight='bold', fontsize=18)
plt.xlabel('Frequency', fontweight='bold', fontsize=16)
plt.show()

Wendet man nun die den FFT Standartbefehl aus Numpy auf diese Zahlenfolge an und plottet den Absolutwert des Ergebnisses erhält man X. Die Achsen sind nicht auf physikalische Einheiten skaliert, weil man sonst nicht sehen kann welche Ergebnisse berechnet wurden.
Exemplarisch sind hier die ersten 6 Absolutwerte von X:

Element: 1 2 3 4 5 6
Wert: 0.010 0.137 49.947 0.248 0.137 0.098

85397e1cac1df73026070b2f1a9adba8

Die DFT lässt sich in Python z.B. so recht einfach implementieren:


def DFT(x):
  N=x.size;
  n = np.arange(N)
  W = np.exp(- 2j * np.pi * n.reshape(100,1) * n / N)
  return np.dot(W, x)

Leave a Reply

Your email address will not be published. Required fields are marked *

*

x

Check Also

An Arduino Unu with the word Assembler superimposed

Arduino Assembler: Einführung

Dieser Beitrag ist Teil 1 von 3 der Serie: Arduino AssemblerArduino AssemblerArduino Assembler: EinführungArduino Assembler: ...