Hallo zusammen, erstmal: erfolgreiches und gesundes 2008 an alle ! folgendes: Ich sample ein Sprachsignal (nicht periodisch) von der Mixerline MIC und erhalte damit ein wert- und zeitdiskretes Signal - ein PCM-Signal. U(k*Ta) = [100 200 150 100 150 100 200 ...] FRAGE: Wie kann ich eine FFT (real) mit C oder C++ implementieren, die aus dieser Folge das Spektrum bildet ?
Was willst du mit dem Spektrum des PCM-Signals? Oder doch das Spektrum des Sprachsignals?
Hallo Hägar, ich möchte natürlich das Spektrum des Sprachsignals - zur grafischen Darstellung. Jedoch steht mir ja nur das abgetastete Sprachsignal, also das PCM Signal, zur verfügung. Also: wie kann ich das Spektrum einer wert- und zeitdiskreten Folge berechnen ? R-FFT: F(i) = SUM[k=0 TO N-1] { u(k) * exp(-j*i*k*2PI/N) } // i = 0...N-1 komme aber irgendwie nicht klar mit dieser Formel. =( eric.
Muss es denn unbedingt C(++) sein? Matlab/Octave/Scilab können das off-the-shelf. Für C gibts auch massig Fertiges (FFTW), aber die Darstellung ist eben schwieriger. Die von dir dargestellte Formel ist eine DFT. Macht außer enormer Rechenzeit aber keinen Unterschied. Schließlich wären da noch fertige Oszi/Spectrum Analyzer Programme.
Karl, leider muss es C(++) sein. Das Programm wird für einen ARM9 Prozesser entwickelt. Ich nutze die WAVE API (win32) um mit waveInStart() die samples in einem Buffer zu schreiben. Nun kann ich diese natürlich grafisch darstellen. Jedoch ist mir das Spektrum lieber. Habe nun auch die Fastest Fourier Transformation in the West entdeckt und werde versuchen eine LIB für Windows CE 6 daraus zu erzeugen. Aber unabhängig davon würde mich interessieren, wie man (bzw. Du) die DFT-Formel von oben interpretierst. Eigentlich erwarte ich: F(i): . . | | . | _|____|____|__...____>Hz ABER : in der DFT-Formel ist keine Frequenz vorhanden. Verstehe die Formel nicht wirklich ! Erkläre mal bitte =) Thx, eric
Die Grundfrequenz in der Formel oben ist f = Sample rate / N Die abgetasteten Frequenzen sind damit f(i) = f * i mit i = [0 .. N / 2] oben ist [0 .. N - 1] angegeben... kommt wir aber falsch vor. Hoffe, ich habe mich jetzt auf die Schnelle nicht vertan - die Schreibweise der obigen Formel ist nicht wirklich toll zu lesen - zumal man raten muss, welche Variable für was steht. Was genau verstehst du nicht an der Formel?
> oben ist [0 .. N - 1] angegeben... kommt wir aber falsch vor.
Ist schon richtig. Bringt nur nix bis N-1 das Spektrum zu berechnen ab
N/2 wiederholen sich die Werte.
Die Formel für die DFT ergibt sich, wenn man die Fourier Transformation (zeit)diskretisiert, sprich das kontinuierliche Signal mit einer Summe aus äquidistanten Impulsen multipliziert. Dabei wird das Integral zur Summe. Wie schon geschrieben, ist der Abstand zwischen den einzelnen "bins" Fs/N, wobei sich das Spektrum alle Fs wiederholt. Der Bezug zu den Frequenzanteilen wird also erst durch Fs gegeben. Ich weiß schon, dass das eine sehr hässliche Formel ist, am besten ein paar Beispiele von Hand rechnen, wenn einem der Weg nicht klar ist.
Kai, für besseres (??) Verständnis die Formel: http://www.kgw.tu-berlin.de/statisch/lehre/skript/ds/node35.html zu deiner Erklärung: ich taste mit 11025 Hz ab. Und zwar eine sec lang. D.h: SampleRate = 11025; //11.025 kHz N = 11025; //Samples per second damit: [quote] f = Sample rate / N f(i) = f * i mit i = [0 .. N / 2] [/quote] f = 1 und f(i) von 0 Hz bis 5512 Hz WAS ICH NICHT VERSTEHE : wie multipliziere ich den Koeffizienten e^(..) ? bzw. wie berechne ich F(2) bei u(k) = [ 1 5 10 ] //N=3 ?? Thx, eric.
Ach ja: Wenn du die FFTW einsetzt, musst du dein Programm entweder unter der GPL veröffentlichen oder dafür zahlen!
>WAS ICH NICHT VERSTEHE : wie multipliziere ich den Koeffizienten >e^(..) ? >bzw. wie berechne ich F(2) bei u(k) = [ 1 5 10 ] //N=3 ?? Err, e^(jx) = cos(x) + jsin(x) Was gibbet da nicht zu verstehen? F(2) = 1*e^(-j2pi*2*0/3) + 5*e^(-j2pi*2*1/3) + 10*e^(-j2pi*2*2/3)
Hier findest du die ganze hässliche Wahrheit ;) http://de.wikipedia.org/wiki/Komplexe_Zahlen#Potenzen Schonmal mit komplexen Zahlen gearbeitet? >bzw. wie berechne ich F(2) bei u(k) = [ 1 5 10 ] //N=3 ?? Genauso wie für alle anderen i - oder wolltest du das jetzt mit konkreten Zahlen sehen? @Karl >Was gibbet da nicht zu verstehen? Für jemanden, der mit der komplexen e-Funktion noch nichts gemacht hat ist das ziemlich harter Tobac finde ich.
>Für jemanden, der mit der komplexen e-Funktion noch nichts gemacht hat >ist das ziemlich harter Tobac finde ich. Stimmt allerdings. Dann ist aber eine DFT oder gar FFT noch härter. Deshalb habe ich daran jetz garnicht gedacht, bei uns wurden komplexe Zahlen und einige ihrer Auswüchse schon im 1. Semester abgehandelt.
Karl wrote: >>WAS ICH NICHT VERSTEHE : wie multipliziere ich den Koeffizienten >>e^(..) ? >>bzw. wie berechne ich F(2) bei u(k) = [ 1 5 10 ] //N=3 ?? > Err, e^(jx) = cos(x) + jsin(x) bekannt. > Was gibbet da nicht zu verstehen? > > F(2) = 1*e^(-j2pi*2*0/3) + 5*e^(-j2pi*2*1/3) + 10*e^(-j2pi*2*2/3) Sorry, aber ist schon etwas her, dass ich mit komplexen Ausdrücken gerechnet habe. Wie ist der skalar bei 5*e^(-j2pi*2*1/3) ?? SORRY !
Wenn dir e^(jx) = cos(x) + jsin(x) bekannt ist, dann müsstest du das doch einfach nur noch einsetzen?! 5*e^(-j2pi*2*1/3) => x= -2pi*2*1/3 = -4,1888 cos(-4,1888) + j * sin(-4,1888) = (0,997; -0,073) oder? Mit Skalaren ist da nichts, das Ergebnis ist natürlich auch komplex. ach ja, *5 nicht vergessen: 5 * (0,997; -0,073) = (4,985; -0,365)
5*e^(-j*4/3*pi) = 5*(cos(4/3*pi) + j*sin(4/3*pi) = 5*(-1/2-j*sqrt(3)/2) Zur Berechnung von F(i) musst du immer über alle N Punkte summieren. Die Leistung auf der jeweiligen Frequenz ist dann der Betrag von F(i), die Phase ist arctan(Im(F(i))/Re(F(i)). exp(-j*i*k*2*pi) kannst du dir als komplexen Drehzeiger vorstellen. Also ein "Vektor", wenn du so willst, in der komplexen Ebene Ursprung in 0,0 und Ende auf einem Kreis mit Radius 1 um den Ursprung. Dieser rotiert mit i und k um den Ursprung, diese Drehung ändert die Phase, aber nicht den Betrag. Der Betrag wäre in deinem Beispiel 5. Denn sin^2(x)+cos^2(x) = 1!! Wie du siehst ist das eine Menge Rechnerei. Die üblichen DFT/FFT Längen liegen bei 128, 256, 512 oder noch höher. Deshalb hat sich die FFT entwickelt, dabei werden bestimmte Symmetrien von sin und cos ausgenutzt und es müssen nicht alle Werte berechnet werden. Ansonsten kannst du dir das alles mit Matlab berechnen lassen. Falls du kein Matlab hast der freie Klone Octave beherrscht das auch.
Leute: Vielen Dank ! jetzt hats klick gemacht. Weiß wieder wie es funktioniert. Danke euch wirklich ! Hat von euch schon mal jemand die FFTW verwendet oder Erfahrung mit FFT und C(++) ?? =)
Hi, habe jetzt eine LIB für WinCE 6 aus der KISS FFT gebaut. Doch die ersten Tests verwundern mich etwas !
1 | #include "kiss_fftr.h" |
2 | #include "_kiss_fft_guts.h" |
3 | |
4 | #define kiss_fft_scalar BYTE
|
5 | |
6 | #define N 4
|
7 | |
8 | kiss_fft_scalar rin[N]; |
9 | kiss_fft_cpx sout[N]; |
10 | |
11 | int main(void) |
12 | {
|
13 | /* rin = {0 1 2 3} */
|
14 | for(BYTE i=0; i<N; i++) |
15 | {
|
16 | rin[i] = (kiss_fft_scalar) i; |
17 | }
|
18 | kiss_fftr_cfg kiss_fftr_state = kiss_fftr_alloc(N,0,0,0); |
19 | kiss_fftr(kiss_fftr_state,rin,sout); |
20 | for(int aa=0; aa<N; aa++) |
21 | {
|
22 | _tprintf(L"rin[%u] = %u\t sout: r = %u i = |
23 | %u\r\n",aa,(BYTE)rin[aa],(BYTE)sout[aa].r,(BYTE)sout[aa].i); |
24 | }
|
25 | }
|
AUSGABE : u(k) = rin = [ 0 1 2 3 ] rin[0] = 0 sout: r = 0 i = 1 rin[1] = 1 sout: r = 2 i = 3 rin[2] = 2 sout: r = 0 i = 0 rin[3] = 3 sout: r = 0 i = 0 ERWARTET : für F(1) = sout[1] : F(1) = 0*e^(...) + 1*e^(-j*2*1*PI/4) + 2*e^(-j*2*2*PI/4) + 3*e^(-j*2*3*PI/4) = 0 + 1*[cos(PI/2) + j*sin(PI/2)] + 2*[cos(PI) + j*sin(PI)] + 3*[cos(3PI/2) + j*sin(3PI/2)] = -2 -2j ABER : -2 -2j != 2 +3j Kann mir jemand auf die Sprünge helfen ? Thx, eric
Nope... außer, dass ich BYTE nicht weit traue, was trigonometrische Funktionen angeht ;) Zum Testen der generellen Funktionen empfehle ich float. Wenn die Funktionen damit funktionieren, kann man die Zahlenformate immer noch anpassen und hat dann gute Vergleichswerte. Wenn die Bibliothek damit allerdings nicht um kann, bringt dir das nicht viel. Ich kenne die ganzen kiss-Dinger nicht. Wenn du das selber rechnen würdest, könnte ich dir wahrscheinlich eher helfen.
Verstehe ich das richtig, dass du das PCM-Signal abtastet und erwartest, das dies dem Spektrum des Audiosignal entspricht?
Also, wenn ich das jetzt mal zu Fuß durchrechne: PCM: 0 1 2 3 Sinusüberlagerung: 0 1 0 -1 Cosinusüberlagerung: 1 0 -1 0 Sinussumme (Imaginärteil): 0*0 + 1*1 + 2*0 + 3*-1 = -2 Cosinussumme (Realteil): 0*1 + 1*0 + 2*-1 + 3*0 = -3 macht: -2 + j*-2 Du befindest dich aber schon auf der obersten Grenzfrequenz der DFT - da gilt glaube ich noch eine Sonderregel, was die anschließende Verrechnung mit Koeffizienten angeht. ps: habe die Rechnung nochmal korrigiert.
> Nope... außer, dass ich BYTE nicht weit traue, was trigonometrische > Funktionen angeht ;) > Zum Testen der generellen Funktionen empfehle ich float. BYTE verwende ich hier nur zum testen. Die Wertefolge aus dem Beispiel führt zu ganzzahligen Ergebnissen für Real und Imaginärteil. Argumente -> PI/2, PI, 3PI/2 Werde aber mit float weiter testen. > Verstehe ich das richtig, dass du das PCM-Signal abtastet und erwartest, > das dies dem Spektrum des Audiosignal entspricht? Nein. Das abgetastete Audiosignal soll durch Transformation in das entsprechende Spektrum überführt werden. > macht: -2 + j*-2 > Du befindest dich aber schon auf der obersten Grenzfrequenz der DFT - da > gilt glaube ich noch eine Sonderregel, was die anschließende Verrechnung > mit Koeffizienten angeht. Dieses Ergebnis bekomme ich auch (siehe Berechnung letzter Beitrag), wenn ich von Hand rechne. Die KISS-FFT ergibt jedoch 2 +3j für F(1). Könnte vielleicht an BYTE liegen, denke ich aber eigentlich nicht wirklich. Ein weitere Test mit FLOAT wird hier Klarheit schaffen. @Kai: > ps: habe die Rechnung nochmal korrigiert. aber nur das Ergebnis (-2 + j*-2) -> siehe Herleitung bzw. Rechenweg obersten Grenzfrequenz der DFT --> wie kommste drauf ? und Verwendet du die KISS FFT ? Kannst du das Beispiel von mir mal bei Dir laufen lassen ? Dabei kannst Du ja auch dein Ergebnis überprüfen =) Thx, eric
>Werde aber mit float weiter testen. Mach das, im schlimmsten Falle kommt dasselbe raus ;) >Dieses Ergebnis bekomme ich auch (siehe Berechnung letzter Beitrag), Ich weiß, wollte nur von meiner Seite aus nochmal bestätigen, dass ich auf dasselbe Ergebnis wie du komme. >> ps: habe die Rechnung nochmal korrigiert. >aber nur das Ergebnis (-2 + j*-2) -> siehe Herleitung bzw. Rechenweg Was genau meinst du? Ich hatte den Beitrag nach dem Absenden noch mal abgeändert, da ich mich verrechnet hatte - daher meine Anmerkung. >obersten Grenzfrequenz der DFT --> wie kommste drauf ? weil i = N/2 (weiter oben schonmal angesprochen) Nach dem Abtasttheorem (Shannon-Nyquist) kannst du bei einer Sample-Rate von f nur Frequenzen kleiner f/2 auslesen. Alle Frequenzen, die darüber liegen liefern die sog. Spiegelfrequenzen. f/2 entspricht immer einer Wellenlänge von 2 Samples. Stell dir einfach mal ein Sinussignal mit einer Wellenlänge von 2 Samples vor. Die Werte sähen wie folgt aus:
1 | s(x) = sin((x / 2) * 2 * Pi) = |
2 | sin(x * Pi) = <0, 0, 0, 0, ...> |
Je nach Phasenverschiebung kämen andere konstante Folgen raus. Nicht wirklich toll um in die Frequenzdomäne transformiert zu werden, oder? ;) Deshalb meinte ich weiter oben ja auch >oben ist [0 .. N - 1] angegeben... kommt wir aber falsch vor. >Verwendet du die KISS FFT ? Kannst du das Beispiel von mir mal bei Dir >laufen lassen ? Dabei kannst Du ja auch dein Ergebnis überprüfen =) Ich kann noch nicht mal wirklich in C programmieren und die DFTs schreibe ich meist kurz selber - dann weiß ich, dass ich's verstanden habe und ich kann mir sicher sein, dass ich alle Fehler in meinem Code finden werde ;) Hoffe, das hilft dir was weiter.
War gestern wohl doch was zu spät...
>weil i = N/2 (weiter oben schonmal angesprochen)
Das stimmt natürlich nicht, da du N = 4 und i = 1 gewählt hast. Du
kommst also nicht an die besagte Grenze.
Gruß
Kai
PCM ist eine digitale Modulation, berechnet man aus diesem das Spektruma, ist dies das Spektrum des Übertragungskanal. Zuerst muss aus PCM die normalen Werte gewonnen werden.
@gast Keine Ahnung, was du uns damit sagen willst... vor allem nicht, was das mit dem Thread zu tun haben soll.
> Werde aber mit float weiter testen.
So, mit
1 | float
|
funktionierts ! ABER : Was habe ich denn jetzt als Ergebnis ? Ich bekomme für F(i) jeweils ne komplexe Zahl. Aus dieser kann ich den Betrag und den Winkel berechen. Und jetzt ? Wie bekomme ich den Bezug von i zur jeweiligen Frequenz ? DENN : Wenn ich den Betrag über i zeichne bekomme ich doch: Betrag . . | | . | | | | . --|--|--|--|-...-|---> i 1 2 3 4 N/2 Das Sprachsignal wird mit Ta abgetastet. D.h. ich bekomme jeweils zu k*Ta einen 8-Bit Wert. Bitte : Erkläre mir was ich mit diesem Ergebnis machen soll. Thx, eric
> Bitte : Erkläre mir was ich mit diesem Ergebnis machen soll.
Also ich meine wie bekomme ich den Bezug zur Frequenz ?
Hi Hägar, danke für deine Antwort. Sorry, aber was ist 'ni' ? Und geh ich richtig in der Annahme, dass ni/N < 1 (eigentlich < 0.5) sein muss ? Denn, fs ist > 2*fm (fm = maximale abzutastende Frequenz). Es sollten ja nur Frequenzen von 100Hz bis 4kHz von Interesse sein. Konkret: F(1) = -3 +5j => |F(1)| = 5.83 N = 6 fs = 8kHz ni = ? Thx und Gruß, eric
Also, hab jetzt mal etwas den Kopf zerbrochen =| Diskretes Sprachsignal: u(k*Ta) Spektrum Sprachsignal: F(i) Frequenz für F(i): (i*SAMPLE_RATE)/NUMBER_OF_SAMPLES Korrekt ? eric
ni ist der Wert der FFT an der Stelle i. So wie du das im letzten Post geschrieben hast. Das ist auch relativ klar: Du berechnest N Werte F(i). d.h. dein Spektrum bewegt sich zwischen 0 und N-1. F(0) ist dann die "DC-Komponente" wenn man mal Aliasing ausser Acht lässt. Mit jedem F(i+1) erhöht sich die Frequenz um einen Bruchteil von fs (Sample-Frequenz) und zwar genau um 1/N. In absoluten Frequenzen ausgedrückt also 1/N*fs.
> So wie du das im letzten Post geschrieben hast. freq(i) = (i*SAMPLE_RATE)/NUMBER_OF_SAMPLES //i = [0....N/2] > ni ist der Wert der FFT an der Stelle i meinst du damit den Betrag von F(i) ? Denn das ist nicht dasselbe wie i. Sorry, aber was verstehst Du unter >der Wert der FFT an der Stelle i Thx, eric
Wie kann ich die Effektivwerte der einzelnen Frequenzen in dB umrechnen ? Also damit die y-Achse wie im Dateianhang aussieht ? Denn ich habe ja keine Bezugsgröße wie ein Eingangssignal. Ich habe nur das Ausgangssignal selber. 10 * LOG(P/Po) // P bekannt (PCM), Po = ? hat jemand ne Idee ? Thx, eric
Bitte 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.