www.mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Power Spectrum und FFT


Autor: Spielername (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo

Ich interessiere mich für die Signalverarbeitung (bin ganz am Anfang 
eines E-Technik Studiums) und wollte mich da vor den eigentlichen 
Vorlesungen schon mit beschäftigen.

Mit Hilfe von Visual Studio und C möchte ich einen Rauschfilter 
programmieren.

Ich hab schon einiges programmiert, stecke jetzt aber etwas bei der 
Rücktransformation fest.

Nach dem Einlesen einer *.wav Datei wird eine FFT ausgeführt.
Aus dem Real- und Imaginärteil bilde ich den Absolutwert.

Es wird eine 1024-Punkte FFT berechnet.
Ich gehe davon aus, dass der erste Frame (also 1024 Werte) reines 
Rauschen sind und bestimme die Rauschleistung des Frames.

Diese subtrahiere ich dann von allen nachfolgenden Frames.
Also erstmal ein sehr einfaches Filter.

Nun habe ich aber das Power Spectrum. Und wenn ich dies zur 
Rücktransformation auf die FFT gebe hört sich das einfach nur grauselig 
an.
Das eigentliche Audiosignal versteht man aber trotzdem.

Das Problem tritt auch auf, wenn ich das Filtern nicht durchführe.
Liegt also am Power Spectrum.

Wie bekomme ich daraus jetzt wieder ein Audiosignal, was vernünftig 
klingt?

Hoffe auf etwas Hilfe.

Gruß

: Verschoben durch Admin
Autor: Fins Bury (finsbury)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
für die Rücktransformation meinst du hoffentlich  IFFT.  nicht die Phase 
wegwerfen. Wenn du nur den Betrag in die IFFT gibst, kanns nix werden

Autor: Spielername (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo und Danke für die Antwort.

Richtig, ich meine die IFFT.

Die berechneten Werte der Phase nach der FFT gebe ich auch wieder mit 
auf die IFFT.

Trotzdem hört es sich leider echt schlimm an.

Gruß

Autor: Spielername (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hat nicht noch jemand eine Idee?

Ich habe die FFT dieser Seite verwendet:
http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/dft/

Habe sie auch als txt im Anhang hinterlegt.

Oder kennt jemand vielleicht eine ähnlich einfach zu handhabende FFT in 
C, mit der das funktioniert?

Gruß

Autor: Xeraniad X. (xeraniad)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Den gezeigten Algorithmus habe ich nicht geprüft, sollte aber OK sein.
Hier ein einige Überlegungen zur Fehlersuche.
° Tritt das Problem auch dann auf, wenn zwischen Lesen und Schreiben der 
.wav-Datei keine FFT durchgeführt wird (d. h. nur .wav-Werte in Vektor 
schreiben und wieder zurück)?
° Gibt es ein Problem bei der Konvertierung der Ganzzahlen aus der 
.wav-Datei in die double-Realteile des FFT-Input und/oder bei der 
Konvertierung der Realteile des iFFT-output in .wav-Integers (Rundung)?
° Wurden die double-Imaginärteile für den FFT-Input richtig auf 0 
gesetzt?
° Tritt das Problem auch dann auf, wenn zwischen FFT und iFFT keine 
Umformungen komplex→polar, polar→komplex durchgeführt werden, d. h. der 
FFT-Output (komplexes Spektrum) gleich so wieder bei iFFT eingegeben 
wird?
° Vermutlich wurden Amplitude r[k] und Phase φ[k] aus den komplexen 
Koeffizienten c[k] ähnlich wie folgt berechnet r[k] = 2·|c[k]|, φ[k] = 
arg(c[k]). Wurde φ[k] richtig mit atan2 (Im(c[k]), Re(c[k])), nicht mit 
atan (Im(c[k]) /Re(c[k])) ausgewertet?
° Wurde die Übersetzung polar→komplex richtig wie folgt durchgeführt 
c[k] = ½·r[k] ·exp(j·φ[k])?
° Wurde auch die rechte Hälfte c[n-k] = conj(c[k]) des komplexen 
Spektrums richtig an die iFFT zurückgegeben?

Autor: Xeraniad X. (xeraniad)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Bei diesem Thread [[Beitrag "Verständnis FFT-Übung"]] habe 
ich eine rekursive Variante angegeben, auch getestet; passt etwa auf 
eine Seite.

Autor: Spielername (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,

danke für die Antwort.

Xeraniad X. schrieb:
> ° Tritt das Problem auch dann auf, wenn zwischen Lesen und Schreiben der
> .wav-Datei keine FFT durchgeführt wird (d. h. nur .wav-Werte in Vektor
> schreiben und wieder zurück)?
Nein, kein Problem.

Xeraniad X. schrieb:
> ° Gibt es ein Problem bei der Konvertierung der Ganzzahlen aus der
> .wav-Datei in die double-Realteile des FFT-Input und/oder bei der
> Konvertierung der Realteile des iFFT-output in .wav-Integers (Rundung)?
Ich denke nicht, da .. (übernächster Punkt)

Xeraniad X. schrieb:
> ° Wurden die double-Imaginärteile für den FFT-Input richtig auf 0
> gesetzt?
Jop, schon als Ausgabe in eine txt getestet.
Der erste Absolutwert berechnet sich auch aus Realteil + 0i = Realteil

Xeraniad X. schrieb:
> ° Tritt das Problem auch dann auf, wenn zwischen FFT und iFFT keine
> Umformungen komplex→polar, polar→komplex durchgeführt werden, d. h. der
> FFT-Output (komplexes Spektrum) gleich so wieder bei iFFT eingegeben
> wird?
.. hier kein Problem auftritt.
Einlesen -> FFT -> IFFT -> Ausgabe erzeugt ein sauberes Audiosignal

Xeraniad X. schrieb:
> ° Vermutlich wurden Amplitude r[k] und Phase φ[k] aus den komplexen
> Koeffizienten c[k] ähnlich wie folgt berechnet r[k] = 2·|c[k]|, φ[k] =
> arg(c[k]). Wurde φ[k] richtig mit atan2 (Im(c[k]), Re(c[k])), nicht mit
> atan (Im(c[k]) /Re(c[k])) ausgewertet?
> ° Wurde die Übersetzung polar→komplex richtig wie folgt durchgeführt
> c[k] = ½·r[k] ·exp(j·φ[k])?
Hier kann ich dir gerade leider nicht folgen.
Wie gesagt, die Vorlesung zu dem Thema fehlt mir leider noch ganz. 
Startet erst in 3 Semestern :-)
Bei der Übersetzung polar->komplex, meinst du die Stelle, bevor ich 
wieder die IFFT ausführe?

Xeraniad X. schrieb:
> ° Wurde auch die rechte Hälfte c[n-k] = conj(c[k]) des komplexen
> Spektrums richtig an die iFFT zurückgegeben?
Ich arbeite im Moment noch mit dem ganzen Spektrum.
Optimierung kommt später :-)

Xeraniad X. schrieb:
> Bei diesem Thread [[Beitrag "Verständnis FFT-Übung"]] habe
> ich eine rekursive Variante angegeben, auch getestet; passt etwa auf
> eine Seite.
Danke, werde ich mir einmal ansehen.

Autor: Xeraniad X. (xeraniad)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
> Einlesen -> FFT -> IFFT -> Ausgabe erzeugt ein sauberes Audiosignal.
Klingt schon mal gut, das Problem dürfte damit eingegrenzt sein. Dies 
spricht auch für die Korrektheit des verwendeten Algorithmus.
Dann geht es vermutlich um die Wandlungen kartesisch→polar und wieder 
zurück (falls meine Vermutung zutrifft, dass so etwas stattfindet).
> Bei der Übersetzung polar->komplex, meinst du die Stelle, bevor ich wieder > die 
IFFT ausführe?
Ja, genau die Stelle. Der von Dir gepostete Algorithmus verwendet die 
beiden double-Arrays x[0..n-1] und y[0..n-1]. Nach der FFT liegen die 
komplexen Koeffizienten des Spektrums in kartesischer Form vor. Für eine 
polare Darstellung holst Du vermutlich aus den (kartesisch) komplexen 
Koeffizienten c[k] (= x[k]+j·y[k]) Amplituden-Werte r[0] = x[k], r[k] = 
2·√({x[k]}²+{y[k]}²) und Phasen-Werte φ[0] = 0, φ[k] = atan2 (y[k], 
x[k]), ... Übrigens sind vom komplexen Spektrum nur die Werte x[0] 
(arithmetischer Mittelwert, Gleichanteil), x[1..½·n], y[1..½·n-1] 
wichtig, die anderen sind Null bzw. redundant (y[0] = 0, y[½·n] = 0, 
x[n-k] = x[k], y[n-k] = -y[k]). Es ist nicht üblich (aber auch nicht 
verboten), Spektral-Linien-Information für Frequenz-Indices k > ½·n 
polar darzustellen (½·n ist der Frequenz-Index für die 
Nyquist-Frequenz).
(Möglicherweise verzichtest Du bei r[k] auf den Faktor 2, wichtig ist 
jedoch, dies bei der Rückumformung nach kartesisch zu berücksichtigen.)
Die Umformung polar→kartesisch muss die Wandlung zuvor rückgängig 
machen, d. h. x[0] = r[0], y[0] = 0, x[k] = ½·r[k]·cos(φ[k]), y[k] = 
½·r[k]·sin(φ[k]).

Falls Du Deine Umformungen posten möchtest, finden wir das Problem 
möglicherweise. (z. B. sicher, dass da nicht unterwegs mit floating 
point geringerer Genauigkeit gerechnet wird?).

Autor: Spielername (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Wenn ich jetzt richtig liege, habe ich gleich zwei Fehler bei mir 
gefunden.

Ich habe den Code teilweise mal als txt angehängt.
Hatte zum Test die Subtraktion schon rausgenommen und nur meine Funktion 
für den Betrag zwischen FFT und IFFT aufgerufen.
Auch da treten die Probleme auf.

Wenn ich das jetzt richtig verstanden habe, dann habe ich den Faktor 2 
vergessen.
Und was noch schlimmer ist: Die Wandlung kartesisch -> polar des 
Imaginärteils vollkommen vergessen.

Autor: Spielername (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Neuere code.txt, hatte beim Kopieren einige Variablennamen übersehen

Autor: Spielername (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Und einen alten Funktionsheader! Tut mir leid. Da habe ich noch mit 
Variablen vom Typ float gearbeitet. Mittlerweile ist es überall double.

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

Bewertung
0 lesenswert
nicht lesenswert
Ich habe mir den den Code nicht angesehen, aber es gibt ein paar 
grundlegende Probleme mit deinem Ansatz. Du kannst das Signal nicht 
einfach in Blöcke aufteilen, jeden einzelnen transformieren, bearbeiten 
und wieder zurücktransformieren. Außerdem kannst du das Rauschspektrum 
nicht schätzen indem du nur von einem Block eine FFT machst, dazu musst 
du über mehrere Blöcke mitteln.

Das im Detail hier im Forum zu erklären würde dir wahrscheinlich nicht 
viel bringen. Wenn es dir ernst ist, besorg dir ein paar Bücher und 
investier ein bisschen Zeit in die Grundlagen (Filter, DFT, 
Spektralschätzung); Stichwörter für die Verarbeitung von Signalen im 
Frequenzbereich sind dann "FFT-Filterbank", für die Rauschreduktion 
"spektrale Subtraktion". Zum Implementieren und Ausprobieren würde ich 
MATLAB (gibt's an der Uni i.d.R. kostenlos) oder Octave (kompatibel, 
kostenlos) statt C empfehlen.

Einfacher wäre es auf jeden Fall bis zur entsprechenden Vorlesung zu 
warten. Vielleicht wird sogar von der Uni eine Übung oder ein Praktikum 
angeboten in dem man genau solche Dinge macht.

Autor: Spieler N. (spielername)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
So, jetzt habe ich es denke ich.

Die Umformung polar -> kartesisch habe ich erst gar nicht 
berücksichtigt.
Ich habe ein in MATLAB erstelltes Modell dazu.
Da nutze ich abs und angle für den Real und Imaginärteil.
Gebe diese Werte dann aber so auf die IFFT. Und da klappt es.
Aus dem Grund habe ich da nicht weiter drüber nachgedacht.

Die Umformung vor der IFFT habe ich jetzt so realisiert, wie in der txt 
zu sehen ist.

Nur sind damit im Ergebnis nach der IFFT alle Werte Null.

edit: @ andreas:
Lesen muss ich zu dem Thema sicher noch einiges, das stimmt.
Im Moment geht es mir auch erst mal nur um die FFT und das Power 
Spectrum.
Wenn das endlich funktioniert, möchte ich mich mehr mit dem Thema 
Rauschreduktion beschäftigen.

EDIT: ES klappt!
Ich hab in meiner Verwirrtheit im Code 0,5 statt 0.5 stehen gehabt.

Vielen Dank für die Hilfe.
Jetzt kann ich mich mehr mit dem Thema Rauschreduktion auseinandersetzen 
:-)

Autor: Xeraniad X. (xeraniad)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Guten Tag

Zu den Umformungen (kart←→polar) habe ich noch ein paar Anmerkungen.

Wenn Du nicht die echten Amplituden der Teilschwingungen, sondern die 
Beträge der komplexen Koeffizienten darstellen willst, kannst Du die 
Funktion "abscplxarr" so verwenden. Dann sollte aber konsequenterweise 
auch kein Faktor ½ bei den Funktionen rueck* verwendet werden.

Gemäss meiner Ansicht könnten die beiden Funktionspaare zusammengefasst 
werden. Ich habe auch etwas die Namen der Argumente geändert, aber sonst 
den Stil beibehalten.

void cart2pol (double const * re, double const * im, double * abs, 
double * arg)
{ for (int k = 0;  SAMPLE > k;  ++ k) {
    abs [k] = sqrt (re[k]*re[k]+im[k]*im[k]);
    arg [k] = atan2 (im[k], re[k]);
  }  # k
  return;
}

void pol2cart (double const * abs, double const * arg, double * re, 
double * im)
{ for (int k = 0;  SAMPLE > k;  ++ k) {
    re [k] = abs [i] * cos (arg [k]);
    im [k] = abs [i] * sin (arg [k]);
  }  # k
  return;
}

Die polare Form gibt, so wie sie hier angegeben wird, genau den Betrag 
der komplexen Koeffizienten wieder. Für die Darstellung im Histogramm 
und um Aussagen bezüglich Leistung zu ermöglichen, sollten diese 
Absolutwerte gedeutet /angepasst werden.

k = 0: abs [0] = Gleichanteil
0 < k < ½·SAMPLE: abs [k] = halbe Amplitude des Frequenz-Anteils k
k = ½·SAMPLE: Amplitude des Nyquist-Frequenzanteils
½·SAMPLE < k < SAMPLE: halbe Amplitude des Frequenz-Anteils SAMPLE-k

Um die echten Amplituden für [1..½·SAMPLE-1] zu erhalten, kann der Teil 
[1+½·SAMPLE..SAMPLE-1] rechts nach links geklappt werden. Aber das wird 
dann in den Vorlesungen behandelt.

Weitere mögliche Entwicklungen:
° Ausnützen der erwähnten Symmetrie c[n] = conj(c[n-k]) etc. ermöglicht 
Halbierung des Speicherplatzbedarfes für die polare Darstellung.
° Es könnten die Paare {re, im} und {abs, arg} in structs 
zusammengefasst und Arrays davon als Argumente übergeben werden 
(reduzierte Anzahl Argumente /Übersicht?).

schönen Tag

Autor: Xeraniad X. (xeraniad)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Korrektur: ... Symmetrie c[ k ] = conj(c[n-k]) ...
Zu den Leistungen.
k = 0: Amplitude: r[0] = c[0] = |c[0]|, Leistungsmass: p[0] = r[0]² = 
|c[0]|² (weil DC).
k = 1..½·n-1:  Amplitude: r[k] = 2·|c[k]|, Effektivwert: r[k]eff = 
½·√2·r[k] = √2·|c[k]| (weil DC), Leistungsmass: p[k] = r[k]eff² = 
2·|c[k]|².
k = ½·n: Amplitude: r[½·n] = c[½·n] = |c[½·n]|, Effektivwert: r[½·n]eff 
= ½·√2·r[½·n] (Annahme: DC), Leistungsmass: p[½·n] = r[½·n]eff² = 
½·r[½·n]² = ½·|c[½·n]|². Hier kann genauso Rechteck angenommen werden, 
dann stimmt auch die Leistungsbilanz.

Autor: Xeraniad X. (xeraniad)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Schon wieder Korrektur. Habe zuvor 2 mal "DC" anstelle von "AC" 
geschrieben.

Autor: Wolfgang M. (womai)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
> Einfacher wäre es auf jeden Fall bis zur entsprechenden Vorlesung zu
> warten. Vielleicht wird sogar von der Uni eine Übung oder ein Praktikum
> angeboten in dem man genau solche Dinge macht.

Warum warten? Zumindest mich persoenlich freut es, wenn jemand so viel 
Initiative hat, sich aus Interesse selber in ein Thema einzuarbeiten, 
anstatt zu warten, bis er/sie es vorgekaut bekommt! Durch Fehler lernt 
man ohnehin mehr, als wenn alles gleich auf Anhieb funktioniert. Und 
spaeter im Beruf kriegt man die Aufgabenstellung (und. ev. Loesung) 
meist eben auch nicht sauber aufbereitet praesentiert...

Viel Erfolg, "Spielername" (melde Dich doch mal an!)

Wolfgang

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

Bewertung
0 lesenswert
nicht lesenswert
Wolfgang M. schrieb:
> Durch Fehler lernt man ohnehin mehr, als wenn alles gleich auf Anhieb
> funktioniert

Durch Fehler die man nicht verstehen kann weil die Grundlagen fehlen 
lernt man nichts. Und dafür dass alles auf Anhieb funktioniert besteht 
auch nach der Vorlesung keine Gefahr; aber man kann verstehen, warum 
etwas funktioniert oder nicht funktioniert. Ansonsten kann man nur durch 
exaktes Verfolgen von Anleitungen etwas erreichen, oder sehr begrenzt 
durch Trial and Error.

> spaeter im Beruf kriegt man die Aufgabenstellung (und. ev. Loesung)
> meist eben auch nicht sauber aufbereitet praesentiert...

Wenn man im Beruf z.B. ohne Wissen über Antennen die Aufgabe bekommt 
eine Antenne zu bauen, dann wird man das Problem aber nicht so angehen 
dass man nach Gutdünken irgendwelche Stäbe zusammenlötet. Und selbst 
wenn das zufällig funktioniert hat man dabei nichts gelernt.

Autor: Wolfgang M. (womai)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Andreas Schwarz schrieb:
> Durch Fehler die man nicht verstehen kann weil die Grundlagen fehlen
> lernt man nichts.

Naja, oder man merkt eben, dass man es noch nicht verstanden hat, und 
graaebt dann aus diesem Grund tiefer. Der Fragesteller macht auf mich 
den Eindruck, dass er ohnehin einigermassen systematisch den Fehler 
einzugrenzen versucht - von blinden Trial-and-Error sehe ich da wenig. 
Fremden Code zu verwenden, ohne dessen Arbeitsweise bis ins Detail zu 
verstehen ist ja keine Schande - sonst duerfte auch niemand z.B. 
Grafik-Rendering-DLLs einsetzen - da verstehen wohl auch nur die 
wenigsten Anwender, was da ganz genau hinter der Buehne ablaeuft.

Es gibt eben unterschiedliche Zugaenge zu einem Problem - manche Leute 
studieren alles ganz genau bevor sie den ersten vorsicvhtigen Schritt 
machen, andere probieren gleich mal was aus und sehen dann weiter. Ich 
denke, keine Variante ist da allgemein richtig oder falsch.

Meine Bemerkung war absolut nicht boese gemeint, mich freut einfach, 
Interesse am Studium zu sehen, das ueber die Pflichterfuellung 
hinausgeht. Sowas werden dann die richtigen Ingenieure!

Autor: Spieler N. (spielername)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo

Also das mir noch eine Menge theoretischer Grundlagen fehlen steht außer 
Frage.
Aber das kommt ^^ Spätestens mit der Vorlesung ;-)

Ich finde Signalverarbeitung einfach nur sehr interessant.

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.