Probleme mit Berechnung der Phase

[3eepoint]

Hallo miteinander!

Ich dreh mich grade ein wenig im Kreis bezüglich der Phase bei meinen Simulationen. Die aus Comsol exportierten Datensätze zeigen nicht wirklich das, was ich erwarte an Phasenverhalten. Erwarten tu ich etwas in diese Richtung:



Das ist, was VituixCAD mir ausgibt. Die Phase in Comsol sieht aber eher so aus:


Nur bis 2 kHz berechnet, blau ist Phase des Schalldruck in 0.5 m, grün ist phase der Beschleunigung der Membrane und ich hatte eigentlich erwartet, dass die sich auf der Distanz ähnlicher sehen würden bzw. näher an Vituix dran ist...

Das Chassis ist in beiden Fällen das Selbe, genau so wie das Volumen. Das sind dann ja schon andere Ergebnisse und ich will Comsol da auch nicht blind vertrauen. Zudem stellt sich mir die Frage, wie sieht die "Natürliche" Pahse eines Chassis im geschlossenen Gehäuse eigentlich aus? Kennt jemand einen Typischen Messfehler der ähnliches verursacht der mir bei der Auswertung unterlaufen sein könnte?

[JFA]

Zudem stellt sich mir die Frage, wie sieht die "Natürliche" Pahse eines Chassis im geschlossenen Gehäuse eigentlich aus? Kennt jemand einen Typischen Messfehler der ähnliches verursacht der mir bei der Auswertung unterlaufen sein könnte?

Kennst du die Polygonapproximation? Da es sich um ein minimalphasiges System (schönen Gruß an den anderen Phasenthread) handelt, kannst du einfach folgende Regel ansetzen:
- Frequenzgang steigt mit 6 dB/Oktave (20 dB/Dekade) => Phase +90°
- Frequenzgang fällt mit 6 dB/Oktave (20 dB/Dekade) => Phase -90°
- Frequenzgang ist konstant => Phase 0°

Jede zusätzliche Ordnung (steigt/fällt 12, 18, 24, ... dB/Oktave) ergibt eine entsprechenden Faktor, also 180, 270, 360, ...

Dazwischen gibt es Übergänge, und da gibt es auch einen Zusammenhang mit der Güte, aber das habe ich mir nie gemerkt, weil es mich nicht interessiert hat.

[3eepoint]

Das hilft mir schonmal weiter, dann ergibt die Phase der Beschleunigung Sinn( ist ein sehr kleines Gehäuse daher wird der Schalldruck erst spät konstant). Jetzt muss ich nur noch die Akustische Phase hinbekommen, danke schonmal =)

[jones34]

Da du die Phase in 0,5m Entfernung mit Bezug zu deiner Anregung berechnest sind da 0,5m Delay drin.
Kannst du die raus machen? Dann gibts auch weniger Sägezahn.

[3eepoint]

Da du die Phase in 0,5m Entfernung mit Bezug zu deiner Anregung berechnest sind da 0,5m Delay drin.
Kannst du die raus machen? Dann gibts auch weniger Sägezahn.

Die Codezeilen dafür sind schon geschrieben, hatte aber noch keine Gelegenheit die zu testen.

[3eepoint]

So, es geht weiter! Die Randbedingungen für den Datenexport haben sich etwas geändert. ICh will an sich ja Spinorama daten auslesen. Das mache ich in 10 m Entfernung. Ich hab mir also eine Zeile geschrieben, die mir vom Mittelpunkt der Messung zu den Messpunkten die Distanz berechnet und auf Basis dessen die Phasenkompensation nach:

360*f*dt

Berechnet, wobei

dt = 343/Distanz

ist. Oder muss ich dafür die Distanz zum Chassismittelpunkt im Raum angeben? Wenn ja wie verhält sich das mit dem rückwertigen Schall der ja erst am Gehäuse vorbei muss?

So sieht das ganze dann in Vituix aus:



Jetzt hab ich spaßeshalber mal den Haken bei minimum Phase gesetzt und bekomme danach folgendes:



Das sieht ja nach dem aus was ich haben will, nur wie kommt vituix da drauf und ist das so richtig? Die Information scheint in den Daten an sich ja drin zu sein. Die Phasenberechnung in Comsol macht mich auch noch stutzig. Hier mal die Phase 10 mm vor der Membrane:



Da sollten Delays bei 200 Hz einen Phasenversatz von 20m° verursachen... das scheint hier ja auch nicht richtig zu sein. Wenigstens sehen sich die daten 10 mm und 10 m soweit ähnlich, dass ich denke das meine Kompensation halbwegs funktioniert.

[fosti]

....
Jetzt hab ich spaßeshalber mal den Haken bei minimum Phase gesetzt und bekomme danach folgendes:



.....

Passt doch: Hochpass 2. Ordnung 180° Phasendrehung!

Zu COMSOL......ein Kollege hat das "Multiphysics" mal in "Multifrustrix" umbenannt.....

[3eepoint]

Ja, passen tut das. Aber erst wenn ich in VituixCAD den Haken bei Minimal Phase setze und ich würde nun gerne wissen, wie es dazu kommt, damit ich es in Matlab selber machen kann. Sonst sind die Datensätze nur dafür brauchbar und das wäre schade und nicht ganz Sinn der Sache.

Frust ist kein Ausdruck für das was ich mittlerweile mit diesem Programm empfinde! Meine nächsten Projekt kommt aber zum glück weitestgehend ohne aus bzw. stützen sich auf das was ich schon habe

[jones34]

Man kann die Minimalphase mit der Hilbert Transformation berechnen.
https://en.wikipedia.org/wiki/Minimu...phase_response

Ob das jetzt die "richtige" Phase ist hängt davon ab ob du davor schon ein Minimalphasensystem hattest.
Bei einem einzelnen Chassis kommt man da in der Regel hin.
Ein Lautsprecher mit mehreren Chassis oder entsprechenden Filtern fällt da nicht mehr drunter.

Deine Comsol Simulation sieht schwer danach aus als wäre da ein Delay drin.
Schau dir mal das Groupdelay an. Das müsste für hohe Frequenzen ja eigentlich gegen 0 gehen.
So kann man das Delay recht einfach raus nehmen und den "richtigen" Phasengang zu erhalten.

Ich kenne Comsol nicht aber ich vermute stark das man bei der Darstellung der Phase irgendwie einen Bezug auswählen kann.
Hier könnte irgendwo das Problem liegen.

[Azrael]

Importiere ich z.B. Messungen, die ich in REW gemacht habe, nach VituixCAD, dann steckt die Laufzeit vom Messobjekt zum Mirkrophon noch drin, was den Phasengang ziemlich wild aussehen lässt. Wenn man sich die Minimalphase einblenden lässt, dann kann man den Phasengang mittels Einfügen eines negativen Delays für den Treiber, das der o.g. Laufzeit entspricht, an die Minimalphase anpassen (das hatte Slaughthammer/Onno hier mal gezeigt)

Mit meinem diesbezüglich allerdings doch recht laienhaften Wissen würde ich auch sagen, dass da bei dir noch Laufzeit welcher Art auch immer drinsteckt.

Viele Grüße,
Michael

[3eepoint]

Die Hilbert Transformation habe ich versuch, nur irgendwie komme ich da auch auf kein vernünftiges Ergebnis:


Die Hilbert trafo in Matlab wurde dafür mit dem komplex dargestellten Werte gefüttert. Matlab( bessergesagt ich) scheitert schon dran, aus den Real und Imaginärteil den Schalldruck wieder her zu stellen:



Gelb = real
Orange = Imaginär
Lila = abs(Complexe Darstellung)

Auch mit dem( deutlich smootheren) Kurven für den Schalldruck ohne komplexe Darstellung sieht der erste graph so aus.


Wie hoch währen denn für dich hoch von der Frequenz her? Das kostet im zweifelsfall alles Rechenzeit ==>Wenn möglich vermeiden =)

Das mit dem Delay leuchtet mir ja alles ein, nur 10 mm vor dem Chassis? Ich hatte ja auch eine Kompensationsfunktion geschrieben, die hat leider auch nicht viel gebracht wie in #6 zu sehen ist oder ist hier die Chassisposition gefragt?.... In den Comsol Beispielen ist auch die Phase mit drin in bezug auf eine Ebene Welle, die Ergebnisse sehen aber einfach gleich aus...

Warum ich da so drauf rumreite ist, dass ich gerne die Filter in Vituix machen würde und villeicht auch mal aus RePhase welche nutzen wollen würde, dass geht aber nur wenn ich die Amplitude und Phase vernünftig vorliegen habe.

[jones34]

Wie schon gesagt: Schaut dir bitte mal das Groupdelay an.
Wenn es da einen konstanten Offset gibt dann hast du irgendwo ein Delay das du nicht willst.
Der Phasengang sieht genau danach aus (linearer Phasenabfall).

Wenn du es nicht schon in Comsol raus bekommst, was imho gehen muss, kannst du das Delay auch nachträglich raus nehmen.
Zum Beispiel in Vituix CAD.

Zur Minimalphase:
Google spuckt für "calculate minimum phase from magnitude" als erstes das hier aus.
https://www.dsprelated.com/thread/26...itude-response

Ich versteh nicht genau was du meinst und kann dir nicht helfen wenn du deinen Code nicht postest.
Da du die Hilberttrafo anscheinend mit komplexen Werten fütterst geht hier schon was schief.



Die Zeile die ich mal frech von Wikipedia geklaut hab ist das was du implementieren musst.
Unglücklicherweise wurde die Systemfunktion H genannt und ist so leicht zu verwechseln mit der Hilberttrafo (Das schmucke H).
Weiterhin ist log() hier der natürlich Logarithmus und nicht der 10er.

[fosti]

.......
Gelb = real
Orange = Imaginär
......

Hat das einen Grund, warum Du das nach Real- und Imaginärteil aufsplittest, statt nach Amplitude und Phase oder wie schon vorgeschlagen nach Amplitude und GLZ?

[3eepoint]

@ Jones


Das Groupdelay, also dphi/dfreq sieht folgendermaßen aus:



Berechnet aus dieser Phasenlage:



Der Code dafür ist:

[phase_data_ver] = Phase_ver(model,x,y,z,angle; //Liest Phasendaten aus Comsol

for i = 1:length(angle)
temp = phase_data_ver(i,: ); //nimmt die erste Zeile die unter 0° gemessen wird
for n = 1:length(temp)-1
dphidf(n) = (temp(n)-temp(n+1))/(freq(n)-freq(n+1)); //berechnet dphi durch dfreq aus den Datenpunkten
end
end

Beides in 0° und 10 m vor dem Gehäuse

Ich würde das Delay gerne vor Vituix raus bekommen um einfach saubere Datensätze zu haben. Das Script von DSP Related ist das, was ich verwendet habe:


[ampl_data_hor] = SPL_hor(model,x,y,z,angle ); //liest die Daten aus Comsol, 73x116 Matrix

for i = 1:length(angle)

temp = ampl_ver(i,: ); // liest eine Zeile aus den ausgelesenen Daten
temp2 = [fliplr(temp) temp]; // Spiegelt diese um 0Hz herrum
min_phase = unwrap(-imag(hilbert(log(abs(temp2))))); // Performt die Hilbert Transformation
phase_data_ver(i,: )= min_phase(length(temp)+1:end)*180/pi; // Schreibt die unsymmetrierte Phase wieder in eine Matrix

end


H(jw)sollte doch die komplexe Übertragungsfunktion sein, daher hatte ich die komplexen Werte genommen. Das Beispiel von DSPrelated füttert da auch komplexe Werte rein.

@Fosti

Amplitude und Phase hatte ich ja zu Anfang, letztere passte aber nicht bzw. ich bekam das Delay nicht raus. Real- und Imaginärteil hatte ich dann genommen, weil ich wie ich oben schrieb dachte, dass die komplexe Übertragungsfunktion gebraucht wird. GLZ siehe oben.

[fosti]

nee groupdelay ist - (also Minus--->negatives Vorzeichen) Änderung der Phase pro Änderung der Kreis-(!!!)frequenz.....nicht nach der Frequenz......man kann das alles leider aus der mangelhaften Achs-Beschriftung nicht herauslesen ..... ja alles Aufwand......ich weiß

[3eepoint]

Also:

dphidfreq = -(temp(i)-temp(i+1))/(2*pi*freq(i)-2*pi*freq(i+1));

?

Dann bekomme ich:

[jones34]

H(jw)sollte doch die komplexe Übertragungsfunktion sein, daher hatte ich die komplexen Werte genommen. Das Beispiel von DSPrelated füttert da auch komplexe Werte rein.

Ja das ist die komplexe Übertragungsfunktion.
In die Hilbert Transformation zur Minimalphasenbestimmung geht aber der Betrag ein.
Dementsprechend muss du nicht erst dein Zahlenformat konvertieren wenn du schon Betrag und Phase hast.

Beim Groupdelay aufpassen.
Hast du denn eine Phase in rad?
Wenn deine Phase in ° ist dann wäre die korrekte Rechnung:

tau = -1/360° * dphi/df

Wenn ich mir deinen Phasengang so anschau sollte man da bei ~3ms raus kommen.

Dass das Groupdelay so zerhackt ist liegt an der Unterabstastung der Phase in der Simulation über 3kHz.

Könntest du zu deinen ganzen Importieren Daten dazu schreiben welche Einheiten hier vorliegen?
Die Hilberttrafo sollte natürlich lineare Werte bekommen und nix in dB.

[3eepoint]

Ok, auf ein neues!

-->Die Hilberttrafo sollte natürlich lineare Werte bekommen und nix in dB.

Gesagt, getan. Habe den Schalldruck in dB zu Pa konvertiert nach p0*10^(spl/20) und das anschließend in die Hilbert Transformation gegeben(code ist der Selbe wie vorher):


Phase: blau
SPL: orange

Das sieht schonmal etwas besser aus was die Werte an sich angeht, der Verlauf passt so aber denke ich noch nicht.

-->Könntest du zu deinen ganzen Importieren Daten dazu schreiben welche Einheiten hier vorliegen?
Die Werte liegen in dB und deg vor.

Bzgl. des Group delay, die Werte für die Phase waren in deg. Hier die Version mit rad bzw. der von dir angemerkten Änderung:



-->Dass das Groupdelay so zerhackt ist liegt an der Unterabstastung der Phase in der Simulation über 3kHz.

Was meinst du mit unterabtastung über 3 kHz? Sie Simulation geht nur bis 2 kHz. Sollten da mehr Stützstellen rein meinst du?

[jones34]

Die Werte deiner Hilberttrafo sehen doch schon gar nicht so schlecht aus.
Der Fehler bei tiefen Frequenzen kommt wohl daher das dein Frequenzgang da zu weniger Stützstellen hat und dann was kaputt geht.
Ich würde vor der Hilberttrafo den Freuquenzgang in eine bekannte Form bringen und die fehlenden Stellen interpolieren.

Da ich das eh demnächst machen hätte müssen hab ich das ganze mal schnell in Python geschmissen.

fr_response = import_response("stf.txt") #Importiert Frequenzgang als Array [frequency, amplitude(dB), phase(°)]
fnew = np.linspace(-fs/2, fs/2, fft) #Neue Frequenzachse

lin = to_lin_complex(fr_response) #Wandelt den Frequenzgang von Betrag/Phase zu komplexen Zahlen (nötig für meine Interpolationsfunktion)
intlin = interpolate(lin, fnew) #Interpoliert den Frequenzgang auf die neue Frequenzachse

hilbert = sig.hilbert(np.log(np.abs(intlin[1]))) #Hilberttrafo

minphase = to_ampl_phase_dB(intlin) #Bringt die interpolierte Darstellung wieder in die Form [frequency, amplitude(dB), phase(°)]
minphase[2] = np.rad2deg((-hilbert.imag + np.pi) % (2*np.pi) -np.pi) #Ersetzt die Phase durch die errechnete Minimalphase

Das liefert die folgenden Ergebnisse, verglichen mit der Beispielhaft verwendeten Messung aus ARTA.
Bitte verzeiht mir die sehr ungleiche Skalierung





Man sieht dass das ganz gut hin kommt aber bei sehr tiefen Frequenzen passiert Unfug.
Bei meinem Beispiel ist das kein Problem da der Bereich in dem das passiert gar nicht relevant ist.
Bei dir passiert das früher, vermutlich geht da was mit den Datenpunkten um die Null schief.
Deine Skalierung ist auch etwas arg und verzerrt das fürs Auge.

Was meinst du mit unterabtastung über 3 kHz? Sie Simulation geht nur bis 2 kHz. Sollten da mehr Stützstellen rein meinst du?

Ich meinte 300Hz, hab mich in der Potenz verguckt. Damit sollte das Groupdelay auch ca. 30ms sein, was bei dir ja jetzt zu passen scheint.
Mehr Stützstellen hilft nur bedingt, da es das Problem nur zu höheren Frequenzen verschiebt.
Im Endeffekt bleibt dir nur das Problem in Comsol zu fixen.
Wenn du das Groupdelay schätzen kannst, dann kannst du die "kaputten" Stützstellen wieder zurück biegen.

Edit: Es lohnt sich evtl. deine ganzen Funktionen an generierten Testdaten zu probieren die in sich schlüssig sind und von denen du die entsprechenden Ergebnisse schon kennst.

Edit2: Die Gruppenlaufzeit kommt doch ganz gut mit deinen 10m Distanz hin, wenn du die raus rechnest sollte das doch passen.

Edit3: Meine Beispielübertragungsfunktion zeigt auch die Beschränkungen der Minimalphase. Gezeigt ist eine gefügte Nahfeldmessung eines 8" PA Topteils.
Bei tiefen Frequenzen ist die Annahme der Minimalphasigkeit korrekt aber sobald man den Bereich in dem alles deutlich kleiner als die Wellenlänge ist verlässt macht man potenziell Fehler.

Hier mal Phase vs. Minimalphase. Die Abweichung nach unten hin ist wohl eine Kombination aus Messfehler und Artefakten bei der Berechnung.

[3eepoint]

Hmmm.... in dem Post von DSP Related schreib der eine ja auch, dass es um 0 Hz rum interpoliert werden muss damit die Hilbert Transformation keine Zicken macht. Ich mach die Simulation nochmal, startend bei 0 Hz und schau ob es dann besser wird.

Die Kompensation mittels group delay und dann wieder rausrechnen ist mir etwas wackelig. Sowas hatte ich ja in Post #6 schon probiert, bin mir da aber mit den Positionen sehr unsicher. Bzw. Welchen Wert aus meinem Group delay nehme ich denn zur Kompensation? einfach 30 ms abziehen, wieder in die Phase umrechnen und gut ist? Das ist mir noch nicht klar...