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.
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.
Danke! Das hilft mir weiter. Wenn du mir mit der Reihenfolge noch helfen kannst, hast du mein Problem gelöst!
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.
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?
Noch was: Wenn ich dann die 40000 Werte habe, sind die immer noch komplex. Wie bekomme ich dann reellwertige Zahlen daraus? Imaginärteil wegschmeissen?
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.
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?
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?
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.
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?
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.
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.
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?
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.
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)...
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.
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...
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.
1 | N = 1024 % Länge der FFT |
2 | fs = 44100 % Abtastrate |
3 | |
4 | a=sin(2*pi/fs*220*[1:130000]); |
5 | |
6 | P = abs(fft(a,N)); |
7 | f = linspace(0,fs/2,N/2); % Punkt N/2 entspricht der Frequenz fs/2 |
8 | 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.
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/multimedia/OctaveTut.pdf
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!
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?
1 | N = 1024 % Länge der FFT |
2 | fs = 44100 % Abtastrate |
3 | dec_fact = 30; % dezimierungs-faktor |
4 | |
5 | fs_dec = fs / dec_fact; |
6 | |
7 | a = sin(2*pi/fs*220*[1:130000]); |
8 | a_dec = decimate(a,dec_fact); |
9 | |
10 | P_dec = abs(fft(a_dec,N)); |
11 | f_dec = linspace(0,fs_dec/2,N/2); % Punkt N/2 entspricht der Frequenz fs_dec/2 |
12 | plot(f_dec, P_dec(1:N/2)) |
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:
1 | fid = fopen('gitarre-tonA-kurz.wav','r'); fread(fid, 44); samples = fread(fid, 1048576, 'int16');
|
2 | N = 256; fs=44100; dec_fact = 55; fs_dec=fs/dec_fact; samples_dec = decimate(samples, dec_fact); |
3 | 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.
Thread beobachten
Seitenaufteilung abschaltenBitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.