Digitale Signale visualisieren

Einleitung

In diesen Artikel werde ich mich damit beschäftigen, wie Geräusche (z.B. Musikstücke, Maschinen und ähnliches) digital analysiert und das Ergebnis visualisiert werden kann.

Der Artikel soll Grundlagen erläutern, die notwendig sind, um diese Geräusche so zu transformieren, dass die für uns interessanten Eigenschaften sichtbar werden. Diese Art der Transformation zählt zu dem Gebiet der Digitalen Signalverarbeitung (im Englischen: Digital Signal Processing – DSP) und wird noch näher erläutert.

Die so gewonnenen Eigenschaften werden dann zum Ende des Artikels mit Hilfe eines Beispielprogrammes grafisch als einfacher Frequenzmonitor dargestellt. Zum Beispiel stellen bekannte Mediaplayer wie Winamp (und andere) so einen Monitor dar, um dem Hörer ein visuelles Feedback des aktuell laufenden Songs zu geben.

Einfacher Equalizermonitor von Winamp

Einfacher Equalizermonitor von Winamp

Diesem Beispielprogramm kann eine WAV oder OGG/Vorbis-Datei übergeben werden und als Resultat wird eine Videodatei (AVI) erzeugt, welche synchron zu dem Inhalt der Eingabedatei ausgewählte Frequenzen darstellt.

Das Beispielprogramm aus dem Artikel in Aktion

Das Beispielprogramm aus dem Artikel in Aktion

Die digitale Signalverarbeitung

Eine kurze Einführung

Im Folgenden werde ich die für das Beispielprogramm notwendigen Grundlagen erläutern. Dabei werde ich dies immer aus Sicht des Anwenders und nicht des Mathematikers darstellen. Für eine detaillierte Einführung in das Thema digitale Signalverarbeitung empfehle ich das Buch „Understanding Digital Signal Processing“. Das Buch stellt auf der einen Seite den praktischen Teil der Anwendung dar, lässt aber nicht die mathematischen Grundlagen außer Acht.

Am Ende des Artikels ist eine Liste von Offline- und Onlinequellen zu dem Thema zu finden.

Abtastfrequenz

Um ein analoges Signal (digital) verarbeiten zu können, muss es von seiner analogen Form in die digitale, diskrete Form überführt werden. Dies geschieht durch Abtasten des analogen Signales in einem vordefinierten festen Zeitabstand und der Umwandlung in diskrete Werte eines vordefinierten Amplitudenbereiches. Weil so ein diskretes Signal im Gegensatz zu dem analogen Original nur ein endlichen Wertebereich hat, kann es nur eine Annäherung an das analoge Vorbild sein.

Dieser zeitliche Abstand zwischen zwei Abtastzeitpunkten des analogen Signals wird als Abtastintervall T bezeichnet. Er bestimmt die sogenannte Abtastfrequenz fA, welche zum Beispiel in Hz, kHz oder MHz angegeben wird. Beträgt die Abtastfrequenz des analogen Signals zum Beispiel 22050 Hz, so errechnet sich das Abtastintervall wie folgt:

equation_01_t_1_div_fa_22050

Als Resultat dieser Abtastung erhält man eine Folge von 22050 Werten pro Sekunde. Ein Signal wird als Sequenz einer Folge von Signalwerten im Intervall [0,n] wie folgt beschreiben:

equation_02_sequence_0_to_n_real

Ein digitales Signal bestehend aus acht Werten könnte demnach wie folgt aussehen:

x(n) = {0.224, 0.555, 0.75, 0.8, 0.75, 0.555, 0.224, -0.224}

Diskrete Fourier Transformation (DFT)

Um das digitale Signal analysieren zu können, benötigen wir die Hilfe der diskreten Fourier Transformation (DFT). Die DFT ist wohl mit eine der am meisten benutzen Methoden in der digitalen Signalverarbeitung. Mit ihrer Hilfe kann eine Signalfolge in ihren Frequenzbereich analysiert, aber auch verändert werden. Ich werde mich aber in diesem Artikel auf die Analyse im Frequenzbereich beschränken.

Die DFT wird benutzt, um die in einer diskreten Signalsequenz enthaltene Frequenzen zu ermitteln und abzubilden. In ihrer Expotentialform stellt sich die DFT Gleichung wie folgt dar:

equation_05_dft_expotential_form

Wobei:

  • x(n) = die diskrete Signalsequenz ist
  • X(m) = die diskrete Sequenz der DFT im Frequenzbereich ist
  • N = die Anzahl der Werte der Ein- und Ausgabe der DFT ist

Wie genau sieht das Ergebnis der DFT nun aber aus? Einfach gesagt, stellt jeder Wert von X(m) den Betrag und die Phase einer bestimmten Frequenz dar. Für welche Frequenz genau der Index m steht, wird durch die Abtastfrequenz des Eingangssignals und die Breite N der DFT bestimmt. Beträgt zum Beispiel die Abtastfrequenz 22050 Hz und die Breite der DFT N=256 Elemente, so errechnet sich die Frequenz des zehnten (n=9) Wertes aus X(m) wie folgt:

equation_06_dft_freq_bin_index9

Der zehnte Wert der Ausgabe der DFT stellt also den Betrag und die Phase der Frequenz von 775,2 Hz dar.

Der Betrag und der Phasenwinkel in X(m) werden als komplexe Zahl dargestellt:

equation_07_Xm_real_imag

Der Betrag einer komplexen Zahl wird wie folgt berechnet:

equation_08_dft_magnitude

Der Phasenwinkel einer komplexen Zahl wird wie folgt berechnet:

equation_09_dft_phase_angle

Eigentlich würde es schon reichen, die Beträge der jeweiligen Elemente der DFT-Ausgabe zu berechnen und für den jeweiligen Elementindex dessen Frequenz zu ermitteln. Die daraus resultierenden Werte sind das Frequenzspektrum des Eingangssignals.

Dennoch gibt es einige besondere Eigenschaften einer DFT, welche bei der Interpretation des Frequenzspektrums beachtet werden müssen.

Symmetrie

Eine DFT kann Eingangswerte aus dem reellen Zahlenbereich oder auch komplexe Zahlen entgegennehmen. In diesem Artikel benutzen wir diskrete Signalsequenzen aus dem reellen Zahlenbereich. Als Konsequenz daraus, bilden die Ausgangswerte der DFT eine Symmetrie in der Form, dass Werte ab dem Index m=1 bis m=(N/2) – 1 für m > (N/2) gespiegelt werden. Da diese gespiegelten Werte keine neuen Informationen über das Frequenzspektrum liefern, werden diese in den meisten Fällen weggelassen. Deshalb liefert FFTW beim Aufruf der DFT als Resultat nur die Werte der Indizes 0 bis (N/2) + 1 zurück.

Das folgende Beispiel zeigt die DFT eines Signals, welches aus einer 50 Hz Sinuswelle (Amplitude = 1,0) und einer 70 Hz Sinuswelle (Amplitude = 0,25) besteht und mit 200 Hz abgetastet wird:

Beispiel für Symmetrie des Resultats einer DFT

Beispiel für Symmetrie des Resultats einer DFT

Deutlich ist die horizontale Symmetrie der X-Achse über der Mitte (N / 2 + 1) des Bildes erkennbar. Die Frequenz von 50 Hz wurde auf 150 Hz und die Frequenz von 70 Hz auf 130 Hz gespiegelt.

Dieses Verhalten hat aber noch einen weiteren Effekt: Hat zum Beispiel eine DFT N=512 reelle Eingangswerte, so sind tatsächlich nur (N / 2) + 1 = (512 / 2) + 1 = 257 Werte für die Analyse verwertbar. Wurde das Eingangssignal beispielsweise mit 22050 Hz abgetastet, so liegt das Frequenzspektrum effektiv nicht von 0 bis 22050 Hz. Stattdessen halbiert sich das Frequenzspektrum auf nur noch 0 bis:

equation_11_effective_frequency_bandwidth

Abtasttheorem

Das Abtasttheorem, auch als Nyquist-Shannon-Abtasttheorem oder Nyquist-Theorem, bekannt beschreibt eine wichtige Eigenschaft in der Signalverarbeitung. Es besagt, dass die Abtastfrequenz eines analogen Signals mehr als doppelt so groß sein muss, als die höchste im Signal enthaltene Frequenz:

equation_10_nyquist_theorem

Nur wenn Abtastfrequenz mehr als doppelt so groß ist, ist sichergestellt, dass dieses Signal vollständig abgetastet und auch wieder rekonstruiert werden kann. Ist die Abtastfrequenz zu gering, führt dies dazu, dass bestimmte Frequenzen sich gegenseitig überdecken und damit nicht mehr erkennbar sind. Dieser Effekt wird auch Aliasing genannt.

Sollen zum Beispiel aus einem Signal die Frequenzen von 0 bis 8 kHz per DFT analysiert werden, so sollte das analoge Signal mit einer Frequenz von z.B.

8 kHz * 2.5 = 20 kHz

abgetastet werden.

Zum Vergleich wurde für die beiden folgenden Bilder dasselbe Eingangssignal erzeugt. Das Signal besteht aus drei Sinuswellen aus je 10, 20 und 50 Hz. Im ersten Beispiel erfolgte eine Abtastung mit 180 Hz (180 Hz > 2 * 50 Hz) und im zweiten Beispiel eine Abtastung mit nur 50 Hz (70 Hz < 2 * 50 Hz).

Frequenzspektrum dreier Sinuswellen (10 Hz, 20 Hz, 50 Hz, alle mit Amplitude = 1,0), Abtastfrequenz = 180 Hz

Frequenzspektrum dreier Sinuswellen (10 Hz, 20 Hz, 50 Hz, alle mit Amplitude = 1,0), Abtastfrequenz = 180 Hz

Frequenzspektrum dreier Sinuswellen (10 Hz, 20 Hz, 50 Hz, alle mit Amplitude = 1,0), Abtastfrequenz = 70 Hz

Frequenzspektrum dreier Sinuswellen (10 Hz, 20 Hz, 50 Hz, alle mit Amplitude = 1,0), Abtastfrequenz = 70 Hz

Im dem Beispiel mit der zu geringen Abtastfrequenz von 70 Hz, ist nur noch die Frequenz von 10 Hz erkennbar. Die anderen beiden Frequenzen (20 Hz, 50 Hz) sind verschwunden. Die gut erkennbare Veränderung des absoluten Betrages des 10 Hz Signals im Gegensatz zu dem korrekt abgetasteten Signal kommt dadurch zustande, weil sich die Abtastfrequenz und damit die Anzahl der Eingangswerte verändert haben.

Leck Effekt

Wie in dem Artikel schon erwähnt wurde, ist für jeden Index der Ausgabewerte der DFT festgelegt, welcher Frequenzanteil in ihm enthalten ist. Hat die DFT einen Breite von 512 Elementen und wurde das Eingangssignal mit 8000 Hz abgetastet, so enthält jeder Index n den Frequenzanteil eines Vielfachen von f_n = 8000 Hz / 512 = 15,625 Hz.

Liegt nun ein Frequenzanteil in dem Signal nicht genau auf einem Vielfachen von f_n (was sehr wahrscheinlich ist), so kommt ein unerwünschter Effekt zum Tragen. Denn in diesem Fall werden Anteile der tatsächlichen Frequenz an Frequenzen in deren Umgebung addiert. Daraus resultiert auch eine Abschwächung des Frequenzanteiles. Dieses Phänomen ist in bei der Anwendung der DFT ein großes Problem.

10 Hz (blau) und 10,4 Hz (grün) Signal mit 100 Hz abgetastet

10 Hz und 10,4 Hz Signal mit 100 Hz abgetastet

Als Beispiel habe ich ein 10 und 10,4 Hz Signal mit 100 Hz abgetastet. Die Frequenzen sind in dem Diagramm grün und blau dargestellt. Die DFT hat eine Breite von 100 Elementen, somit liegt jeder Index auf einem Vielfachen von 1 Hz. Sehr gut erkennbar ist, dass das grüne Signal (10,4 Hz) seine Frequenzanteile an dessen Umgebung verteilt und dadurch dessen Betrag verringert wird. Durch diese neu entstandenen (Pseudo)-Frequenzen, wird das tatsächliche Signal verfälscht.

Diesen Effekt versucht man in der digitalen Signalverarbeitung mit Anwendung von sogenannten Fensterfunktionen (Windows) abzuschwächen.

Fensterfunktionen

Um die unerwünschten Leckeffekte abzuschwächen, wird auf dem Eingangssignal vor Durchführung der DFT, jeweils eine Fensterfunktion angewandt. Dabei werden die Amplituden des Eingangssignal verändert, um die Leckeffekte zu minimieren.

Es gibt verschiedene Fensterfunktionen – die bekanntesten dürften aber die Rechteck-, die Dreieck- die Hamming- und die Hanning-Fensterfunktion sein. Für den Softwareequalizer verwende ich die letztere Fensterfunktion.

Als Beispiel soll eine 10 Hz Frequenz in einem Eingangssignal dienen. Die DFT soll 200 Elemente breit sein.

10 Hz Eingangssignal

10 Hz Eingangssignal

Dasselbe Eingangssignal nach Anwendung der Hanning-Fensterfunktion

Dasselbe Eingangssignal nach Anwendung der Hanning-Fensterfunktion

Das erste Diagramm zeigt das Original-Eingangssignal – im zweiten Diagramm wurde die Hanning-Fensterfunktion auf die 200 Elemente angewandt.

In dem nächsten Diagramm wurde die Hanning-Fensterfunktion auf das Eingangssignal aus dem Abschnitt über die Leckeffekte angewandt. Das Eingangssignal enthält Frequenzanteile von 10 Hz und 10,4 Hz. In dem Diagramm ist der grüne Verlauf die 10,4 Hz Frequenz, welche ohne Fensterfunktion Leckeffekte verursacht. Der blaue Verlauf wurde mit der Hanning-Fensterfunktion erzeugt und unterbindet das Abgeben von Frequenzanteilen an die Umgebung von 10,4 Hz. Wiederum ist eine weitere Abschwächung der 10,4 Hz Frequenz durch die Fensterfunktion zu erkennen. Durch Anwendung einiger Fensterfunktionen (z.B. Hanning, Hamming), werden die Frequenzamplituden abgeschwächt, was nach der DFT zu verminderten Frequenzbeträgen führt.

Hanning-Fensterfunktion auf 10,4 Hz Signal angewandt

Hanning-Fensterfunktion auf 10,4 Hz Signal angewandt (grün = Leck Effekt, rot = abgeschwächter Leck Effekt)

FFT – die schnelle Version der DFT

Die DFT in ihrer am Anfang des Artikels dargestellten Form wäre zu langsam um damit in Echtzeit das Frequenzspektrum eines Signals analysieren und auch modifizieren zu können. Aus diesem Grund wird in der Praxis eine modifizierte Version der DFT angewandt, welche als schnelle Fourier Transformation (FFT) bezeichnet wird. In der Praxis gibt es mehrere Versionen der FFT.
Genaue Beschreibungen der einzelnen Umsetzungen (z.B. Radix-4) gibt es in der einschlägigen Literatur und füllt dort teilweise mehrere Kapitel. Ich wollte hier nur einen kurzen Hinweis auf das Thema geben.

weiterführende Links
http://de.wikipedia.org/wiki/Diskrete_Fourier-Transformation
http://de.wikipedia.org/wiki/Abtasttheorem
http://de.wikipedia.org/wiki/Leck-Effekt
http://de.wikipedia.org/wiki/Fensterfunktion
http://de.wikipedia.org/wiki/Schnelle_Fourier-Transformation

Beispiel: Der Frequenzmonitor

Wie schon zum Anfang des Artikels erläutert, habe ich für die Umsetzung des Frequenzmonitors mehrere 3rd-Party-Bibliotheken verwendet, welche im Folgenden näher erklärt werden.

FFTW („Fastest Fourier Transform in the West“)

Das Herzstück des Monitors ist die DFT. Auch wenn die Umsetzung einer DFT als Programmcode nicht so schwierig erscheint, kostet es sehr viel Zeit, diesen Code fehlerfrei, speichersparend und vor allem schnell zu machen. Die Bibliothek FFTW stellt unter anderen genau diese DFT mit den angesprochenen Anforderungen zur Verfügung. Sie wurde am MIT von Matteo Frigo entwickelt und steht unter GNU General Public License (GPL).

Die Bibliothek bietet eine einfach zu verwendende C-API, welche bei Bedarf aber auch angepasst werden kann. FFTW unterteilt die Schnittstellen in Basic, Advanced und Guru Interfaces. Weiterhin bietet FFTW die Möglichkeit, die DFT während der Laufzeit auf maximale Performanz zu trimmen. Es ist aber auch möglich, stattdessen die Verarbeitung auf eine geringe Speicherverwendung zu optimieren. Dies geht dann aber mit einem Performanzverlust einher.

FFTW kann unter Ubuntu mit folgendem Befehl installiert werden:

$ sudo apt-get install libfftw3-dev

Die Headerdateien sind unter Ubuntu im Verzeichnis

/usr/include

und die Libraries im Verzeichnis

/usr/lib

zu finden.

Link: http://www.fftw.org/

libsndfile

Um Sounddateien einlesen zu können, verwende ich die libsndfile-Bibliothek. Die Bibliothek kann viele Soundformate wie WAV, AIFF, VOC, FLAC und OGG/Vorbis einlesen. Sie bietet zwar keine Funktionalität um MP3-Dateien einzulesen, dies ist für diesen Artikel aber auch nicht notwendig. Im Notfall können über Audiocity MP3-Dateien nach OGG/Vorbis konvertiert werden.

Das OGG/Vorbis-Format sollte ohnehin MP3 vorgezogen werden, weil es eine bessere Soundqualität bietet und geringe Dateigröße hat. Für Signale von kurzer Dauer bieten sich WAV-Dateien, für längere Signale eher OGG/Vorbis-Dateien an.

Ähnlich wie FFTW bietet libsndfile eine einfach C-API mit verschiedenen Abstufungen im Detailgrad an. Auf der untersten Ebene erlaubt libsndfile das Lesen einzelner Bytes. Auf einer abstrakteren Ebene ist das Lesen von Frames möglich. Die Bibliothek passt die eingelesenen Daten auch automatisch auf den gewünschten Datentyp (short, int, float, double) an.

libsndfile bietet eine Lizensierung unter der GNU Lesser General Public License (LGPL) in Version 2.1 und Version 3 an.

LIBSNDFILE kann unter Ubuntu mit folgendem Befehl installiert werden:

$ sudo apt-get install libsndfile1-dev

Die Headerdateien sind unter Ubuntu im Verzeichnis

/usr/include

und die Libraries im Verzeichnis

/usr/lib

zu finden.

Link: http://www.mega-nerd.com/libsndfile/

OpenCV (Open Source Computer Vision Library)

OpenCV (OCV) ist eine umfassende Bibliothek um Probleme aus den Bereichen der Bildverarbeitung, Maschinenlernen, der linearen Algebra und vielen mehr lösen zu können. Die Bibliothek wird auch im Bereich der Robotik eingesetzt, um Mustererkennung durchzuführen. OCV ist in C/C++ geschrieben und auf Echtzeit-Anwendungen hin optimiert. Sie kann in C, C++, Java und Python eingebunden werden und unterstützt Linux, Windows, Mac und sogar Android.

Im Gegensatz zu FFTW und libsndfile ist OCV eine recht große API mit vielen Modulen, welche teilweise nur C oder C++-Schnittstellen oder beide anbietet. Sie ist also schwergewichtiger, kann aber durch weglassen von nicht benötigten Modulen zugeschnitten werden.

In diesem Artikel werde ich OCV verwenden, um die grafische Darstellung der Frequenzen als Bildabfolge zu erzeugen und in einer Videodatei zu speichern. OCV bietet primitive Operationen zum Zeichnen von gefüllten Polygonen an (mehr braucht es für den Artikel nicht) und kann AVI-Dateien erzeugen. Es mag andere geeignete Bibliotheken geben, die diese Funktionalitäten bieten, aber da ich mich viel OCV beschäftigt habe, fiel die Wahl auf diese Bibliothek.

OpenCV wird unter BSD lizensiert.

Die Installation von OCV ist leider nicht trivial und bietet viele Möglichkeiten der Anpassung. Auch kann OCV je nach Konfiguration von vielen anderen Bibliotheken abhängen. Für das Beispielprogramm ist FFMPEG Voraussetzung, damit die Videodatei erzeugt werden kann. Für eine detaillierte Installationsanweisung kann ich diese beiden Wiki-Links empfehlen: InstallGuide – OpenCV Wiki, FFMPEG – OpenCV Wiki.

Link: http://opencv.org

Sourcecode

Der Sourcecode ist auf der Plattform sourceforge.net gehostet und kann einfach als tar.gz-Archiv heruntergeladen werden.

Kompiliert wird das kleine Projekt durch Aufruf des make-Kommandos. Vorher müssen aber die Pfade zu den Header-Dateien von FFTW, libsndfile und OpenCV in der Datei makefile in Zeile 3 angepasst werden.

Ein Beispiel:

INCFLAGS=-I../OpenCV-2.4.3/modules/core/include -I../OpenCV-2.4.3/modules/highgui/include -I../OpenCV-2.4.3/modules/imgproc/include

Zu guter Letzt müssen auch die korrekten Pfade zu den Bibliotheken in Zeile 4 angegeben werden.

Ein Beispiel:

LDFLAGS=-L../OpenCV-2.4.3/debug/lib -L/usr/lib -lsndfile -lfftw3 -lopencv_core -lopencv_highgui

main.cpp

Der Monitor ist sehr einfach aufgebaut und besteht aus vier cpp-Dateien und drei Klassen (GraphEqualizer, DftUtils, Utils). Die main.cpp nimmt die Kommandozeilenparameter entgegen und wertet diese aus. Danach wird die Quell-Audiodatei eingelesen, die Videowriter-Klasse instanziiert und dann die Methode Utils::calcFrequencyBandHistograms(..) aufgerufen. Diese Methode führt die eigentliche Arbeit der Frequenzanalyse und der grafischen Aufbereitung durch.

// process commandline arguments
...
// read the audio input file and return a real number samples vector
Utils::readAudioFile(audioInputFile, &sfInfo, samples);

// create the video file for writing
VideoWriter vidWriter(videoOutputFile, CV_FOURCC('D', 'I', 'V', '3'), videoFps, Size(videoWidth, videoHeight), true);
if (!vidWriter.isOpened()) {
  throw runtime_error("VideoWriter failed to open output file!");
}

// iterate over the input sample vector and generate images of the frequency spectrum
Utils::calcFrequencyBandHistograms(samples, sfInfo, vidWriter, videoFps, videoWidth, videoHeight);
..

Leider ist das Einlesen von Stereo-Audiodateien zur Zeit noch etwas fehlerbehaftet. Ich konnte bis jetzt nicht ermitteln, ob der Fehler beim Einlesen über libsndfile oder später passiert. Der Fehler stellt sich so dar, dass das Eingangssignal immer zeitlich doppelt verzögert analysiert wird. Die Frequenzanalyse funktioniert aber dennoch korrekt. Nur die Eingangsdaten sind verzögert.

Utils::calcFrequencyBandHistograms

Am interessantesten ist die Methode calcFrequencyBandHistograms. Diese Methode iteriert über das Eingangssignal und errechnet für jeweils 512 Werte mit Hilfe der DFT ein Frequenzspektrum. Über dieses nun 257 Werte breite Spektrum wird ein Histogramm von 32 Werten gebildet. Dieses Histogramm wird dann normalisiert und über OpenCV in ein Bild gerendert. Dieses Bild wird dann der Instanz der Klasse VideoWriter übergeben.

Die Hauptschleife sieht wie folgt aus:

unsigned int numOfFpsSlices = <number of slices to analyze>
vector samples = <audio input vector of real numbers>

for (unsigned int sliceIdx = 0; sliceIdx < numOfFpsSlices; sliceIdx++) {
  // DFT generates vector of magnitudes (length = (dftSize / 2) + 1 )
  vector magVec = DftUtils::calcDftWithOffset(samples, sliceIdx * numOfSamplesPerFrame, dftSize);

  // generate a normalised histogram from the calced frequency magnitudes
  generateHistogram(magVec, magVecLenToDisplay, magnitudeHistogram);

  // generate the graphical frequency representation and store the frame in the video stream
  vidWriter << graphEqualizer.generateEqualizerImageFromHistogram(magnitudeHistogram);
}

Auf ein wichtiges Detail in Bezug auf den Aufruf der Generierung des Histogramms möchte ich noch eingehen. Die Beispiele und Tests habe ich mit Musikstücken durchgeführt. Dabei stellte sich heraus, dass Amplitudenveränderungen im Frequenzspektrum im Bereich von 0 – 6000 Hz stattfinden – meistens sogar nur von 0 bis 4000 Hz. Deshalb ignoriert das Histogramm alle Frequenzen welche höher als 4000 Hz liegen. Dieser Grenzwert kann aber einfach über eine Konstante in der Methode calcFrequencyBandHistograms angepasst werden.

DftUtils::calcDftWithOffset

Dank FFTW ist die Anwendung der DFT so einfach gestaltet, so dass die Methode calcDftWithOffset sehr unspektakulär aussieht. Zuerst wird Speicher für die Aufnahme der Eingangswerte mit Hilfe der FFTW-Methoden angelegt. Dies ist notwendig, damit FFTW den Zugriff auf den Speicher optimieren kann.

unsigned int N = numOfSamples;
double *dftIn = fftw_alloc_real(N);

unsigned int outN = N / 2 + 1;
fftw_complex* dftOut = fftw_alloc_complex(outN);

Als nächstes wird ein sogenannter Plan angelegt. Dieser Plan beschreibt die durchzuführende Operation und die Länge der Eingangswerte. FFTW bietet auch Methoden an, welche statt reeller Zahlen, komplexe Zahlen als Eingangswerte entgegennehmen. Auch die Ausführungsgeschwindigkeit der DFT kann optimiert werden, indem statt FFTW_ESTIMATE zum Beispiel FFTW_MEASURE verwendet wird (mehr Details dazu hier).

Als weitere Optimierungsmaßnahme sollte ein Plan nur einmal erzeugt und dann wiederverwendet werden. Dies ist aber nur dann möglich, wenn sich die Anzahl der Eingangswerte und dessen Wertetyp (reell, komplex) nicht ändert. In dem Beispielcode ist die Performance aber unwichtig, so dass deshalb der Plan bei jedem Aufruf erzeugt und am Ende zerstört wird:

fftw_plan dftPlan = fftw_plan_dft_r2c_1d(N, dftIn, dftOut, FFTW_ESTIMATE);

Um die Leckeffekte zu minimieren wird als nächstes die Hanning-Fensterfunktion auf die Eingangswerte angewandt. Für Testzwecke ist es möglich per #define USE_HANNING_WINDOW die Fensterfunktion auszukommentieren um die Auswirkungen auf das Resultat der DFT sehen zu können.

#define USE_HANNING_WINDOW
#ifdef USE_HANNING_WINDOW    
    applyHanningWindow(samples, startOffset, N, dftIn);
#else    
    for (unsigned int i=0;i<N;i++) {
        dftIn[i] = samples[startOffset + i];
    }
#endif

Als nächstes wird die DFT über den vorher erzeugten Plan ausgeführt und dann aus dem Vektor komplexer Ausgangswerte ein Vektor von Frequenzbeträgen generiert.

fftw_execute(dftPlan);

vector magVec(outN);
for (unsigned int i=0;i<outN;i++) {

    double real = dftOut[i][0];
    double imag = dftOut[i][1];
    double mag = sqrt((real * real) + (imag * imag));
        
    magVec[i] = mag;
}

Zu guter Letzt wird der Plan zerstört und der benutze Speicher zurückgegeben:

fftw_destroy_plan(dftPlan);

fftw_free(dftIn);
fftw_free(dftOut);

Utils::generateHistogram

Um das per DFT errechnete Frequenzspektrum darstellen zu können, wird ein Histogramm erstellt. Dazu wird ein Histogramm mit einer bestimmten Größe (z.B. 32) definiert und die Werte der jeweiligen Frequenzanteile zu dem entsprechenden Index addiert. Nachdem das vollständige Frequenzspektrum akkumuliert wurde, wird es für die grafische Darstellung normalisiert. Der Normalisierungs-Faktor ist unter anderen von der Breite der DFT und der Parameter der Fensterfunktion abhängig.

Normalerweise würde jeder Histogrammwert normalisiert werden, indem er durch den größten im Histogramm vorkommenden Wert dividiert wird. So wäre es aber nicht möglich, leise Passagen im Eingangssignal von lauten unterschiedlich darstellen zu können. Deshalb wird als Kompromiss ein fester Divisor verwendet (hier: 30), welcher durch Testen verschiedener Musikstücke ermittelt wurde.

// accumulate histogram of magnitudes for each frequency band
std::fill(histogram.begin(), histogram.end(), 0.0);
for (unsigned int i=0;i<magVecSize;i++) {
    int idx = i * ((double)histogram.size() / magVecSize);
    histogram[idx] += fabs(magVec[i]);
}

for (unsigned int i=0;i<histogram.size();i++) {
    histogram[i] /= 30.0;
}

GraphEqualizer::generateEqualizerImageFromHistogram

Die grafische Darstellung des Histogramms übernimmt die Methode generateEqualizerImageFromHistogram. Unter Verwendung der Methode fillPoly von OpenCV wird für jeden Wert des Histogramms eine gefärbte „Säule“ dargestellt.

Anhand von mehreren konstanten Farbwerten wird per lineare Interpolation ein Farbwert der jeweiligen Säule ermittelt und diese damit gefüllt. Die Säule hat ihren Y-Ursprung jeweils bei der Hälfte der Höhe des Bildes und dehnt sich nach oben und unten mit gleicher Länge aus.

Um einen schöneren Effekt des Frequenzverlaufs zu erzielen, wird zwischen hohen und darauf folgenden niedrigen Frequenzwerten nicht sofort gewechselt, sondern sich zeitverzögert vom höheren Frequenzwert zum niedrigen angenähert.

// generate an empty (black colored) image
Mat image = Mat::zeros(imageHeight, imageWidth, CV_8UC3);
..    
Point polyPntArr[1][4];
const Point * polyPnt[1] = {polyPntArr[0]};
const int numPoints[] = {4};

for (unsigned int i=0;i<histogram.size();i++) {
    double barHeight;
    if (histogram[i] > 0.05) {
        barHeight = maxBarHeight * histogram[i];
    } else {
        barHeight = minBarHeight;
    }
    barHeight /= 2.0;
        
    // smooth the bar height based on the previous bar height
    barHeights[i] -= 4.0;
    if (barHeights[i] < barHeight) {
        barHeights[i] = barHeight;
    } else {
        barHeight = barHeights[i];
    }

    // calc the bar color
    ..
    const Scalar& predColor = barColors[predColorIdx];
    const Scalar& succColor = barColors[predColorIdx + 1];

    int diffR = (predColor[0] - succColor[0]) / histogram.size();
    int diffG = (predColor[1] - succColor[1]) / histogram.size();
    int diffB = (predColor[2] - succColor[2]) / histogram.size();

    Scalar col(predColor[2] + diffB, predColor[1] + diffG, predColor[0] + diffR);

    // draw filled polygon bar
    double x0 = (i * barWidth) + (barWidth / 2.0) - finalBarHalfWidth;

    polyPntArr[0][0] = Point(x0, halfHeight - barHeight);
    polyPntArr[0][1] = Point(x0, halfHeight + barHeight);
    polyPntArr[0][2] = Point(x0 + finalBarWidth, halfHeight + barHeight);
    polyPntArr[0][3] = Point(x0 + finalBarWidth, halfHeight - barHeight);

    fillPoly(image, polyPnt, numPoints, 1, col);
}

Bücher und URLs zum Thema „digitale Signalverarbeitung“

Tagged with: , , , , ,
Veröffentlicht in Digitale Signalverarbeitung

Kommentar verfassen

Trage deine Daten unten ein oder klicke ein Icon um dich einzuloggen:

WordPress.com-Logo

Du kommentierst mit Deinem WordPress.com-Konto. Abmelden / Ändern )

Twitter-Bild

Du kommentierst mit Deinem Twitter-Konto. Abmelden / Ändern )

Facebook-Foto

Du kommentierst mit Deinem Facebook-Konto. Abmelden / Ändern )

Google+ Foto

Du kommentierst mit Deinem Google+-Konto. Abmelden / Ändern )

Verbinde mit %s

%d Bloggern gefällt das: