The Mectron Load Line

In this notebook, we'll explore load line placement using extracted data from the 25L6G plate characteristic curves. The interactive example allows changing $I_a, V_a, n$ and $R_l $ (speaker impedance) while drawing and calculating output performance estimations.

The code then determines the transfer characteristic from the selected operating point and extracted data. It then calculates the circuit transient response highlighting second harmonic distortion. Finally, the user can play 3 normalized audio levels with varying levels of second harmonic to demonstrate what a distorted sine wave sounds like.


Mectron MPA-II Schematic


25L6GT - Tung-Sol 1957


Data Extraction


In [1]:
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import sys
from scipy import interpolate
import math
from ipywidgets import *
from scipy.fftpack import fft
from scipy.optimize import fsolve
from scipy import signal


# used engauge to extract plot data from datasheet
fn = "25L6GT - Triode.csv"
datasheetCurveData = pd.read_csv(fn)

colnames = datasheetCurveData.columns.values
colcount = len(colnames)
rowcount = len(datasheetCurveData[colnames[0]])

# do a little clean up
# remove negative values of plate current
for i in range(1,colcount):
    for j in range(rowcount):
        if datasheetCurveData[colnames[i]][j] < 0.00:
            datasheetCurveData[colnames[i]][j] = 0.0
            
# engauge adds bogus points on curves where plate voltage is greater than curve
# scan through array data column at a time, find point where engauge starts duplicating data
# calculate slope and y-axis intercept (m, b) then fill data past that point with a line
for i in range(1,colcount):
    m = 0
    b = 0
    for j in range(rowcount):
        if datasheetCurveData[colnames[i]][j] > 0.01:
            try:
                if m == 0:
                    if datasheetCurveData[colnames[i]][j] == datasheetCurveData[colnames[i]][j+2]:
                        m = (datasheetCurveData[colnames[i]][j] - datasheetCurveData[colnames[i]][j-5])/(datasheetCurveData['PlateVoltage'][j] - datasheetCurveData['PlateVoltage'][j-5])
                        b = datasheetCurveData[colnames[i]][j] - m*datasheetCurveData['PlateVoltage'][j]
                        # print j,datasheetCurveData[colnames[i]][j],m,b
                if m != 0:
                    datasheetCurveData[colnames[i]][j] = m*datasheetCurveData['PlateVoltage'][j] + b
            except KeyError:
                pass # j+2 is now > rowcount for the higher curves

In [2]:
datasheetCurveData.head(5)


Out[2]:
PlateVoltage 0V -5V -10V -15V -20V -25V -30V -35V -40V -45V -50V
0 0.30 0.00030 0.0 0.0 0.0 0.000004 0.0 0.0 0.000001 0.000007 0.0 5.893480e-07
1 0.89 0.00030 0.0 0.0 0.0 0.000004 0.0 0.0 0.000001 0.000007 0.0 0.000000e+00
2 1.48 0.00030 0.0 0.0 0.0 0.000004 0.0 0.0 0.000000 0.000007 0.0 0.000000e+00
3 2.08 0.00043 0.0 0.0 0.0 0.000004 0.0 0.0 0.000000 0.000007 0.0 0.000000e+00
4 2.67 0.00071 0.0 0.0 0.0 0.000004 0.0 0.0 0.000000 0.000007 0.0 0.000000e+00

In [3]:
datasheetCurveData.tail(5)


Out[3]:
PlateVoltage 0V -5V -10V -15V -20V -25V -30V -35V -40V -45V -50V
467 277.45 0.477849 0.311147 0.227536 0.186323 0.115034 0.084500 0.056021 0.038267 0.025745 0.009730 0.00346
468 278.04 0.479063 0.311994 0.228229 0.187020 0.115495 0.084931 0.056389 0.038641 0.026188 0.009970 0.00355
469 278.63 0.480276 0.312840 0.228922 0.187717 0.115956 0.085362 0.056756 0.039014 0.026630 0.010110 0.00365
470 279.23 0.481511 0.313700 0.229627 0.188426 0.116425 0.085801 0.057130 0.039394 0.027080 0.010329 0.00375
471 279.82 0.482725 0.314547 0.230321 0.189124 0.116886 0.086232 0.057498 0.039767 0.027523 0.010544 0.00385


Load Line Placement


In [4]:
#initial values
Ia = 0.0484  #plate current mA
Va = 118.3   #plate voltage V
Rl = 4       #speaker impedance
n  = 33      #pri/sec turns ratio

VaMAX = 300.0
IaMAX = 0.120
GraphWidth = 840
GraphHeight = 390

# later, we find intersection of loadline with plate current curves by resampling
# so all have the same x values.
# http://stackoverflow.com/questions/17928452/find-all-intersections-of-xy-data-point-graph-with-numpy
# http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html

PlateVoltages = np.arange(0,VaMAX,1.0)

saturationCurveVoltage = '0V'
cutoffCurveVoltage = '-40V'

# creating 1D interpolation functions from the datasheet extracted curves
iaf = {}
for i in range(1,colcount):
    iaf[colnames[i]] = {'valueAt': None,'loadLineIntersectionV':0,'loadLineIntersectionI':0}
    iaf[colnames[i]]['valueAt'] = interpolate.interp1d(datasheetCurveData['PlateVoltage'].tolist(), 
                                            datasheetCurveData[colnames[i]].tolist())

def plot(_Ia,_Va,_Rl,_n):
    global i0intersect, i16intersect,Ia,Va,Rl,n,b,m
    Ia = _Ia # set the slider values to global
    Va = _Va
    Rl = _Rl
    n  = _n

    # plot the csv colums versus plate/anode voltage
    fig = plt.figure(figsize=(15, 8))
    null = [plt.plot(datasheetCurveData['PlateVoltage'],
                     datasheetCurveData[colnames[x]],label='') for x in range(1,colcount)]
    plt.grid(linestyle='--', linewidth=0.5)
    null = plt.xticks(np.arange(0,VaMAX,10))
    null = plt.yticks(np.arange(0,IaMAX,0.01))

    # plot power dissipation limit curve
    Pd = 10.0 # 10W
    null = plt.plot(PlateVoltages, Pd/PlateVoltages,label='Power Limit',linestyle='--')

    null = plt.xlim([0,VaMAX])
    null = plt.ylim([0,IaMAX])

    def placeLabel(plt,text,x,y,angle):
        null = plt.annotate(s=text,
                        rotation=angle,
                        xy=(x,y),
                        xycoords='data',
                        xytext=(0,0),
                        textcoords='offset points',
                        bbox=dict(boxstyle="round", fc="1.0"),
                        size=12,
                        # arrowprops=dict(arrowstyle="->",connectionstyle="angle,angleA=0,angleB=70,rad=10")
                           )

    powerLimitVoltage = 100.0
    slope = -Pd/(powerLimitVoltage*powerLimitVoltage)
    graphSlope = slope*(GraphHeight/IaMAX)/(GraphWidth/VaMAX)
    angle = (180.0/np.pi) * np.arctan(graphSlope)
    placeLabel(plt,"Power\n%.1fW"%Pd,powerLimitVoltage,Pd/powerLimitVoltage,angle)

    for i in range(1,colcount):
        l = colnames[i]
        for j in range(rowcount):
            if datasheetCurveData['PlateVoltage'][j]*datasheetCurveData[colnames[i]][j] > (Pd+1.0):
                placeLabel(plt,l,datasheetCurveData['PlateVoltage'][j],datasheetCurveData[colnames[i]][j],0)
                break

    plateImpedance = float(Rl * n**2)

    m = -1/plateImpedance
    b = Ia + Va/plateImpedance
    ll = m*PlateVoltages+b

    null = plt.plot(PlateVoltages,ll,label='%d ohm Load'%Rl)
    null = plt.plot(Va,Ia, 'or',label='Op Point',color='g')

    for i in range(1,colcount):
        mindiff = 10
        for v in PlateVoltages:
            try:
                ia = iaf[colnames[i]]['valueAt'](v)
                iall = m*v+b
                diff = abs(ia - iall)

                if diff < mindiff:
                    vinter = v
                    iinter = iall
                    mindiff = diff
            except ValueError:
                pass

        iaf[colnames[i]]['loadLineIntersectionV'] = vinter
        iaf[colnames[i]]['loadLineIntersectionI'] = iinter

        if colnames[i] == cutoffCurveVoltage:
            break

    vsat = iaf[colnames[1]]['loadLineIntersectionV']
    isat = iaf[colnames[1]]['loadLineIntersectionI']
    null = plt.annotate(s="%.1fV@%.1fmA"%(vsat,isat*1000),
                        xy=(vsat,isat),
                        xycoords='data',
                        xytext=(-120,20),
                        textcoords='offset points',
                        bbox=dict(boxstyle="round", fc="1.0"),
                        size=12,
                        arrowprops=dict(arrowstyle="->",
                                        connectionstyle="angle,angleA=0,angleB=110,rad=10"))
    null = plt.annotate(s="%.1fV@%.1fmA"%(Va,Ia*1000),
                        xy=(Va,Ia),
                        xycoords='data',
                        xytext=(20,20),
                        textcoords='offset points',
                        bbox=dict(boxstyle="round", fc="1.0"),
                        size=12,
                        arrowprops=dict(arrowstyle="->",
                                        connectionstyle="angle,angleA=0,angleB=70,rad=10"))
    vcut = iaf[cutoffCurveVoltage]['loadLineIntersectionV']
    icut = iaf[cutoffCurveVoltage]['loadLineIntersectionI']
    null = plt.annotate(s="%.1fV@%.1fmA"%(vcut,icut*1000),
                        xy=(vcut,icut),
                        xycoords='data',
                        xytext=(20,20),
                        textcoords='offset points',
                        bbox=dict(boxstyle="round", fc="1.0"),
                        size=12,
                        arrowprops=dict(arrowstyle="->",
                                        connectionstyle="angle,angleA=0,angleB=70,rad=10"))


    # from RCA RC-22 P19
    distortion = abs((((icut+isat)/2 - Ia)/(icut-isat))*100)

    dvlower = Va - vsat
    dilower = isat - Ia
    dvhigher = vcut - Va
    dihigher = Ia - icut

    if dvlower < dvhigher:
        # closer to saturation
        Pout =  dvlower*(1/math.sqrt(2))*dilower
        title = "~%.2fW@%0.1f%% - closer to saturation\n%.1fV@%.1fmA, %.1fV@%.1fmA, %.1fV@%.1fmA, dV=%.1fV, dI=%.1fmA\n%d ohms - %d ohms"%(Pout,
                                                                            distortion,
                                                                            vsat,isat*1000,
                                                                            Va,Ia*1000,
                                                                            vcut,icut*1000,
                                                                            dvlower,dilower*1000,
                                                                            Rl,
                                                                            plateImpedance)
        null = plt.plot((vsat,Va),(isat,Ia),linewidth=2,color='b')
    else:
        # closer to cutoff
        Pout =  (vcut-Va)*(1/math.sqrt(2))*(Ia-icut)
        title = "~%.2fW@%0.1f%% - closer to cutoff\n%.1fV@%.1fmA, %.1fV@%.1fmA, %.1fV@%.1fmA, dV=%.1fV, dI=%.1fmA\n%d ohms - %d ohms"%(Pout,
                                                                        distortion,
                                                                        vsat,isat*1000,
                                                                        Va,Ia*1000,
                                                                        vcut,icut*1000,
                                                                        dvhigher,dihigher*1000,
                                                                        Rl,
                                                                        plateImpedance)
        null = plt.plot((Va,vcut),(Ia,icut),linewidth=2,color='b')



    for i in range(1,colcount):
        if iaf[colnames[i]]['loadLineIntersectionV']:
            vinter = iaf[colnames[i]]['loadLineIntersectionV']
            iinter = iaf[colnames[i]]['loadLineIntersectionI']
            null = plt.plot(vinter,iinter,'or',color='#EEEEEE')


    null = plt.suptitle(title,fontsize=14, fontweight='bold')
    null = plt.xlabel('Plate Voltage')
    null = plt.ylabel('Plate Current')
    null = plt.legend()

null = interact(plot,
             _Ia=widgets.FloatSlider(min=0.01,max=0.1,step=0.0025,value=Ia),
             _Va=widgets.FloatSlider(min=50,max=300,step=5,value=Va),
             _Rl=widgets.FloatSlider(min=2,max=16,step=2,value=4),
             _n=widgets.FloatSlider(min=10,max=40,step=1,value=33))


From RCA RC-22, p19, $$THD\% = \large{\frac{\frac{I_{max}+I_{min}}{2}-I_O}{I_{max}-I_{min}} * 100}$$
The THD shown above uses this equation, but its unrealistic because the circuit's operating point isn't near the load line center. Large vg will drive the circuit into saturation easily. Let's move on to show more detailed explanation.


Moving on to the transfer characteristic

We'll take the load line/curve intersections shown in graph above and plot those, $v_g$ versus $i_a$.


In [5]:
xlimit = -40.0

fig = plt.figure(figsize=(15, 10))

x = []
y = []
for i in range(1,colcount):
    if iaf[colnames[i]]['loadLineIntersectionV'] > 0:
        x.append(float(colnames[i][:-1]))
        y.append(iaf[colnames[i]]['loadLineIntersectionI'])

plt.plot(x,y,label="Va = %.0fV"%Va,marker='o')

plt.tick_params(axis='y', which='both', labelleft='off', labelright='on')
plt.grid(True)

'''
for i in range(len(x)):
    print "%.6f,%0.6f"%(x[i],y[i])

'''

coeff = np.polyfit(x, y, 2) # changing the '2' has implications for the rest of notebook and js

def func(x):
    global coeff,Ia
    return coeff[0]*x*x+coeff[1]*x+coeff[2] - Ia

title = "extracted transfer curve\npolynomial curve fit"

x = np.linspace(xlimit,0)
null = plt.plot(x,coeff[0]*x*x+coeff[1]*x+coeff[2],marker='x',label="Curve Fit",linestyle='None')

# using the curve-fitted curve, find Vg at Ia
Vg = fsolve(func,-8.0)

null = plt.annotate(s="%.1fV@%.1fmA"%(Vg,Ia*1000),
                        xy=(Vg,Ia),
                        xycoords='data',
                        xytext=(-120,20),
                        textcoords='offset points',
                        bbox=dict(boxstyle="round", fc="1.0"),
                        size=12,
                        arrowprops=dict(arrowstyle="->",
                                        connectionstyle="angle,angleA=0,angleB=100,rad=10"))

null = plt.plot(Vg,Ia, 'or',label='Op Point',color='g')
null = plt.suptitle(title,fontsize=14, fontweight='bold')
null = plt.xlim([xlimit,0])
null = plt.legend(loc='upper left')
null = plt.xlabel('Grid Voltage')
null = plt.ylabel('Plate Current')


Transient Response and Harmonic Content


In [6]:
ylowerlimit = -40.0

def plotTransient(vgamplitude):
    ra = Rl*n*n
    vaMax = -b/m

    N = 2000
    T = 1.0 / 5000.0
    f = 100.0

    t = np.linspace(0.0, N*T, N)
    vi = Vg + vgamplitude*np.sin(f*2.0*np.pi*t)
    ia = coeff[0]*vi*vi + coeff[1]*vi + coeff[2]
    vo = vaMax - ia*ra

    yf = fft(vo)
    xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

    plt.figure(figsize=(15, 10))

    ax = plt.subplot(311)
    null = plt.xlim([0,0.02])
    null = plt.ylim([ylowerlimit,0])
    ax.plot(t,vi)
    ax.grid(True)
    ax.tick_params(axis='x', which='both', labeltop='on',labelbottom='off')
    ax.set_xlabel("Grid Voltage")
    
    ax = plt.subplot(312)
    null = plt.xlim([0,0.02])
    null = plt.ylim([0,vaMax])
    ax.plot(t,vo)
    ax.grid(True)
    ax.set_xlabel("Plate Voltage")
   
    ax = plt.subplot(313)
    yfl = 2.0/N * np.abs(yf[:N/2])
    ax.plot(xf[1:], 20*np.log10(yfl[1:])) # take the dc component off
    null = plt.xlim([0,500])
    ax.grid(True)
    ax.set_xlabel("Frequency")
    ax.set_ylabel("dB")
     
    mags = []
    for i in range(len(xf)):
        if xf[i]%100 < 1:
            mags.append(2.0/N*np.abs(yf[i]))
            # print xf[i],2.0/N*np.abs(yf[i])

    thdsum = 0
    for i in range(2,8): # all but the fundamental
        thdsum += mags[i]*mags[i]/ra

    # print "sum",sum
    dist = 100*math.sqrt(thdsum/(mags[1]*mags[1]/ra)) # root mean square
    
    thdsum += mags[1]*mags[1]/ra # add power from fundamental 
    power = thdsum/math.sqrt(2)

    null = plt.suptitle("Anode Voltage AC\nvg amplitude = %.2fV, Vg=%.1fV\n%.2f%% %.2fW"%(vgamplitude,Vg,dist,power),fontsize=14, fontweight='bold')
    
    try:
        peakind = signal.find_peaks_cwt(20*np.log10(yfl), np.arange(1,10))
        # print peakind, xf[peakind], yfl[peakind]

        # the peaks are off by 1 for some reason
        i1 = peakind[0]-1
        p1 = yfl[i1]
        p1l = 20 * np.log10(p1)

        i2 = peakind[1]-1
        p2 = np.abs(yfl[i2])
        p2l = 20 * np.log10(p2)

        null = plt.annotate(s="%.1fdB, %.1f"%(p1l,p1),
                    xy=(xf[i1],p1l),
                    xycoords='data',
                    xytext=(30,0),
                    textcoords='offset points',
                    bbox=dict(boxstyle="round", fc="1.0"),
                    size=12,
                    arrowprops=dict(arrowstyle="->",
                                    connectionstyle="angle,angleA=0,angleB=10,rad=10"))
        null = plt.annotate(s="%.1fdB, %.1f"%(p2l,p2),
                    xy=(xf[i2],p2l),
                    xycoords='data',
                    xytext=(30,0),
                    textcoords='offset points',
                    bbox=dict(boxstyle="round", fc="1.0"),
                    size=12,
                    arrowprops=dict(arrowstyle="->",
                                    connectionstyle="angle,angleA=0,angleB=10,rad=10"))

        null = plt.annotate(s="peak diff %.1fdB"%(abs(p1l-p2l)),
                    xy=((xf[i2]+xf[i1])/2,p2l),
                    xycoords='data',
                    xytext=(-40,-30),
                    textcoords='offset points',
                    bbox=dict(boxstyle="round", fc="1.0"),
                    size=12)

    except IndexError:
        pass # no peaks

    plt.show()
    
null = interact(plotTransient, vgamplitude=widgets.FloatSlider(min=0,max=30,step=0.25,value=1))



In [7]:
from IPython.display import HTML

# references
#   https://jakevdp.github.io/blog/2013/06/01/ipython-notebook-javascript-python-communication/
#   https://github.com/mdn/audio-buffer/blob/gh-pages/index.html
    
javascript = """
<script>
    var audioCtx = new (window.AudioContext || window.webkitAudioContext)();
    var pre = document.querySelector('pre');
    var myScript = document.querySelector('script');

    var channels = 2;
    var frameCount = audioCtx.sampleRate * 5.0;

    var myArrayBuffer = audioCtx.createBuffer(channels, frameCount, audioCtx.sampleRate);
    var play = function(Vg,vgamplitude,Ia,a0,a1,a2,gain) {
      console.log(Vg,vgamplitude,Ia,a0,a1,a2)
      for (var channel = 0; channel < channels; channel++) {
       // This gives us the actual array that contains the data
       var nowBuffering = myArrayBuffer.getChannelData(channel);
       min = 100;
       max = -100;
       for (var k = 0; k < frameCount; k++) {
         // audio needs to be in [-1.0; 1.0]
         vi = Vg + vgamplitude*Math.sin(50*3.14*k);
         i = (a0*vi*vi+a1*vi+a2);
         val = i*gain;
         if (val < min) {
             min = val;
         }
         if (val > max) {
             max = val;
         }
         nowBuffering[k] = val;
       }
       ave = (max+min)/2;
       for (var k = 0; k < frameCount; k++) {
           nowBuffering[k] -= ave;
       }
       console.log(nowBuffering);
      }

      var source = audioCtx.createBufferSource();
      source.buffer = myArrayBuffer;
      source.connect(audioCtx.destination);
      source.start();
    }
  </script>
"""

HTML(javascript)


Out[7]:

In [8]:
# these 3 js buttons call the play function with Vg,vgamplitude,Ia,a0,a1,a2,gain where gain
# is set so fundamental has same amplitude, i.e. close to same loudness, so the user speaker
# amplifer doesn't distort and taint resulting audio.
# I use https://play.google.com/store/apps/details?id=com.zephyr.soundAnalyserPRO for spectrum
# display

# print Vg,Ia,coeff

input_form = """
<button onclick='play(%f,%f,%f,%f,%f,%f,%f)' style="width:200px">Play vg=0.1Vpp</button><br>
<button onclick='play(%f,%f,%f,%f,%f,%f,%f)' style="width:200px">Play vg=4pp</button><br>
<button onclick='play(%f,%f,%f,%f,%f,%f,%f)' style="width:200px">Play vg=8pp</button><br>
<button onclick='play(%f,%f,%f,%f,%f,%f,%f)' style="width:200px">Play vg=16pp</button><br>
<div id='buf'><div>
"""%(Vg,0.1,Ia,coeff[0],coeff[1],coeff[2],80,
     Vg,4,Ia,coeff[0],coeff[1],coeff[2],2,
     Vg,8,Ia,coeff[0],coeff[1],coeff[2],1.1,
     Vg,16,Ia,coeff[0],coeff[1],coeff[2],0.6
    )
HTML(input_form)


Out[8]:





In [ ]: