Forum: Digitale Signalverarbeitung / DSP / Machine Learning FFT auf PCM Signal anwenden


von Daniel K. (Firma: southpark SOP) (eric2000)


Lesenswert?

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 ?

von Hägar (Gast)


Lesenswert?

Was willst du mit dem Spektrum des PCM-Signals? Oder doch das Spektrum 
des Sprachsignals?

von Daniel K. (Firma: southpark SOP) (eric2000)


Lesenswert?

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.

von Karl (Gast)


Lesenswert?

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.

von Daniel K. (Firma: southpark SOP) (eric2000)


Lesenswert?

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

von Kai G. (runtimeterror)


Lesenswert?

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?

von Hägar (Gast)


Lesenswert?

> 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.

von Karl (Gast)


Lesenswert?

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.

von Daniel K. (Firma: southpark SOP) (eric2000)


Lesenswert?

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.

von Karl (Gast)


Lesenswert?

Ach ja: Wenn du die FFTW einsetzt, musst du dein Programm entweder unter 
der GPL veröffentlichen oder dafür zahlen!

von Karl (Gast)


Lesenswert?

>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)

von Kai G. (runtimeterror)


Lesenswert?

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.

von Karl (Gast)


Lesenswert?

>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.

von Daniel K. (Firma: southpark SOP) (eric2000)


Lesenswert?

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 !

von Kai G. (runtimeterror)


Lesenswert?

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)

von Hägar (Gast)


Lesenswert?

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.

von Daniel K. (Firma: southpark SOP) (eric2000)


Lesenswert?

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

=)

von Andreas S. (andreas) (Admin) Benutzerseite


Lesenswert?

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

von Daniel K. (Firma: southpark SOP) (eric2000)


Lesenswert?

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

von Daniel K. (Firma: southpark SOP) (eric2000)


Lesenswert?

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

von Kai G. (runtimeterror)


Lesenswert?

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.

von gast (Gast)


Lesenswert?

Verstehe ich das richtig, dass du das PCM-Signal abtastet und erwartest, 
das dies dem Spektrum des Audiosignal entspricht?

von Kai G. (runtimeterror)


Lesenswert?

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.

von Daniel K. (Firma: southpark SOP) (eric2000)


Lesenswert?

> 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

von Kai G. (runtimeterror)


Lesenswert?

>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.

von Kai G. (runtimeterror)


Lesenswert?

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

von gast (Gast)


Lesenswert?

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.

von Kai G. (runtimeterror)


Lesenswert?

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

von eric (Gast)


Angehängte Dateien:

Lesenswert?

> 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

von eric (Gast)


Lesenswert?

> Bitte : Erkläre mir was ich mit diesem Ergebnis machen soll.

Also ich meine wie bekomme ich den Bezug zur Frequenz ?

von Hägar (Gast)


Lesenswert?

fi = ni/N * fs

von eric (Gast)


Lesenswert?

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

von eric (Gast)


Angehängte Dateien:

Lesenswert?

Suche Ergebnis ähnlich wie im Dateianhang ! =)

von eric (Gast)


Lesenswert?

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

von Hägar (Gast)


Lesenswert?

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.

von eric (Gast)


Lesenswert?

> 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

von Hägar (Gast)


Lesenswert?

Sorry. War blöd ausgedrückt. ni = i

von eric (Gast)


Angehängte Dateien:

Lesenswert?

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
Noch kein Account? Hier anmelden.