mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Verständnisfrage Fouriertrafo


Autor: Markus R. (maggggus)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo!

Ich habe eine Frage zur Fouriertrafo. Ich habe ein C++ Programm 
geschrieben (läuft auf dem PC) das eine (mono) WAVE-Datei einliest und 
auf die Werte die FFT anwendet. Das ganze funktioniert auch, ich bekomme 
die (komplexen) Fourierkoeffizienten berechnet, das Ergebnis sieht in 
etwaa so aus wie in der angehängten Datei zu sehen ist.

Was muss ich nun machen, um ein Histogramm mit den einzelnen Frequenzen 
zu bekommen? Wenn meine Eingabedatei aus 130000 Abtastpunkten besteht, 
dann bekomme ich 130000 komplexe Zahlen als Ergebnis.

Was ich aber haben möchte, ist für Frequenz von 1 Hz bis 20000 Hz wie 
oft sie vorkommt.

Autor: mki (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich denke was du braust ist so was wie eine Kurzzeit FFT. Lies die 
ersten 40000 Werte aus und wende darauf die FFT an. Sie liefert dir dann 
die Frequenzen von -20000 bis +20000 - aber Vorsicht, nicht unbedingt in 
dieser Reihenfolge. Davon dann den Betrag und du hast schon ein 
Histogramm. Das gleiche machst du dann mit den nächsten 40000 Werten und 
addierst das dann zu deinem Teilhistogramm hinzu. Ich denke so sollte es 
gehen.

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Danke! Das hilft mir weiter. Wenn du mir mit der Reihenfolge noch helfen 
kannst, hast du mein Problem gelöst!

Autor: mki (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Das mit der Reihenfolge ist glaub ich je nach Implementation 
unterschiedlich. Bei der vorgefertigten FFT von Matlab ist es glaub in 
etwa so - zuerst die Werte von Null bis 20000, danach von -20000 bis -1. 
Frag mich jetzt nicht wieso es gerade so, das hat vielleicht 
irgentwelche Geschwindigkeitsvorteile. Ich vermute wenigstens mal, dass 
es bei dieser FFT auch so sein wird.

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich verwende die Implementierung "FFTW" von fftw.org. Ich kann mal 
nachschauen wie es da ist, bin mir aber nicht sicher, ob ich das finde. 
Wenn du dir bei MATLAB sicher bist, kann ich es ja mal gegenchecken und 
so zumindest sehen, ob es bei FFTW so ist wie bei MATLAB.

Kann ich dann die Werte von -20000 bis -1 einfach wegwerfen?

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Noch was: Wenn ich dann die 40000 Werte habe, sind die immer noch 
komplex. Wie bekomme ich dann reellwertige Zahlen daraus? Imaginärteil 
wegschmeissen?

Autor: mki (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Du nimmst den Betrag von jedem Wert. Dann sind die Werte alle reell. Und 
was die Reihenfolge angeht: Matlab benutzt auch die FFT von fftw. Ich 
habe mir nochmal die Reihenfolge angeschaut und es kommen zuerst die 
Frequenzen eins bis 20000 und danach -20000 bis Null. Die negativen 
Frequenzen kannst du denk ich mal einfach weglassen.

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Danke erst mal!

OK Lasse alle negativen Frequenzen weg. Mit Betrag meinst du in MATLAB 
die abs() Funktion oder? Also die für a= 1+i dann abs(a) = 1.4142 
liefert.

Noch eine Frage: Ich möchte ja die ganze Wave-Datei analysieren, und 
deren Abtastpunkte sind ja nichtimmer ein ganzzahliges Vielfaches von 
40000. Soll ich dann bei der letzten Anwendung der FFT die Daten einfach 
mit Nullen auffüllen, so dass ich wieder auf 40000 Abtastpunkte für die 
Eingabe der FFT komme?
Oder verändere ich dann das Signal und bekomme ein anderes Ergebnis?

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich habe jetzt mal mit Audacity eine Datei generiert, die aus einem 220 
Hz Sinuston besteht.

Ich erwarte in meinem Frequenzdiagramm also einen Peak bei 220, die 
anderen Werte sollen also 0 sein.

Was ich kriege ist schon fast das, nur das der Peak bei genau 400 liegt!


Liegt das daran, das ich keine Fensterfunktion (siehe 
http://de.wikipedia.org/wiki/Fensterfunktion ) verwende?

Autor: mki (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ja ich meine die abs() Funktion. Und der Rest wird einfach mit Nullen 
aufgefüllt. Was aber noch wichtig ist: Ich hoffe du beachtest die 
Abtastrate. Die muss zwingend bei 40000Hz liegen. Das Stichwort dazu 
lautet Nyquist-Shannon-Abtasttheorem. Ich hoffe ich stell dich damit 
nicht vor noch mehr Problemen.

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Abtasttheorem: Wenn ich den Bereich von 0-20kHz Untersuchen will, muss 
ich mit MINDESTENS 40kHz abtasten, oder (also nicht unbedingt exakt 
40000)? Ich habe bei meinen Versuchen bis jetzt immer 44100 Hz 
eingestellt gehabt.

Liegts jetzt an der Fensterfunktion?

Autor: mki (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
OK Abtasttheorem ist beachtet. Daran kann es nicht liegen. Wieso der 
Peak aber bei 400 ist das weiss ich jetzt auch nicht, da bin ich 
momentan auch überfragt. Ich kann mir nur nicht vorstellen, dass es an 
der Fensterfunktion liegt. Die kannst du glaub ich ausschliessen.

Autor: mki (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hier steht was wie du die transormierte Funktion in Frequenzen umrechnen 
kannst. Das liegt an der fft von Matlab wieso der Peak bei 400 liegt.

http://www.phys.ufl.edu/~aaronm/fft.html

Und vielleicht brauchst du auch gar nicht mal die 40000 Punkte 
ausschneiden und darauf die fft machen. Wenn du die fft auf den ganzen 
Vektor machst, bekommst du zwar 130000 Werte, aber die Frequenzen größer 
20000 sind dann alle Null. Du kannst dann also später ausschneiden.

Ich hoffe das hilft dir alles weiter.

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Danke für deine Antwort, deine Beiträge sind sehr hilfreich!

Habe mir angegebenen Link durchgelesen, was muss ich für D einsetzen? 
Das habe ich nicht verstanden. 1/44100 , also die Zeit die zwischen zwei 
Abtastpunkten verstreicht?

Autor: mki (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich denke eher nicht. Die Zahlen werden ja sonst ziemlich klein - 
jedenfalls nicht größer als eins. Ich denke es ist eher so eine Art 
Strecke. Versuch es mal mit der Strecke, die eine Schallwelle innerhalb 
zwei t-Werten zurücklegt. Also so was wie c/44100 
(c=Lichtgeschwindigkeit).
Ansonsten weiss ich es jetzt auch nicht.

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hm, die Zahlen werden zwar klein, das ist aber eigentlich egal. Es 
werden ja alle Werte gleich skaliert. Ich kann auch 330/44100 nehmen 
(Schallgeschwindigkeit durch Abtastrate), dann werden sie halt etwas 
größer.


Habe aber grade festgestellt, dass plot(eingabedaten) keinen schönen 
Sinus liefert. Oder liegt das an der Abtastung?

Ich werde mal versuchen die Eingangsdaten direkt in Matlab zu erzeugen 
(da muss ich aber erst heruasfinden wie man in Matlab ein Signal 
abtastet)...

Autor: mki (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Das kein schöner sinus herauskommt liegt daran, dass Matlab zwei Punkte 
einfach verbindet.

Das abgetastete Signal kannst du in Matlab folgendermaßen erzeugen:
a=sin(2*pi/44100*220*[1:130000]);

Naja die Nulldurchgänge sind hier recht dicht beieinander. Man sieht 
nicht wirklich was. Probier mal mit größeren Frequenzen.

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hm.

Also hier mal meine Matlab Befehle:

a=sin(2*pi/44100*220*[1:130000]);
P = abs(fft(a));
L=130000;
f = (330/44100)*P(1:L/2)/L;
plot(f);


Ich habe in der vorletzten Zeile das P ergänzt, das fehlt auf der Seite, 
die du vorhin verlinkt hast. Außerdem habe statt 0:L/2 jetzt 1:L/2 
geschrieben, weil MATLAB sonst meckert. Das Problem bleibt: Auf dem 
Graphen ist der Peak bei ca. 650 und nicht bei 220 wie ich erwarten 
würde.
Außerdem skaliert die vorletzte Zeile die Werte ja nur (d.h. streckt 
oder staucht sie), verschiebt sie jedoch nicht nach links oder rechts.


Mit dem doppelt logarithmischen Plot kann ich gar nichts anfangen...

Autor: mki (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ja OK stimmt auch wieder. Bin momentan auch überfragt. Vielleicht später 
noch. Schreib mal rein wenn du es geschafft hast, würde mich jetzt auch 
interessieren.

Autor: Andreas Schwarz (andreas) (Admin) Benutzerseite Flattr this
Datum:

Bewertung
0 lesenswert
nicht lesenswert
N = 1024                  % Länge der FFT
fs = 44100                % Abtastrate

a=sin(2*pi/fs*220*[1:130000]);

P = abs(fft(a,N));
f = linspace(0,fs/2,N/2);  % Punkt N/2 entspricht der Frequenz fs/2
plot(f, P(1:N/2))

Mit Schallgeschwindigkeit oder Histogramm hat das alles nichts zu tun.

Wenn du ein "reales", stochastisches Signal (Musik, Sprachsignal) 
analysieren willst, dann musst du die FFTs von mehreren (gefensterten) 
Stücken des Signals mitteln. Schau dir die Hilfe zur Funktion pwelch an, 
die zeigt wie's gemacht wird.

Autor: Rüdiger Knörig (sleipnir)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Die N FFT-Koeffizienten decken den Frequenzbereich von 0 bis zur 
Abtastfrequenz f_T ab. Somit ist die Zuordnung FFT-Index k -> Frequenz f 
durch die folgende Verhältnisgleichung gegeben:
f/f_T= k/N
Daraus folgt
f=k*f_T/N

Zu beachten ist dabei, daß gemäß dem Abtasttheorem ab der Frequenz f_T/2 
die erste spektrale Wiederholung beginnt.
Ich habe das damals mal für meine Studenten in einer Kurzanleitung 
zusammengefaßt:
http://www.nue.tu-berlin.de/teaching/praktika/mult...

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Erstmal vielen Dank für die fundierten Beiträge aller 3 Poster!

Also der Code von Andreas hat einwandfrei funktioniert. Der Peak war bei 
ca. 216 Hz, als ich N = 16384 gesetzt habe, war er dann exakt bei 220 
Hz.

Ich habe dann Abtastpunkte aus meiner 220-Hz-Sinuston Wavedatei 
importiert, dort war der Peak bei 440 Hz zu sehen. Ich habe dann im 
linspace und im plot-Befehl das N/2 durch N ersetzt und jetzt passts 
(mir gehts nicht darum ob es theoretisch richtig ist, wenn ich in der 
Praxis die dominanteste Frequenz ermitteln kann reicht mir das).

Außerdem habe ich mal mit einem "echten" Signal getestet 
(Akustik-Gitarrensaite mit 220 Hz, per Mic aufgenommen), auch hier sehe 
ich schön die Peaks da wo sie sein sollen (wenn N hoch ist und N/2 durch 
N ersetzt wurde).

Ich werde die nächsten Tage versuchen, das ganze rein in C/C++ zu 
übersetzen.


Vielen Dank nochmal an alle Helfer!

Autor: Markus R. (maggggus)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Inzwischen sind einige Tage vergangen und mein C++ Programm läuft. Ich 
habe festgestellt, dass ich recht viel Speicher brauche, denn N muss 
mind. 16384 sein.

Eigentlich interessieren mich aber nur die Frequenzen aus dem Bereich 
50-500 Hz. Was muss ich machen, damit ich nur diesen Bereich analysiere, 
und so mit dem N runterkomme?

Autor: Andreas Schwarz (andreas) (Admin) Benutzerseite Flattr this
Datum:

Bewertung
0 lesenswert
nicht lesenswert
N = 1024                  % Länge der FFT
fs = 44100                % Abtastrate
dec_fact = 30;            % dezimierungs-faktor

fs_dec = fs / dec_fact;

a = sin(2*pi/fs*220*[1:130000]);
a_dec = decimate(a,dec_fact);

P_dec = abs(fft(a_dec,N));
f_dec = linspace(0,fs_dec/2,N/2);  % Punkt N/2 entspricht der Frequenz fs_dec/2
plot(f_dec, P_dec(1:N/2))

Autor: Markus R. (maggggus)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Danke für deine Antwort!

Wenn ich das recht verstehe, basiert das darauf, dass man einfach mit 
der Abtastfrequenz weiter runtergeht.

Ich habe jetzt mal selber etwas rumprobiert und bin zu folgendem Code 
gekommen:
fid = fopen('gitarre-tonA-kurz.wav','r'); fread(fid, 44); samples = fread(fid, 1048576, 'int16');
N = 256; fs=44100; dec_fact = 55; fs_dec=fs/dec_fact; samples_dec = decimate(samples, dec_fact);
P_dec = abs(fft(samples_dec, N)); f_dec = linspace(0,fs_dec/2,N/2); plot(f_dec, P_dec(1:N/2));

Die Datei habe ich mal angehängt, es ist der Ton A einer Gitarrensaite. 
Peaks erwarte ich also bei 55 und 110 und 220 Hz.

Das funktioniert recht gut. Die FFT läuft also nur mehr über 256 Punkte, 
statt wie vorher über 16384 Punkte, und die Peaks sind noch schön zu 
erkennen.


Was könnte ich noch tun, um den Speicherbedarf zu drücken?
Ich brauche ich immer noch 256 Werte für den Eingabevektor und 256 Werte 
für den Ausgabevektor.

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.