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


von Zombie (Gast)


Angehängte Dateien:

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)
1
0,-22.8513,0
2
689.063,-0.361587,-23.8445
3
1378.13,-2.39234,-39.5591
4
2067.19,-4.48345,-48.1488
5
2756.25,-6.32989,-52.5538
6
3445.31,-7.91214,-54.5978
7
4134.38,-9.27027,-55.2451
8
4823.44,-10.4469,-55.0161
9
5512.5,-11.4765,-54.2072
10
6201.56,-12.3857,-52.9975
11
6890.63,-13.1947,-51.4999
12
7579.69,-13.9191,-49.7896
13
8268.75,-14.5708,-47.9178
14
8957.81,-15.1593,-45.9208
15
9646.88,-15.6923,-43.825
16
10335.9,-16.1758,-41.65
17
11025,-16.6149,-39.4108
18
11714.1,-17.0136,-37.1189
19
12403.1,-17.3754,-34.7833
20
13092.2,-17.7032,-32.4113
21
13781.3,-17.9993,-30.0088
22
14470.3,-18.2657,-27.5806
23
15159.4,-18.5043,-25.1307
24
15848.4,-18.7164,-22.6625
25
16537.5,-18.9032,-20.1788
26
17226.6,-19.0659,-17.6821
27
17915.6,-19.2053,-15.1746
28
18604.7,-19.3221,-12.6583
29
19293.8,-19.4169,-10.1348
30
19982.8,-19.4902,-7.60583
31
20671.9,-19.5423,-5.07277
32
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

von Dergute W. (derguteweka)


Angehängte Dateien:

Lesenswert?

Moin,

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

Gruss
WK

von Zombie (Gast)


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

von Dergute W. (derguteweka)


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

von Zombie (Gast)


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
1
coef = exp(-2 * PI * fc / fs)

und dann für jedes Sample
1
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

von Dergute W. (derguteweka)


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

von Helmut S. (helmuts)


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
von Zombie (Gast)


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:
1
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

von Dergute W. (derguteweka)


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

von Zombie (Gast)


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:
1
double phase_diff(double a, double b) {
2
  double c = std::fmod(a - b + 180.0, 360.0);
3
  if (c < 0.0) {
4
    c += 360.0;
5
  }
6
  return c - 180.0;
7
}
und dann folgender Code für die Umrechnung
1
  std::ofstream file("fft.txt");
2
  double frequency = 0.0;
3
  double prev_frequency = 0.0;
4
  double amplitude = 0.0;
5
  double phase = 0.0;
6
  double prev_phase = 0.0;
7
  double groupdelay = 0.0;
8
  for (size_t i = 0; i < fft.getSize(); i += 2) {
9
    std::complex<float> c(arr[i], arr[i + 1]);
10
    frequency = 44'100.0 / double(fft.getSize()) * double(i / 2);
11
    amplitude = std::log10(std::abs(c)) * 20.0;
12
    phase = std::arg(c) * 360.0 / 2.0 / juce::MathConstants<double>::pi;
13
    groupdelay = -1.0 * phase_diff(phase, prev_phase) / ((frequency - prev_frequency) * 360.0);
14
15
    file << frequency << "," << amplitude << "," << phase << "," << groupdelay << "\n";
16
17
    prev_frequency = frequency;
18
    prev_phase = phase;
19
  }

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

viele grüsse
Zombie

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.