mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP FFT Analyse Phasenlage und Gruppenlaufzeit berechnen


Announcement: there is an English version of this forum on EmbDev.net. Posts you create there will be displayed on Mikrocontroller.net and EmbDev.net.
Autor: Zombie (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Forum,

ich würde gerne etwas mit digitaler Signalverarbeitung 
herumexperimentieren, leider happert es bereits bei der Basis. Mag 
jemand helfen?

Mein Idee:
ich erzeuge mittels iFFT ein Testsignal mit bekannter Amplitude (1.0) 
und Phasenlage (0.0) pro Frequenz (oder als Komplexe-Zahl: real=1.0; 
imag=0.0). Dieses Testsignal verarbeite ich mit meiner zu testenden 
Funktion. Anschliessend führe ich eine normale FFT durch, und zeige das 
Ergebnis an.

Ich habe das mal als Beispiel angehängt. Einmal nur den Code und einmal 
das ganze Juce-Projekt.

Für meinen ersten Versuch verwende ich ein 1-pol Tiefpassfilter mit 
einer Grenzfrequenz von fc = (fs * 0.03125) bei einer virtuellen 
Samplingfrequenz von fs = 44'100Hz.

Das Ergebnis in grober Auflösung (Frequenz, Amplitude, Phase)
0,-22.8513,0
689.063,-0.361587,-23.8445
1378.13,-2.39234,-39.5591
2067.19,-4.48345,-48.1488
2756.25,-6.32989,-52.5538
3445.31,-7.91214,-54.5978
4134.38,-9.27027,-55.2451
4823.44,-10.4469,-55.0161
5512.5,-11.4765,-54.2072
6201.56,-12.3857,-52.9975
6890.63,-13.1947,-51.4999
7579.69,-13.9191,-49.7896
8268.75,-14.5708,-47.9178
8957.81,-15.1593,-45.9208
9646.88,-15.6923,-43.825
10335.9,-16.1758,-41.65
11025,-16.6149,-39.4108
11714.1,-17.0136,-37.1189
12403.1,-17.3754,-34.7833
13092.2,-17.7032,-32.4113
13781.3,-17.9993,-30.0088
14470.3,-18.2657,-27.5806
15159.4,-18.5043,-25.1307
15848.4,-18.7164,-22.6625
16537.5,-18.9032,-20.1788
17226.6,-19.0659,-17.6821
17915.6,-19.2053,-15.1746
18604.7,-19.3221,-12.6583
19293.8,-19.4169,-10.1348
19982.8,-19.4902,-7.60583
20671.9,-19.5423,-5.07277
21360.9,-19.5734,-2.53705

Die Amplitude könnte stimmen, aber die Phasenlage meiner Meinung nach 
nicht. Diese müsste sich doch zu höheren Frequenzen immer in die selbe 
Richtung bewegen und nicht plötzlich umdrehen und nach 0. Was mache ich 
falsch?

Und als nächstes möchte ich noch die Gruppenlauf brechnen. Aber das 
Thema macht eigentlich erst sinn, wenn die Phasenberechnung 
funktioniert, darum werde ich später darauf zurück kommen.

Vielen Dank und grüsse

Autor: Dergute W. (derguteweka)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Moin,

Doch das kann schon so hinhauen. So sollte der Phasen- und 
Amplitudenverlauf ausschauen.

Gruss
WK

Autor: Zombie (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Weka,

Danke für das Validieren!

Aber ich verstehe es trotzdem nicht. Ich erwartete eigentlich so einen 
Verlauf 
https://de.wikipedia.org/wiki/Tiefpass#/media/Datei:Bodediagramm_Tiefpass.svg

Ist mein Tiefpassfilter-Code falsch?

grüsse

Autor: Dergute W. (derguteweka)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Moin,

Das von dir verlinkte Bild passt zu einem analogen Tiefpass; wenn du das 
mit einem digitalen System nachbaust, dann gibt's ein paar 
"Zusatzeffekte" durch Aliasspektren, etc.
Die Frequenz ist in deinem Bild auch logarithmisch aufgetragen und nicht 
linear. Es koennte auch "noch aehnlicher" werden, wenn du das 
Ausgangssignal noch um 1 Takt verzoegerst.

Zombie schrieb:
> Ist mein Tiefpassfilter-Code falsch?
Hm - keine Ahnung, ich fuercht' ich kann zu wenig c++, um da was 
rauszulesen.
Aber ich glaub' eher, dass deine Erwartungen "falsch" sind :-)

Gruss
WK

Autor: Zombie (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Dergute W. schrieb:
> Die Frequenz ist in deinem Bild auch logarithmisch aufgetragen und nicht
> linear.

Ja, aber es geht mir jetzt nur um den Phasenverlauf, der mich irritiert.

> Es koennte auch "noch aehnlicher" werden, wenn du das
> Ausgangssignal noch um 1 Takt verzoegerst.

Ja, von der Form her schon, aber dann ist die Phase -180° bei fs/2. 
Vermutlich wäre eine Verzögerung um ein halbes Sample richtig ;)

> Hm - keine Ahnung, ich fuercht' ich kann zu wenig c++, um da was
> rauszulesen.

im Prinzip ist es
coef = exp(-2 * PI * fc / fs)

und dann für jedes Sample
out[n] = out[n - 1] * coef + in * (1 - coef)

Dergute W. schrieb:
> Aber ich glaub' eher, dass deine Erwartungen "falsch" sind :-)

Sieht fast so aus ;) Ich habe vorallem erwartet, dass die Phasenlage bei 
fs/2 -90° oder so beträgt, wie bei einem analogen LPF und nicht 0°.

Ich habe das Filter von https://www.dspguide.com/ch19/2.htm

Er schreibt dort auch:
> They are easy to program, fast to execute, and produce few surprises.
Vielleicht ist der Phasenverlauf eines der "few surprises".

Wie wäre denn der Code um ein simples, analoges RC-Glied 1:1 zu 
simulieren (in der Zeitebene)?

Aber es ist schon mal gut zu wissen, dass der Synthese-/Analysepart 
funktioniert.

grüsse

Autor: Dergute W. (derguteweka)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Moin,

Von einem analogen auf ein digitales Filter zu kommen, was in Bezug auf 
Frequenz-,Phasengang, Impulsantwort identisch ist, haut normalerweise 
nicht hin.
Man kann, wenn man sich auf bestimmte Eigenschaften beschraenkt, mehr 
oder weniger gute Naeherungen finden, z.B. mit der bilinearen 
Transformation. Oder Annaeheren der Impulsantwort. Aber das haut nie 
100% hin.

Im Analogen ist halt der olle RC-Tiefpass eines der simpelsten Filter 
und das Dingens, was du gebaut hast, eine Art der Naehrung im Digitalen. 
Haut fuer die Impulsantwort gut hin, fuer den Frequenzgang eher nur bei 
niedrigen Grenzfrequenzen und - wie du gemerkt hast - fuer die Phase nur 
sehr wenig.

Gruss
WK

Autor: Helmut S. (helmuts)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
> ich erzeuge mittels iFFT ein Testsignal mit bekannter Amplitude (1.0)

Wozu bneötigst du die iFFT? Das ist doch schon mal die erste Möglichkeit 
zu scheitern.

Definiere einfach ein Signal cos(w*T) und fülle damit den Datenvektor 
über einen Zeitraum von 110*T (110 Perioden). Verarbeite diese 
Datenvektor mit deinem Tiefpassfilter. Wirf die Daten der ersten 10*T 
Sekunden weg da du den eingeschwungenen Zustand analysieren willst.
Auf diesen verkürzten Datenvektor wendest du dann die DFT bzw die FFT 
an. Damit solltest du jetzt nur bei einer Frequenz einen nennenswerten 
Beitrag haben mit -45°. Alle anderen Beiträge im Spektrum sind im 
Bereich 10^-6? (numerisches Rauschen). Das wird nur deshalb so ideal, 
weil bei dir auch die Abtastfrequenz exakt ein Vielfaches der 
Signalfrequenz ist.

: Bearbeitet durch User
Autor: Zombie (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Danke WK für deine hilfreichen Beiträge. Das Thema Phasenlage ist somit 
für mich (im Moment) hinreichend geklärt.

Kommen wir zum zweiten Bereich: die Berechnung der Gruppenlaufzeit.
Bei Wikipedia steht:
> Die Gruppenlaufzeit ergibt sich aus dem Phasengang durch:
> 
https://wikimedia.org/api/rest_v1/media/math/render/svg/9acaeeeeaae8a09a597469f611b83dc721cffaef 
(<- klick)

Was ist d? Und wenn d sich als Faktor auf beiden Seiten der Division 
befindet, kann man es auch einfach wegkürzen, nicht? Dann bliebe:
groupdelay = -1.0 * (phase / kreisfrequenz)

Das kann doch aber nicht stimmen... das groupdelay kann ja auch mehr als 
einen Zyklus gross sein.

Helmut S. schrieb:
>> ich erzeuge mittels iFFT ein Testsignal mit bekannter Amplitude (1.0)
>
> Wozu bneötigst du die iFFT? Das ist doch schon mal die erste Möglichkeit
> zu scheitern.

Warum scheitern? Mir ist schon klar, dass ich das Testsignal auch anders 
und wesentlich ressoursenschonender generieren könnte, aber für mich war 
es für den Anfang die einfachste und sicherst Möglichkeit ein 
definiertes, synchrones, breitbandiges Signal zu generieren. Die iFFT 
ist schliesslich eine getestete Funktion (die FFT Klasse gehört zum Juce 
Framework).

grüsse

Autor: Dergute W. (derguteweka)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Moin,

Zombie schrieb:
> Was ist d? Und wenn d sich als Faktor auf beiden Seiten der Division
> befindet, kann man es auch einfach wegkürzen, nicht? Dann bliebe:

Nee, das d kannst du natuerlich nicht weggkuerzen, der Ausdruck bedeutet 
die Ableitung der Funktion phi(omega) nach omega; also die Steigung.

Gruss
WK

Autor: Zombie (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Dachte mir schon fast, dass d keine Variable oder Konstante ist, 
trotzdem war das etwas verwirrend.
Mir ist in der Vergangenheit auch schon aufgefallen, dass die 
Gruppenlaufzeit jeweils dort am höchsten ist, wo die Phase sich am 
stärksten ändert.

Ich hab' das mal so implementiert:

Helfer Funktion:
double phase_diff(double a, double b) {
  double c = std::fmod(a - b + 180.0, 360.0);
  if (c < 0.0) {
    c += 360.0;
  }
  return c - 180.0;
}
und dann folgender Code für die Umrechnung
  std::ofstream file("fft.txt");
  double frequency = 0.0;
  double prev_frequency = 0.0;
  double amplitude = 0.0;
  double phase = 0.0;
  double prev_phase = 0.0;
  double groupdelay = 0.0;
  for (size_t i = 0; i < fft.getSize(); i += 2) {
    std::complex<float> c(arr[i], arr[i + 1]);
    frequency = 44'100.0 / double(fft.getSize()) * double(i / 2);
    amplitude = std::log10(std::abs(c)) * 20.0;
    phase = std::arg(c) * 360.0 / 2.0 / juce::MathConstants<double>::pi;
    groupdelay = -1.0 * phase_diff(phase, prev_phase) / ((frequency - prev_frequency) * 360.0);

    file << frequency << "," << amplitude << "," << phase << "," << groupdelay << "\n";

    prev_frequency = frequency;
    prev_phase = phase;
  }

Getestet mit einer Verzögerung um 16 Samples war die berechnete 
Gruppenlaufzeit in etwa korrekt. :)

viele grüsse
Zombie

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.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.