FFT auf PCM Signal anwenden

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.

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?

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.

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.

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)

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(++) ??

=)

Alternative zu FFTW (mit freierer Lizenz):
http://sourceforge.net/projects/kissfft/

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

Jungs, ihr habt doch sicher nen Tip, gell. =)

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.

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

@gast
Keine Ahnung, was du uns damit sagen willst... vor allem nicht, was das 
mit dem Thread zu tun haben soll.