Error al identificar correctamente Ls y Rs de Vin y Vo medidos en LTSpice

Como ejercicio (usando Python para procesar los datos exportados de LTSpice), estoy tratando de determinar los valores de Rs y Ls que se usaron en un circuito, a partir de las medidas de Vin, Vout y el adelanto/retraso de fase entre los dos.

Simulé un circuito usando un valor específico para Rs y Ls, y luego traté de determinar cuáles eran esos valores a partir de los datos medidos. El circuito de muestra se muestra a continuación, así como las formas de onda obtenidas de él.

LTSpice no muestrea a intervalos uniformes, por lo que un truco para recuperar datos muestreados uniformemente es tomar una FFT de las señales y luego tomar otra FFT de las señales FFT para recuperar los datos del dominio del tiempo, ahora con muestreo uniforme (LTSpice no interpolación para la FFT). Hago esto y luego exporto estos datos a una secuencia de comandos de Python donde trato de determinar las amplitudes de Vo y Vi, y el retraso de fase entre ellos, y luego uso esos valores para identificar las R y las L a través de las siguientes fórmulas :

desde ,

entonces
donde V1 y Vo son las amplitudes de las sinusoides, y

Por lo tanto, resolviendo la magnitud y la fase de Z, puedo determinar los valores de R y L:

Desde

entonces

.

El problema es que la respuesta que obtengo es incorrecta, y dado que no hay "ruido" en las mediciones, supuse que obtendría fácilmente la respuesta exacta. Suponiendo que mis cálculos anteriores sean correctos, creo que el problema se debe al hecho de que parece que no puedo determinar con precisión las amplitudes de Vo, Vi y el retraso de fase entre ellos a partir de las muestras procesadas.

Intenté la diferenciación para encontrar los picos y el retraso de fase: diferencio las formas de onda y luego interpolo los valores de tiempo y voltaje cuando ocurre un cambio en el signo que indica el máximo (y el mínimo) de las formas de onda. Luego promedio los valores de los picos y la diferencia entre los picos de Vin y los picos de Vo sobre todas las muestras para tratar de reducir el error.

También probé la correlación cruzada circular para encontrar el retraso de fase, pensando que quizás el error en la determinación del retraso de fase era demasiado grande.

Finalmente, intenté ajustar las diferentes muestras de Vo y Vi a sinusoides utilizando el método de mínimos cuadrados.

Los tres casos no me dan la solución correcta para Ls y Rs, como se muestra en la tabla (también incluye la solución calculada ideal)

Puedo compartir los datos, el código del cuaderno de Python, etc., si alguien está interesado en ayudarme a descubrir por qué este ejercicio aparentemente sencillo no funciona para mí.

Gracias.

[EDITAR] Se agregó algo de código....

Configuración:

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]

Método Derivado:

# 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)

Método XCorr

desde numpy.fft importar 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

Estimador LS para encontrar Amplitud

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]

Solución FFT

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)

Resultados en:

Solución ideal

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)

Resumen de la solución Después de la recomendación de @ocspro acerca de limitar el paso máximo en la simulación LTSpice a 1n, los resultados son mejores, aunque no 100 % correctos en la identificación de las Rs. Quizás esto se deba a la alta sensibilidad numérica del cos (fase Z) alrededor del pi/2... no estoy seguro de cómo abordar esto (ver la solución xcorr). De todos modos, parece que todos los enfoques tomados arrojan soluciones similares (excepto para Xcorr donde Rs es negativo debido a que la fase Z calculada es ligeramente mayor que pi/2):

Método de diferenciación

Método Xcorr

Método de mínimos cuadrados

Método FFT

Solución ideal