Unterschiede LTspice - Octave

[Andreas B. (abeld)]

Hallo,

als Maschinenbauingenieur habe ich nur lückenhafte Kenntnisse der 
E-Technik. Als Rentner habe ich begonnen mich mit Gammaspektrometrie / 
Szintillation zu beschäftigen. Einen Adapter, der die Impulse aus dem 
Photomultiplier eines Detektors so in die Länge streckt, das sie mit 
einer Soundkarte digitalisiert werden können, habe ich mir selbst 
zusammengelötet. Vorlage war ein PMT-Adapter von 
https://www.theremino.com/de/blog/gamma-spectrometry (open Hardware). 
Von Theremino stammt auch die Open-Source-Software Theremino MCA, mit 
der aus den Daten ein Histogramm erzeugt wird.

Das Hardware-Design ist nun aber auch schon über 10 Jahre alt und 
enthält Bauteile, die nicht mehr aktuell sind und nur mit Problemen zu 
bekommen sind (BC237, BSP300, MJE13003). Daher habe ich begonnen mich 
mit dem Design zu beschäftigen, um es evtl. aktualisieren zu können.

Anbei 2 LTspice-Simulationen (Transient und AC Analysis). In der Datei 
Audacity_Impuls_kurz_32PCM.wav befindet sich ein digitalisierter Impuls. 
Die Parameter für PULSE sind die beste Näherung an diesen Impuls, die 
ich finden konnte.

Ich möchte daher aus dem digitalisierten Output-Impuls den dazugehörigen 
Input-Impuls berechnen.

Das geht wohl mit LTspice nicht!? Daher habe ich die Simulationsdaten im 
Zeit- und im Frequenzbereich als Text abgespeichert und mit dem Editor 
so bearbeitet, das ich sie als csv-Datei in Octave einlesen kann.

Mit Octave ließen sich dann aus den LTspice-Werten Pole und Nullstellen 
für Übertragungsfunktionen herausrechnen. Dazu sind die Funktionen zpk 
und bode aus dem Octave-Paket „Control“ nötig.

Für den Output ergibt sich damit als Näherung:
G(s) = (2.4e+18 s^2) / (s^6 + 1.39e+06 s^5 + 1.073e+11 s^4 + 1.343e+15 
s^3 + 3.328e+18 s^2 + 1.355e+20 s + 4.004e+20)

Die Pole liegen dabei bei den Frequenzen:
1.308321137056221e+06
6.637385825993693e+04
1.155793960442228e+04
3.260724322521371e+03
3.812998621639556e+01
3.208606224259006e+00

Um aber die Octave-Berechnungen an die LTspice-Werte heran zu bringen, 
musste ich tricksen.

1. um eine Steigung von 20dB/Dekade hinzubekommen musste ich die 
Octave-Werte mit 20*log10() anrechnen. Eigentlich müsste es doch mit 
10*log10() richtig sein?!

2. An den ersten beiden Netzpunkten (n1 uns n2) konnte ich noch mit 
einer Verstärkung von 1 arbeiten, ab dann waren große Verstärkungen 
nötig um an die LTspice-Werte heran zu kommen: n3: 1500, n4: 2300000, 
n5: 7.75*10^13 n6 und output: 2.4*10^18.

3. Bei n5, n6 und für den output waren die Phasenwerte von Octave um 
180° höher, als bei LTspice.

4. Um aus den Übertragungsfunktionen den Output im Zeitbereich zu 
berechnen musste ich ab Netzpunkt „n3“ eine 8-fach längere Zeit 
eingeben, sie aber gerafft ausgeben.

5. Ab „n5“ stimmt das Vorzeichen nicht mit LTspice überein.
Kann mir Jemand erklären, woher die Unterschiede zwischen Octave und 
LTspice kommen?

Gruß Andreas

[Mark S. (voltwide)]

Auch wenn ich im Einzelnen hier überfragt bin - auf jedenfall
chapeau, monsieur!

[Achim M. (minifloat)]

Andreas B. schrieb:
> um eine Steigung von 20dB/Dekade hinzubekommen musste ich die
> Octave-Werte mit 20log10() anrechnen.

Da es um Leistungswurzelgrößen bzw. Feldgrößen 
(https://de.m.wikipedia.org/wiki/Leistungsgr%C3%B6%C3%9Fe) wie z.B. 
Spannung geht, ist das mit der 20 schon korrekt.

mfg mf

[Rainer W. (rawi)]

Andreas B. schrieb:
> 3. Bei n5, n6 und für den output waren die Phasenwerte von Octave um 180°
> höher, als bei LTspice.

180° höhere Phasenwerte bedeutet eine Invertierung. Ohne jetzt im Detail 
nachgerechnet zu haben, hört sich das für mich nach einem 
Vorzeichenthema an.

zu 5:
Das Vorzeichen kann an der Einbaurichtung hängen.

[Andreas B. (abeld)]

Hallo,

auf diesen Unterschied von Feld- und Energiegrößen war ich noch nie 
gestoßen. Ich war nur verwundert, das mir die Diagramme von Octave 
Steigungen von 10dB/Dekade zeigten, wo ich doch in allen Lehrbüchern 
immer Diagramme mit 20dB/Dekade gesehen hatte.

Frage nun: Spielt dieser Unterschied auch in anderen Octave-Berechnungen 
eine Rolle?

Zu den Berechnungen im Zeitbereich (Punkt 4 in meinem ersten Beitrag):

So sieht das im Skript aus:
fn1 = lsim(xsys1,fin,t2(t2max));
fn2 = lsim(xsys2,fin,t2(t2max));
fn3 = lsim(xsys3,fin,8*t2(t2max));
fn4 = lsim(xsys4,fin,8*t2(t2max));
fn5 = lsim(xsys5,fin,8*t2(t2max));
fn6 = lsim(xsys6,fin,8*t2(t2max));
fou = lsim(xsys7,fin,8*t2(t2max));

t2(t2max)=6e-4 also 0,6ms oder 600us

Warum ist das für n1 und n2 richtig und ab n3 nicht?

Alle fnx werden mit dem Zeitvektor t2 geplottet. fn1 und fn2 also mit 
der „richtigen“ Zeitbasis, ab fn3 werden dann 8*0,6ms=4,8ms auf 0,6ms 
zusammengestaucht. Anbei nochmal die Zeitbereichsdiagramme im png-Format 
direkt aus Octave abgespeichert (die im ersten Beitrag waren als jpeg 
abgespeichert und dann konvertiert worden).

Gruß Andreas

[Andreas B. (abeld)]

Hallo,

ich weiß leider immer noch nicht, wie ich mit einer Übertragungsfunktion 
und einem Output den Input errechnen kann.

Es gibt im Octave-Paket „Control“ zwar die Funktion inv(), die eine 
Übertragungsfunktion invertiert. Die vertauscht Zähler und Nenner, macht 
also aus Polen Nullstellen und umgekehrt. Mit dieser invertierten 
Übertragungsfunktion lässt sich aber nicht mit lsim() weiterrechnen. 
Lsim() bricht mit einer Fehlermeldung ab, weil:
     LTI model. System must be proper, i.e. it must not have more zeros 
than poles.

Also vorwärts geht es:
sys = zpk([0 -zero1],[-pol1 -pol2],1); # bilde Übertragungsfunktion aus
                                         Polen und Nullstellen
[mag, pha] = bode (sys,freq); # berechne Magnitude und Phase
                                (Frequenzbereich)
fout = lsim(sys,fin,tmax); # berechne Output aus Input (Zeitbereich)

Auch das noch:
isys=inv(sys); # Übertragungsfunktion invertieren

Aber das nicht mehr:
fin=lsim(isys,fout,tmax); # berechne Input aus Output (Zeitbereich)

Es gibt zwar ein Octave-Paket „symbolic“, das neben der 
Laplace-Transformation auch die inverse Laplace-Transformation 
beherrscht, das aber nur symbolisch und eben leider nicht für diskrete 
Funktionen.

Hat Jemand einen Tipp für mich, wie ich mit der Lösung meines Problems 
weiter komme?

Gruß Andreas

[Christoph M. (mchris)]

Hast Du schon mal ChatGPT gefragt?:

> "dein text" + wie kann man das Problem in "Octave" lösen

pkg load control;

% Beispiel: Übertragungsfunktion G(s)
num = [1];      % Zähler von G(s)
den = [1, 2, 1]; % Nenner von G(s)
G = tf(num, den);

% Output-Signal (angenommen)
t = 0:0.01:10;
y = sin(t);  % Beispiel-Output

% Inversion von G(s)
Gi = tf(den, num);  % inverse Übertragungsfunktion

% Proper machen (optional)
epsilon = 1e-6;     % kleiner Term zur Stabilisierung
Gi_proper = tf(den, num + epsilon);

% Berechnung des Inputs
u = lsim(Gi_proper, y, t);

% Plotten
figure;
subplot(2,1,1);
plot(t, y); title('Output y(t)');
subplot(2,1,2);
plot(t, u); title('Input u(t) berechnet');

[Giovanni (sqrt_minus_eins)]

Hier mein Stand der Ermittlungen:

Habe versucht aus der Simulation im Zeitbereich eine 
Übertragungsfunktion (1->6) zu ermitteln.

Bild: Kleine Frequenzen sind falsch, da Zeitbereich nur 600µs

G = TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
    -1.4825867447070777s^2 + 3.915381717456055e7s - 4.5138075195136e14
---------------------------------------------------------------------------
1.0s^3 + 354999.27666149335s^2 + 2.74055777345083e10s + 5.38145309638144e14

zu Inversion: Die Ordnung im Zähler von G(s) ist um 1 höher als die 
Ordnung im Nenner.
Dreht man G(s) um, dann muss man differenzieren. Geht das? Vielleicht 
braucht man dann noch ein PT1?

z.B. Derivative aus MODELICA hat noch ein T - sonst geht es nicht.

This blocks defines the transfer function between the input u and the output y as approximated derivative:

             k * s
     y = ------------ * u
            T * s + 1

[Giovanni (sqrt_minus_eins)]

Giovanni schrieb:
> zu Inversion: Die Ordnung im Zähler von G(s) ist um 1 höher als die
> Ordnung im Nenner.
ist nicht richtig.

KORREKTUR:
zu Inversion: Die Ordnung im Zähler von G(s) ist um 1 **kleiner** als 
die Ordnung im Nenner.

[Andreas B. (abeld)]

Hallo Christoph,

bei mir antwortet lsim() mit der Fehlermeldung
   this descriptor system cannot be converted to regular state-space 
form
?!

Gruß Andreas

[Giovanni (sqrt_minus_eins)]

jetzt der Versuch mit Filter 2. Ordnung.

wie man die optimale Filter-Zeitkonstante bestimmt, ist mir aktuell 
nicht klar. Aber das rekonstruierte Eingangssignal hat zumindest mit 
eine Ähnlichkeit mit dem Original.

    filter = tf(1,[T^2,T,1])
    G1 = 1/G
    Ginv = G1 * filter
    Gd = c2d(Ginv,deltaT)
    GdSS = ss(Gd)
    result = lsim(GdSS, y')  # y' = output

Details dazu:

   T = 1.5e-5
   deltaT = 4.995304547350856e-7  # sample
   800 samples im Bereich 0.0µs<=time<=400µs


G1 = 
1.0s^3 + 2.7920053056435944e6s^2 + 1.6294937913486914e11s + 3.6924968356864e15
------------------------------------------------------------------------------
    -170.02790203411132s^2 - 1.3058370031801758e9s - 2.946599056677888e15

filter = 
              1.0
------------------------------------
2.250000e-10s^2 + 1.5000e-5s + 1.0


Ginv = 
                    1.0s^3 + 2.7920053056435944e6s^2 + 1.6294937913486914e11s + 3.6924968356864e15
----------------------------------------------------------------------------------------------------------------------
-3.825627795767506e-8s^4 - 0.29636374424605133s^3 - 682742.3707022617s^2 - 4.55048228533485e10s - 2.946599056677888e15

[Andreas B. (abeld)]

Hallo Giovanni,

anliegendes zu Deinem vorletzten Beitrag, den letzten muss ich mir erst 
noch genauer ansehen.

für den Netzpunkt “n6“ ermittelt mein Skript die Pole
   1.308314590210393e+06
   6.636831434701991e+04
   1.155467770258950e+04
   3.256208624016595e+03
   3.098415825682461e+01

Mit einer Nullstelle bei s=0 ergibt sich die Übertragungsfunktion

                           2.4e+18 s
y1: ------------------------------------------------------------------
    s^5+1.39e+06 s^4+1.073e+11 s^3+1.341e+15 s^2+3.308e+18 s+1.012e+20

Wie kommst Du zu Deiner Übertragungsfunktion?

Mit 600us und 1048 Stützstellen ergibt sich eine Samplingrate von 
1,746MHz. Zur Zeit arbeite ich mit einem LTspice-Datensatz mit 6010 
Stützstellen und bekomme keine erkennbar anderen Ergebnisse. Wo siehst 
Du eine Begrenzung durch die 600us? Hilft es diese Zeit in LTspice zu 
verlängern?

Für mein Problem spielt die Filterwirkung im Frequenzbereich eigentlich 
keine Rolle. Ich muss ja meinen Input-Impuls in die Länge ziehen, damit 
ich ihn mit der Soundkarte digitalisieren kann. Ich habe aber leider 
noch keine Kenngröße gefunden, die dieses „in die Länge ziehen“ 
beziffert. Die Phase im Bode-Diagramm beschreibt ja so etwas wie einen 
zeitlichen Verzug, ich brauche so etwas aber für meinen Impuls im Ganzen 
und nicht für einzelne Frequenzen. Mit einer Fourier-Zerlegung meines 
Impulses bin ich auch nicht weitergekommen, mit einer solchen hätte ich 
dann ja einige diskrete Frequenzen bei denen die Amplitude hoch ist.

Hat Jemand eine Idee, wie ich da weiterkomme?

Gruß Andreas

[Andreas B. (abeld)]

Hallo,

ich hab Mal mit den Skripten aus dem „Control“- und dem „Symbolic“-Paket 
gespielt.

Bei Übertragungsfunktionen mit vielen Polen erhalte ich Fehlermeldungen 
und Abbrüche. Drumm hier Mal mit der Übertragungsfunktion am Netzpunkt 
„n2“:

Nullstellen: 0.0 und 2685.839972740719
Pole: 3.227696673321165e+01 und 4.688862094194152e+04

        s^2 + 2686 s
y1: -------------------------
    s^2+4.692e+04 s+1.513e+06

Diese Übertragungsfunktionen lassen sich nicht direkt weiterverwenden, 
sondern müssen händisch umgeformt werden.

pkg load symbolic
syms s t
iG=(s^2+4.692e+04*s+1.513e+06)/(s^2+2686*s);
ilaplace(iG)
                   -2686*t
  44500   3449986*e
  ----- + ----------------
    79           79

Das muss dann wieder händisch umgesetzt werden um damit weiterarbeiten 
zu können.

iGt= 44500/79+3449986*exp(-2686.*t)/79;

Das ergibt eine recht einfache Zeitfunktion, die aber kein erwartetes 
Ergebnis liefert. Mit:

G=(s^2+2686*s)/(s^2+4.692e+04*s+1.513e+06)
ilaplace(G)

ergibt sich eine recht undurchsichtige Formel, die ich so interpretieren 
würde:

fout=(-3570325193*exp(20*sqrt(5488586).*t‘)+1523848*sqrt(5488586)*exp(20 
*sqrt(5488586).*t‘)-3570325193*1523848*  sqrt(5488586))* 
exp(-10*(sqrt(5488586)+2346).*t)/1611429

Das Ganze führt aber – auch nach einer Multiplikation mit der 
Input-Funktion – zu keinen vernünftigen Ergebnissen.

Oder mach ich da was falsch?

Gruß Andreas

[Giovanni (sqrt_minus_eins)]

Andreas B. schrieb:
> Wie kommst Du zu Deiner Übertragungsfunktion?

Ich bin etwas in Verzug. Musste erst LTspice installieren um die 
Schaltung einmal zu sehen, dann habe ich versucht das Control PKG für 
octave zu installieren, bin aber dann wieder bei den mir vertrauten 
Werkzeugen gelandet.

Muss gestehen, dass ich noch nicht alle Beiträge gelesen habe.
Die Aufgabenstellung: Ein Signal, das durch eine Schaltung geht an Hand 
vom Ausgang wieder zu rekonstruieren. Dazu gibt es Simulation im Zeit- 
und im Frequenzbereich.

ich verwende jetzt JULIA mit 
https://baggepinnen.github.io/ControlSystemIdentification.jl/dev/

Bisher:
* Aus dem Zeitbereich die Übertragungsfunktion (in->out) ermittelt
* Versuch einer Umkehr (out->in), siehe oben

coming soon:
* Ermittlung Übertragungsfunktion aus Messungen im Frequenzbereich
* Übertragungsfunktion der einzelnen Blöcke (n1->n2), (n2->n3), ... mit 
Dokumentation
* Versuche mit octave ??

Andreas B. schrieb:
> Hallo,
>
> als Maschinenbauingenieur habe ich nur lückenhafte Kenntnisse der
> E-Technik. Als Rentner habe ich begonnen ....

PS: als Energietechniker habe ich nur lückenhafte Kenntnisse der
Regelungstechnik. Als Rentner habe ich begonnen mich mit JULIA zu 
beschäftigen.

giovanni

[Giovanni (sqrt_minus_eins)]

hier der 1. Teil mit den Übertragungsfunktionen V1 -> Vn[3,4,5,6].
Der gesamte Zeitbereich war 500µs mit 800 samples
Ich habe auch höhere Ordnungen versucht - bringt keine Verbesserung.

im PDF sind alle Daten für Continuous, Discrete und StateSpace - mit der 
angegebenen Schrittweite!!

[Andreas B. (abeld)]

Hallo Giovanni,

anhand des Bodediagramms von LTspice hatte ich schon vermutet, das es 
sich bei der Übertragungsfunktion der gesamten Schaltung um eine 
Funktion mit mehreren Nullstellen und mehreren Polen handeln wird. Zur 
Herleitung wollte ich Schritt für Schritt durch die Schaltung gehen um 
dabei immer komplexer zu werden. daher Input->n1 dann Input->n2 ... bis 
Input->Output. Die Zwischenschritte n1->n2 hab ich mir nicht angeschaut.

Überrascht war ich, das die Pole/Nullstellen aus der ersten 
Übertragungsfunktion nicht genau so in der zweiten wieder vorkommen.

In Octave habe ich dann nur reale Pole/Nullstellen gefunden. Ich hatte 
zwar schon gesehen, das man auch mit imaginären Polen/Nullstellen 
arbeiten kann. Meine Versuche die Octave-Funktionen durch imaginäre 
Anteile näher an die LTspice-Funktionen heranzubekommen waren aber nicht 
von Erfolg gekrönt.

Gruß Andreas

[Andreas B. (abeld)]

Hallo Giovanni,

hier Mal die Übertragungsfunktionen V(n2)/V(n1), V(n3)/V(n2), 
V(n4)/V(n3), V(n5)/V(n4) und V(n6)/V(n5)

V(n2)/V(n1)
Nullstelle:   2.771e+3
Pol:     45.29e+3
V(n3)/V(n2)
Nullstelle:   keine
Pol:     2.690e+3
V(n4)/V(n3)
Nullstelle:   keine
Pol:     7.245e+3
V(n5)/V(n4)
Nullstellen:  31e+3, 900e+6, 4.5e+9, 600e+9
Pole:    30.95e+9, 1.108e+06, 6.5e+03, 1.0e+12
Alle Nullstellen, sowie die beiden hintersten Pole habe ich „per Hand“ 
iterativ ermittelt!
V(n6)/V(n5)
Nullstelle:   178.1e+6
Pole:    27.11e+3, 2.043e+09

Die Ermittlung der Übertragungsfunktionen erfolgte so:
xsys1 = zpk([-zer(1,1)],[-pol(1,1)],1);
[mag1, pha1] = bode (xsys1,freq);
xsys2 = zpk([],[-pol(2,1)],1496);
[mag2, pha2] = bode (xsys2,freq);
xsys3 = zpk([],[-pol(3,1)],1469);
[mag3, pha3] = bode (xsys3,freq);
xsys4 = zpk([-zer(4,1) zer(4,2) -zer(4,3) -zer(4,4)], ...
            [-pol(4,1) -pol(4,2) -pol(4,3) -pol(4,4)],0.3722);
[mag4, pha4] = bode (xsys4,freq);
xsys5 = zpk([-zer(5,1)],[-pol(5,1) -pol(5,2)],3.042e+5);
[mag5, pha5] = bode (xsys5,freq);

Wobei zu Fragen ist: Wieso musste zer(4,2) ein anderes Vorzeichen 
bekommen?

Gruß Andreas

[Giovanni (sqrt_minus_eins)]

Hallo Andreas,

Ich denke das "Problem" liegt darin, dass die Eigenwerte des Systems tw. 
conj-Complex sind. Der Grund dafür liegt in der Schaltung (muss ich noch 
anschauen).
Beispiel n1->n6:  Die Näherung (ermittelt im Frequenzbereich) passt 
recht gut.

TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}

                          -3.34880e21
-------------------------------------------------------------------
1.0s^4 + 7.69952e6s^3 + 2.952016e12s^2 + 1.801918e17s + 4.997285e21

Eigenwerte  (EW's): eigen(ss(G).A).values/2pi  = 
 -1.161580611485247e6 + 0.0im
   -53080.16624805282 + 0.0im
   -5378.688450414837 - 4803.464208327636im
   -5378.688450414837 + 4803.464208327636im

Dies ist auch der Grund warum man mit einer Näherung mit realen Werten 
für die Pole nicht ans Ziel kommt.
Es gibt nur eine Schaltung. D.h. Die EW's müssen mehr oder weniger in 
allen Signalen zu sehen sein.

PS: die Funktion "tfest" gibt es NUR in MATLAB-Control. OCTAVE-Control 
und Python-Control kennt diese Funktion nicht.

https://de.mathworks.com/help/ident/ref/tfest.html

grüße giovanni

[Andreas B. (abeld)]

Hallo,

ich hab in Octave bisher nur mit Nullstellen und Polen gearbeitet, die 
keinen Imaginärteil hatten. Kann man aus einem Bodediagramm auch 
Nullstellen/Pole mit Imaginärteil herauslesen?

Ich hatte schon mit Nullstellen/Polen mit Imaginärteil experimentiert. 
Meine Vermutung war, das der Absolutwert einer imaginären Stelle dem 
Wert der Frequenz der Stelle im Bodediagramm gleichen müsse. Das ist 
nicht so, dabei wandert die Phase auf der Frequenzachse! Wie ist es 
dann?

An der Übertragungsfunktion V(n5)/V(n4) habe ich mich mit einer 
Kurvendiskussion versucht. Weil die Frequenzen des Bodediagramms nicht 
äquidistant verteilt sind, habe ich mir extra eine Octave-function 
geschrieben, die diskrete Funktionen mit nicht-äquidistanten 
Stützstellen ableiten/differenzieren kann (Anlage). Das Ergebnis (siehe 
Anlage) sieht mir nicht plausibel aus, weil an horizontalen Stellen der 
Funktion deren Ableitung nicht Null ist. Was mach ich da falsch? Ich 
hatte auch schon versucht aus den dB-Werten der Magnitude mit 
10^(dB-Wert/20)=Verstärkung zu arbeiten. Das klappt auch nicht.

Gruß Andreas

[Andreas B. (abeld)]

Hallo Giovanni,

ich hab Dein Beispiel mal in ein Octave-Skript umgesetz
pol1 = -1161580.611485247
pol2 = -53080.16624805282
pol3 = -5378.688450414837 - 4803.464208327636i
pol4 = -5378.688450414837 + 4803.464208327636i
mit
xsys1 = zpk([],[pol1 pol2 pol3 pol4],-3.34880e21);
[mag1, pha1] = bode (xsys1,freq);

Ergib sich bei mir die Übertragungsfunktion zu:
                               -3.349e+21
 y1:  -------------------------------------------------------------
      s^4 + 1.225e+06 s^3 + 7.478e+10 s^2 + 7.264e+14 s + 3.206e+18

Wenn ich meine Polsuche darauf ansetze, ergeben sich bei mir folgende 
nicht-imaginäre Polstellen
Polstelle  1:  1.28325e+06
Polstelle  2:  65876.8
Polstelle  3:  11380.6
Polstelle  4:  3332.99

In der Anlage ein Vergleich beider Übertragungsfunktionen.

Gruß Andreas

[Andreas B. (abeld)]

Hallo Giovanni,

hatte in meinem Skript noch einen Fehler. Anbei ein korrigiertes 
Zeit-Diagramm.

Und hier auch noch die Übertragungsfunktion mit den rein realen Polen.

                               -3.349e+21
 y1:  -------------------------------------------------------------
      s^4 + 1.364e+06 s^3 + 1.044e+11 s^2 + 1.295e+15 s + 3.207e+18

Gruß Andreas

[Andreas B. (abeld)]

Hallo Giovanni,

hab mir Mal die Unterschiede zwischen Deiner und meiner 
Übertragungsfunktion angesehen:

>> pol(1,:)'
ans =
  -1.161580611485247e+06 -                     0i
  -5.308016624805282e+04 -                     0i
  -5.378688450414837e+03 + 4.803464208327636e+03i
  -5.378688450414837e+03 - 4.803464208327636e+03i
>> pol(1,1)*pol(1,2)*pol(1,3)*pol(1,4)
ans =  3.206377748151160e+18 + 2.560000000000000e+02i
>> abs(pol(1,1))*abs(pol(1,2))*abs(pol(1,3))*abs(pol(1,4))
ans = 3.206377748151160e+18
>> pol(1,1)+pol(1,2)+pol(1,3)+pol(1,4)
ans = -1225418.154634130
>> abs(pol(1,1))+abs(pol(1,2))+abs(pol(1,3))+abs(pol(1,4))
ans = 1229083.476211367
>> pol(1,1)*pol(1,2)+pol(1,1)*pol(1,3)+pol(1,1)*pol(1,4)+pol(1,2)*pol(1,3)+ 
pol(1,2)*pol(1,4)+pol(1,3)*pol(1,4)
ans = 74775459318.73108

Für mich ist die meine plausibel. Kannst Du den Unterschied zu Deiner 
erklären?

Gruß Andreas

[Andreas B. (abeld)]

Hallo,

nun hab ich auch das mit den Ableitungen hinbekommen (Anlage). Zum 
Plotten müssen die Ableitungen mit der Frequenz multipliziert werden und 
bei der 2. Ableiten aus der 1. muss diese auch mit der Frequenz 
multipliziert werden.

ma1=diffqt(magnitude,freq)';
ma2=diffqt(freq.*ma1,freq)';
ph1=diffqt(phase,freq)';
ph2=diffqt(freq.*ph1,freq)';

semilogx(freq,magnitude,"linewidth",1,"r", ...
         freq,freq.*ma1,"linewidth",1,"linestyle",":","m", ...
         freq,freq.*ma2,"linewidth",1,"linestyle","-.","b")

semilogx(freq,phase,"linewidth",1,"r", ...
         freq,freq.*ph1,"linewidth",1,"linestyle",":","m", ...
         freq,freq.*ph2,"linewidth",1,"linestyle","-.","b")

freq sind die Frequenzen der Stützstellen, die nicht-äquidistant sind.
magnitude, phase und freq kommen aus dem Datenexport einer 
LTspice-Simulation.

Gruß Andreas

[Giovanni (sqrt_minus_eins)]

Hallo Andreas,

sorry, hatte noch keine Zeit um mir die OCTAVE-scripts anzusehen. 
OCTAVE+Control habe ich aber schon installiert.

Habe versucht die Schaltung zu simulieren um die Übertragungsfunktion 
(ABCD) zu bekommen. Leider habe ich keine korrekten Transistormodelle, 
aber das Ergebnis ist soweit OK [Modelica-Fn0n6.png]. Es treten durch 
den Verstärker tatsächlich komplexe Eigenwerte auf. Das wollte ich 
bestätigt haben.

ALLE Simulationen jetzt vom Eingang (n0) zum Ausgang (n6).

Die komplexen EW's kann man mit Variation der Transistorparameter (die 
ich nicht kenne) verändern.

Eigenwerte aus MODELICA Übertragungsfunktion: 
7-element Vector{ComplexF64}:
 -4.65332334802824e6 + 0.0im
 -52353.003162961424 + 0.0im
  -8857.108506193686 - 3493.9192360511533im
  -8857.108506193686 + 3493.9192360511533im
 -1123.6518096924808 + 0.0im
 -31.300341804948534 + 0.0im
    -3.9598292898365 + 0.0im

Dann noch eine verbesserte Näherung mit OPTI n0->n6 [Fn0n6.png]. Für 
mich OK, um hier weiterzumachen.

Eigenwerte auf Basis LTspice Übertragungsfunktion + Näherung.
5-element Vector{ComplexF64}:
 -1.1795034187874736e6 + 0.0im
   -52513.231798263354 + 0.0im
    -5668.055982862397 - 5140.836713850675im
    -5668.055982862397 + 5140.836713850675im
   -29.886406763793637 + 0.0im

Werde jetzt versuchen die Daten/Parameter in einem Format zu speichern, 
um mit OCTAVE weitermachen zu können.

To whom it may concern: Wäre interessant das Problem mit MATLAB+Control 
(tfest) zu simulieren. Muss aber nicht sein.


giovanni

[Andreas B. (abeld)]

Hallo Giovanni,

der BC237 wird nicht mehr produziert. Ich konnte ihn nur noch aus China 
bestellen (Restposten?). Für Octave habe ich keine Parameter von ihm. In 
den aktuellen Standardbibliothek von LTspice gibt es ihn auch nicht 
mehr, es lassen sich aber ältere Bibliotheken im Internet finden, die 
ihn enthalten. In meinem LTspice-Schema verende ich die Zeile:

.MODEL BC237B NPN(IS=1.8E-14 ISE=5.0E-14 ISC=1.72E-13 XTI=3 BF=400 
BR=35.5 IKF=0.14 IKR=0.03 XTB=1.5 VAF=80 VAR=12.5 VJE=0.58 VJC=0.54 
RE=0.6 RC=0.25 RB=0.56 CJE=13p CJC=4p XCJC=0.75 FC=0.5 NF=0.9955 
NR=1.005 NE=1.46 NC=1.27 MJE=0.33 MJC=0.33 TF=0.64n TR=50.72n EG=1.11 
VCEO=45V ICRATING=100M MFG=ZETEX)

Ich hab aber keine Ahnung was das alles bedeutet. Ein Datenblatt gibt es 
von onsemi.

Die LTspice-Schemata gibt es im zip-Archiv in meinem ersten Beitrag. Ich 
leg Mal Bilder der Schemata bei.

Gruß Andreas

[Andreas B. (abeld)]

Hallo Giovanni,

ich hab noch Mal einen LTspice-Datensatz bis 100THz erzeugt, um auch die 
Pole bei höheren Frequenzen sehen zu können. Für n6 findet mein 
Octave-Skript nun die Pole:

   3.943523323388408e+10
   1.308313771881181e+06
   6.636830580309634e+04
   1.155466989750578e+04
   3.256208464970547e+03
   3.098414922897306e+01

Hier die Übertragungsfunktion, den Verstärkungsfaktor hab ich noch nicht 
bestimmt:
                                             s
------------------------------------------------------------------------
s^6 + 3.944e+10 s^5 + 5.48e+16 s^4 + 4.23e+21 s^3 + 5.289e+25 s^2
                                   + 1.305e+29 s  + 3.992e+30

Wenn ich wüsste, wie ich aus den LSspice-Bode-Daten auch imaginäre 
Anteile ermitteln kann, würde ich das auch versuchen. Wo kann ich 
nachlesen, wie das geht?

Gruß Andreas

[Andreas B. (abeld)]

Hallo,

ich hab mir Mal 6 Übertragungsfunktionen aus Nullstellen und Polen 
zusammengestellt, auch mit komplexen Stellen. Ich habe aber noch keinen 
Weg gefunden um aus den diskreten Funktionsdaten die komplexen Stellen 
wieder herauszuholen.

Gruß Andreas

[Kai D. (Firma: CAD- und CAE-Konstruktion) (robokai)]

Mark S. schrieb:
> Auch wenn ich im Einzelnen hier überfragt bin - auf jedenfall
> chapeau, monsieur!

Ist sehr interessant, leider ist es auch nicht mein Thema. Ich bin 
soweit es Mechanik angeht, eher im FEM-Bereich fit. Da hat es ähnliche 
Fragestellungen.

Allerdings scheint mir das Thema hier 2-geteilt:

1) inhaltliche Fragen zum richtigen Ansatz

2) Unterschiede bei Octave.

Zu 2 kann ich sagen, dass Octace einfach vielfältiger ist, allerdings 
die E-Modelle schlecht unterstüzt sind.

Ich würde das mal im dortigen Forum aufwerfen!

[Michi S. (mista_s)]

Andreas B. schrieb:
> Ich möchte daher aus dem digitalisierten Output-Impuls
> den dazugehörigen Input-Impuls berechnen.

Ich möchte mal behaupten, das ist theoretisch wohl nicht möglich, denn 
jede reale Übertragungsfunktion ist zwangsweise verlustbehaftet, weil 
Bandbreitenbegrenzt; obendrein ist auch die Digitalisierung natürlich 
nochmals verlustbehaftet.

Damit enthält Dein digitalisierter Output-Impuls schlicht nicht 
ausreichende Informationen, um den Input-Impuls eindeutig zu bestimmen.

Ob das jetzt aber die wesentliche Ursache für Deine Probleme ist, da bin 
ich überfragt.

[Giovanni (sqrt_minus_eins)]

Folgende Erkenntnisse mit Ende 2024 (und gleichzeitig auch das Ende 
meiner Aktivitäten).

Auf Basis der ursprünglich übermittelten Daten habe ich versucht die 
Übertragungsfunktion aus der LTspice Simulation im Frequenzbereich zu 
ermitteln. Das Modell zeigt ein sehr gute Übereinstimmung (Daten in der 
MAT Datei).

Dann der Versuch die Funktion umzudrehen *) und aus dem Ausgangssignal - 
Punkt n7 das Eingangssignal zu rekonstruieren.
*) dazu notwendig Filter um die Ordnung Zähler/Nenner für Simulation 
richtig zu stellen.
Ein Versuch deshalb, weil ich das gezeigte Ergebnis durch Versuche mit 
verschiedenen Parameter gefunden habe (Bild: Versuch1.png).

Dies ist eine Aufgabe für Experten in DSP (digital signal processing) - 
der ich nicht bin.