Fehler bei der korrekten Identifizierung von Ls und Rs aus gemessenem Vin und Vo in LTSpice

Als Übung (unter Verwendung von Python zur Verarbeitung von LTSpice-exportierten Daten) versuche ich, die Werte von Rs und Ls zu bestimmen, die in einer Schaltung verwendet wurden, aus den Messungen von Vin, Vout und der Phasenvoreilung/-verzögerung zwischen den beiden.

Ich habe eine Schaltung mit bestimmten Werten für Rs und Ls simuliert und dann versucht, diese Werte anhand der gemessenen Daten zu bestimmen. Die Beispielschaltung ist unten gezeigt, ebenso wie die daraus erhaltenen Wellenformen.

LTSpice tastet nicht in gleichmäßigen Intervallen ab, daher besteht ein Trick, um gleichmäßig abgetastete Daten wiederherzustellen, darin, eine FFT der Signale zu nehmen und dann eine weitere FFT der FFT-Signale zu nehmen, um die Zeitbereichsdaten wiederherzustellen, jetzt mit gleichmäßiger Abtastung (LTSpice tut dies Interpolation für die FFT). Ich mache das und exportiere diese Daten dann in ein Python-Skript, wo ich versuche, die Amplituden von Vo und Vi und die Phasenverzögerung zwischen ihnen zu bestimmen, und dann diese Werte verwende, um die Rs und Ls über die folgenden Formeln zu identifizieren :

Seit ,

dann
wobei V1 und Vo die Amplituden von Sinuskurven sind, und

Daher kann ich durch Auflösen nach Größe und Phase von Z die Werte von R und L bestimmen:

Seit

Dann

.

Das Problem ist, dass die Antwort, die ich bekomme, falsch ist, und da es kein "Rauschen" in den Messungen gibt, nahm ich an, dass ich leicht die genaue Antwort bekommen würde. Unter der Annahme, dass meine obige Mathematik korrekt ist, rührt das Problem meiner Meinung nach von der Tatsache her, dass ich die Amplituden von Vo, Vi und die Phasenverzögerung zwischen ihnen aus den verarbeiteten Samples nicht genau bestimmen kann.

Ich habe versucht, die Spitzen und die Phasenverzögerung zu differenzieren: Ich differenziere die Wellenformen und interpoliere dann die Werte für Zeit und Spannung, wenn eine Vorzeichenänderung auftritt, die das Maximum (und das Minimum) der Wellenformen anzeigt. Ich mittele dann die Spitzenwerte und die Differenz zwischen den Vin-Spitzen und den Vo-Spitzen über alle Proben, um zu versuchen, den Fehler zu reduzieren.

Ich habe auch die zirkuläre Kreuzkorrelation versucht, um die Phasenverzögerung zu finden, weil ich dachte, dass der Fehler bei der Bestimmung der Phasenverzögerung vielleicht zu groß sei.

Schließlich habe ich versucht, die verschiedenen Vo- und Vi-Samples mit der Methode der kleinsten Quadrate an Sinuskurven anzupassen.

Alle drei Fälle geben mir nicht die richtige Lösung für Ls und Rs, wie in der Tabelle gezeigt (sie enthält auch die ideale, berechnete Lösung).

Ich kann die Daten, den Python-Notebook-Code usw. teilen, wenn jemand daran interessiert ist, mir dabei zu helfen, herauszufinden, warum diese scheinbar einfache Übung für mich nicht funktioniert.

Danke schön.

[BEARBEITEN] Etwas Code hinzugefügt ....

Aufstellen:

vi=data[1:,1]
vo=data[1:,2]
time=data[1:,0]
plt.plot(time,vi)
plt.plot(time,vo)
plt.show()

 deltaT= time[2]-time[1]

Ableitungsmethode:

# Use the derivatve to find the slope change from positive to negative. Select the point half way between previous and current sample for the corresponding values of vo, vi, and time.

# In[28]:


tmpderVo = 0
tmpderVi = 0

peakVotimearray=[]
peakVoarray    =[]
peakVitimearray=[]
peakViarray    =[]
derVoArr       =[]
derViArr       =[]

for i in range(len(time)-1):
    derVo = (vo[i+1]-vo[i])/deltaT  # derivative
    derVoArr.append(derVo)
    if np.sign(tmpderVo)==1:
        if np.sign(derVo) != np.sign(tmpderVo):
            # interpolate time and Vo 
            newtime= time[i]+deltaT/2
            peakVotimearray.append(newtime)
            peakVo = vo[i]+(vo[i+1]-vo[i])/2
            peakVoarray.append(peakVo)

    derVi=(vi[i+1]-vi[i])/deltaT  # derivative
    derViArr.append(derVi)
    if np.sign(tmpderVi)==1:
        if  np.sign(derVi) != np.sign(tmpderVi):
            # interpolate time and Vi  -- half way point 
            newtime= time[i]+deltaT/2
            peakVitimearray.append(newtime)
            peakVi = vi[i]+(vi[i+1]-vi[i])/2
            peakViarray.append(peakVi)

    tmpderVo = derVo
    tmpderVi = derVi

plt.plot(derVoArr[10:100],'-*k')
plt.show()
# Average Vo and Vi peaks
peakVoave= np.mean(np.array(peakVoarray)[10:-10])
stdVoave = np.std(np.array(peakVoarray)[10:-10])
peakViave= np.mean(np.array(peakViarray)[10:-10])
stdViave = np.std(np.array(peakViarray)[10:-10])

# Average Time delay
timedlyarray=np.array(peakVitimearray)-np.array(peakVotimearray)
timedlyave = np.mean(timedlyarray[10:-10])
timedlystd  = np.std(timedlyarray[10:-10])

print('time delay average= ', timedlyave, '\t Time Delay STD = ',timedlystd )
print('Coefficient of Variability of Delay Measurement = ', timedlystd/timedlyave)
print('\nAverage Vo Amplitude = ' , peakVoave,'\t Average Vi Amplitude = ' , peakViave)
print('Vo Amplitude STD = ', stdVoave, '\t Vi Amplitude STD = ', stdViave)
print('\nCoefficient of Variability of Vo Peak Measurement = ', stdVoave/peakVoave)
print('Coefficient of Variability of Vi Peak Measurement = ', stdViave/peakViave)

print('\nSkipped the first 10 values in array for average calculation to avoid any data edge effects\n')


plt.plot(time[periodstart:periodend],vo[periodstart:periodend], time[periodstart:periodend],vi[periodstart:periodend])
# indices for peak arrays no longer match original data indices, need to recover them below
frac = periodstart/len(time)  # what fraction of whole time array are we starting at
offset=int(len(peakVotimearray)*frac) # determine offset into peaktime array, one peak per period
plt.vlines(peakVotimearray[offset:int(offset+len(time[periodstart:periodend])*deltaT/1e-6)], -1, 1, colors='r', linestyles='dashed',linewidth=1)
plt.vlines(peakVitimearray[offset:int(offset+len(time[periodstart:periodend])*deltaT/1e-6)], -1, 1, colors='r', linestyles='dashed',linewidth=1)
plt.title('Sketch of Vi and Vo and their phase relationship')
plt.legend(['Vin','Vo'])
plt.show()

# ### Determine Time Delay using Peaks found via Derivatives
peakdly=timedlyave
peakdlyrad=timedlyave/T*2*np.pi
print(peakdlyrad)

XCorr-Methode

aus numpy.fft importiere fft, ifft

def periodic_corr(x, y):
    """Periodic correlation, implemented using the FFT.

    x and y must be real sequences with the same length.
    """
    return ifft(fft(x) * fft(y).conj()).real

xcorr=periodic_corr(vi,vo)
dlyndx= np.argmax(xcorr)
print('Index: ',dlyndx, '\tTime delay: ',dlyndx*deltaT)

plt.plot(time,vi)
plt.plot(time+dlyndx*deltaT,vo)
plt.show()
timedly=dlyndx*deltaT/T*2*np.pi

LS-Schätzer zum Ermitteln der Amplitude

D0 =np.array([np.cos(2*np.pi*f*time),np.sin(2*np.pi*f*time),np.ones(time.size)],'float').transpose()

vin=np.array([vi]).T
vout=np.array([vo]).T
print(np.concatenate([vin,vout], axis=1))

from numpy.linalg import inv

s_hat_vin = np.matmul(inv(np.matmul(D0.T,D0)),np.matmul(D0.T,vin))
s_hat_vo  = np.matmul(inv(np.matmul(D0.T,D0)),np.matmul(D0.T,vout))

vinMag = np.sqrt(s_hat_vin[0]**2+s_hat_vin[1]**2)
vinPh  = -np.arctan2(s_hat_vin[1],s_hat_vin[0])

voMag = np.sqrt(s_hat_vo[0]**2+s_hat_vo[1]**2)
voPh  = -np.arctan2(s_hat_vo[1],s_hat_vo[0])

print(vinMag,vinPh)
print(voMag, voPh)

#plt.plot(time,vo)
plt.plot(time,vinMag*np.cos(2*np.pi*f*time + vinPh) + s_hat_vin[2])
plt.plot(time,voMag *np.cos(2*np.pi*f*time + voPh) + s_hat_vo[2])
plt.plot(time,voMag *np.cos(2*np.pi*f*time + voPh+(vinPh-voPh)) + s_hat_vo[2])
plt.show()

lsm_dly = (vinPh-voPh)

def Zmagphase(vipeak,vopeak,theta,R1):
    """Returns the magnitude and phase of Z."""

    magZ  = (vopeak*R1)/(np.sqrt(vipeak**2 - 2*vipeak*vopeak*np.cos(theta)+ vopeak**2))
    phaseZ = theta - np.arctan2(-vopeak*np.sin(theta),(vipeak-vopeak*np.cos(theta)))

    return [magZ,phaseZ]

def Z2LsRs(mag,ph,f):
    """Determines Ls and Rs from Z in polar form"""
    w= 2*np.pi*f
    Rs = mag*np.cos(ph)
    Ls = mag*np.sin(ph)/(w)

    return [Rs, Ls]

FFT-Lösung

Fs = 1/deltaT
T= deltaT
N =  len(vi)
freq = Fs/N*np.arange(1,int(N/2)+1)

y=np.fft.fft(vi)
vimagfft=2/N*np.abs(y)[0:int(N/2)+1]
vimagfft=vimagfft[1:]
viphase = np.angle(y)[1:int(N/2)+1]

x=np.fft.fft(vo)
vomagfft=2/N*np.abs(x)[0:int(N/2)+1]
vomagfft=vomagfft[1:]
vophase = np.angle(x)[1:int(N/2)+1]

plt.plot(freq,vimagfft,'-*k', freq, vomagfft, '-*r')
plt.axis([0, 10000000, 0, 10])
plt.autoscale(True, 'y')
plt.show()


viFFT = np.max(vimagfft)
voFFT = np.max(vomagfft)

thetaFFT = vophase[np.argmax(vomagfft)]-viphase[np.argmax(vimagfft)]
print('ViampFFT = ', viFFT, '\t VoampFFT = ' , voFFT)
print('Phase Delay =',thetaFFT)

Ergebnisse in:

Ideale Lösung

from numpy import exp, abs, angle

def polar2z(r,theta):
    return r * exp( 1j * theta )

def z2polar(a,b):
    z = a + 1j * b
    return ( abs(z), angle(z) )

Vin = 1*exp(0j)


Vo=1*exp(0j)*(.118+(2*np.pi*1e6*204e-6j))/(1e3+.118+(2*np.pi*1e6*204e-6j))

Vomag=np.abs(Vo)
Votheta=np.angle(Vo)

magZideal= (Vomag*R1)/(np.sqrt(abs(Vin)**2 - 2*abs(Vin)*Vomag*np.cos(Votheta)+ Vomag**2))
print('Z_magIdeal = ', magZideal)


phZideal = Votheta - np.arctan2(-Vomag*np.sin(Votheta),(abs(Vin)-Vomag*np.cos(Votheta)))
print(phZideal)


R = magZideal*np.cos(phZideal)
L = magZideal*np.sin(phZideal)/(w)
print('R = ',R,'\t', 'L = ',L)

Zusammenfassung der Lösung Nach der Empfehlung von @ocspro, den maximalen Schritt in der LTSpice-Simulation auf 1n zu begrenzen, sind die Ergebnisse besser, wenn auch nicht 100 % korrekt bei der Identifizierung der Rs. Vielleicht liegt dies an der hohen numerischen Empfindlichkeit des cos (phaseZ) um den Pi / 2 ... nicht sicher, wie dies angegangen werden soll (siehe xcorr-Lösung). Wie auch immer, es scheint, dass alle gewählten Ansätze ähnliche Lösungen ergeben (mit Ausnahme von Xcorr, wo Rs negativ ist, da die berechnete phaseZ etwas größer als pi / 2 ist):

Differenzierungsmethode

Xcorr-Methode

Methode der kleinsten Quadrate

FFT-Methode

Ideale Lösung