3
\$\begingroup\$

As an exercise (using Python to process LTSpice exported data), I am trying to determine the values of Rs and Ls that were used in a circuit, from the the measurements of Vin, Vout, and the phase lead/lag between the two.

I simulated a circuit using specific value for Rs and Ls,and then tried to determine what those values were from the measured data. The sample circuit is shown below, as well as the waveforms obtained from it.

circuit diagramwaveforms

LTSpice does not sample at uniform intervals, so a trick to recover uniformedly sampled data is to take an FFT of the signals, and then take another FFT of the FFT'd signals to recover the time domain data, now with uniform sampling (LTSpice does interpolation for the FFT). I do this, and then export this data to a Python script where I try to determine the amplitudes of Vo and Vi, and the phase delay,, between them, and then use those values to identify the Rs and Ls, via the following formulas:

Since vout,

then
enter image description here where V1 and Vo are the amplitudes of sinusoids, and enter image description here

Therefore, solving for the magnitude and phase of Z, I can determine the values of R and L:
enter image description here

enter image description here

Since enter image description here

then

eq10.

The problem is that the answer I get is wrong, and since there is no "noise" in the measurements, I assumed I would easily get the exact answer. Assuming my math above is correct, I think the problem stems from the fact that I can't seem to accurately determine the amplitudes of Vo, Vi, and the phase delay between them from the processed samples.

I've tried differentiation to find the peaks and phase delay: I differentiate the waveforms, and then interpolate the values for time and voltage when a change in sign occurs indicating the max (and mins) of the waveforms. I then average the peak values and the difference between the Vin peaks and the Vo peaks over all the samples to try and reduce the error.

I also tried circular crosscorrelation to find the phase delay, thinking perhaps the error in the phase delay determination was too great.

Finally, I tried fitting the different Vo and Vi samples to sinusoids using Least Squares method.

All three cases fail to give me the correct solution for Ls and Rs, as shown in the table (it includes the ideal, computed solution as well) enter image description here

I can share the data, the python notebook code, etc, if anyone is interested in helping me figure out why this seemingly straightforward exercise is not working for me.

Thank you.

[EDIT] Added some code....

Setup:

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]

Derivative Method:

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

waveform sketch

XCorr Method

from numpy.fft import 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

xcorr waveforms

LS Estimator to find 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)

enter image description here

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 Solution

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)

Results in:

enter image description here

fft_sol

Ideal Solution

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)

Solution Summary After recommendation from @ocspro about limiting the max step in the LTSpice simulation to 1n, the results are better, although not 100% correct in the identification of the Rs. Perhaps this is due to the high numerical sensitivity of the cos(phaseZ) around the pi/2...not sure how to address this (see xcorr solution). Anyway, it seems all approaches taken yield similar solutions (except for Xcorr where Rs is negative due to calculated phaseZ being slightly greater than pi/2):

Differentation Method diff_method_sol

Xcorr Method xcorr_method_sol

Least Square Method ls_sqr_sol

FFT Method FFT_method_sol

Ideal Solution

Ideal Soln

\$\endgroup\$
20
  • \$\begingroup\$ Step 1: Can you read the LTspice output file into your Python code and just re-plot it to verify you are reading it correctly. Step 2: Can you estimate the magnitude and phase of \$V_i\$ and \$V_o\$ correctly? Only after making sure you can do steps 1 and 2 should you worry about step 3: Use \$V_i\$ and \$V_o\$ to estimate the circuit variables \$R_s\$, \$L_s\$, and \$R_1\$. \$\endgroup\$
    – The Photon
    Mar 3, 2020 at 20:51
  • \$\begingroup\$ Also, does your Python code assume the frequency of the signal is known a priori, or do you also have to estimate the frequency from the data? \$\endgroup\$
    – The Photon
    Mar 3, 2020 at 20:52
  • \$\begingroup\$ How many point FFT? Have you experimented with different FFT lengths and seen any trend betweenFFT length and accuracy? \$\endgroup\$
    – user16324
    Mar 3, 2020 at 21:00
  • \$\begingroup\$ @ThePhoton, yes -- the estimates are all good...it is the accuracy that seems to be the issue. Frequency is known ahead of time. \$\endgroup\$
    – jrive
    Mar 3, 2020 at 21:29
  • 1
    \$\begingroup\$ "LTSpice does not sample at uniform intervals, so a trick to recover uniformedly sampled data is to take an FFT of the signals, and then take another FFT of the FFT'd signals to recover the time domain data" - that seems like a ridiculously complicated way to solve a simple problem. \$\endgroup\$ Mar 4, 2020 at 2:03

2 Answers 2

2
\$\begingroup\$

Maybe the accuracy error lies in the number of simulation points? Using the built-in fft-function in LTspice (.four 1000k 1 -1 V(Vi) V(Vo)), I got large discrepancies in magnitude and phase by running your simulation with standard simulation settings.

By decreasing the maximum time step, I managed to get magnitude and phase shift almost identical to the ideal solution

Comparison of max step of 1 ns vs none:

  • Vi = .9999 vs .9803
  • Vo = .7883 vs .7722
  • phase shift = .6625 vs .6475

Screenshot of circuit and log

Alternatively, one can decrease the relative tolerance to force shorter time steps. The automatic stepping is probably rather large due to the simple and linear circuit.

\$\endgroup\$
1
  • \$\begingroup\$ thank you. Much better. I still get an error in the identified Rs, but much closer than before: Rs_ident =0.204, Ls_ident =.204., using the Least squares method:Vi = [0.99996636] Vo = [0.78834297] Phase Dly = [0.6624794] Zmag = [1281.7480386] Zph = [1.5706335] Rs = [0.20870581] Ls = [0.000204]. Xcorr method doesn't work with the new data, I'm not sure why yet...the phase delay goes crazy. \$\endgroup\$
    – jrive
    Mar 4, 2020 at 14:41
1
\$\begingroup\$

In LTspice, a commonly forgotten aspect with using inductors is its default series resistance of 1 mΩ.
Make sure to explicitly set it to zero to get even better matching results!

enter image description here

\$\endgroup\$
2
  • \$\begingroup\$ Also add the bit about Rpar, which defaults to 1e12 times the inductance value. For low inductances, Rpar gets to be noticed, too. \$\endgroup\$ Mar 4, 2020 at 19:19
  • \$\begingroup\$ yes...I had done that.... thank you. \$\endgroup\$
    – jrive
    Mar 4, 2020 at 21:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.