In [29]:
from IPython.display import display
from IPython.display import HTML
from IPython.display import Image
import IPython.core.display as di # Example: di.display_html('<h3>%s:</h3>' % str, raw=True)
import warnings
warnings.filterwarnings('ignore')
# This line will hide code by default when the notebook is exported as HTML
di.display_html('<script>jQuery(function() {if (jQuery("body.notebook_app").length == 0) { jQuery(".input_area").toggle(); jQuery(".prompt").toggle();}});</script>', raw=True)
# This line will add a button to toggle visibility of code blocks, for use with the HTML export version
di.display_html('''<button onclick="jQuery('.input_area').toggle(); jQuery('.prompt').toggle();">Toggle code alt-t</button>''', raw=True)
di.display_html('''<script>document.addEventListener("keyup",function (e) {if (e.keyCode=84 && e.altKey) {jQuery(".input_area").toggle();}})</script>''',raw=True)
display(HTML("<style>.container { width:80% !important; }</style>"))
In [30]:
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import scipy
import scipy.signal as signal
import sys
from scipy import interpolate
import math
from ipywidgets import *
from scipy.fftpack import fft
from scipy.optimize import fsolve
# used engauge to extract plot data from datasheet
fn = "6f6g.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
display(HTML("<table><tr><td align='center'><img src='schematic.png'><br><a href='http://www.valvewizard.co.uk/se.html' target=_blank>Valve Wizard's The Single Ended Output Stage</a></td><td align='center'><img src='specs.png'><br><a href='https://frank.pocnet.net/sheets/127/6/6F6.pdf' target=_blank>Frank's Electron Tube Data</a></td></tr><tr><td colspan=2 align='center'><img src='6f6g.engauge.png'></a><br><a href='https://markummitchell.github.io/engauge-digitizer/' target=_blank>Engauge Digitizer</a></td></tr></table>"))
In [31]:
datasheetCurveData.head(5)
Out[31]:
In [32]:
datasheetCurveData.tail(5)
Out[32]:
In [33]:
#initial values
PSVoltage = 418.0
PSResistance = 3000.0
Ropt = 760.0
Ia = 0.034 #plate current mA
ScreenRatio = 5.33 #average ratio of plate/screen current
Va = 300.0 #anode voltage V
Rl = 4.0 #speaker impedance
n = 33.0 #pri/sec turns ratio
Pd = 11.0 #plate dissipation max
Screen = 250.0 #datasheet plate characteristic screen voltage
Vc = -16
VaMAX = 600.01
IaMAX = 0.1001
GraphWidth = 840 # get these from jpg size
GraphHeight = 560
# 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 = '-35'
scaledCurveData = []
# creating 1D interpolation functions from the datasheet extracted curves
iaf = {}
def plot(Plate_I,Screen_V,Speaker_R,Turns_Ratio,PSVoltage,PSResistance):
global i0intersect, i16intersect,Ia,Va,Vs,n,b,m,Pd,ScreenRatio,Vc,Ropt
global gridvoltages,platevoltages
global colcount, datasheetCurveData, scaledCurveData
Ia = Plate_I
Is = Plate_I / ScreenRatio
Vb = PSVoltage - PSResistance * (Ia + Is)
Va = Vb - Ropt*Ia
Vs = Screen_V
Rl = Speaker_R
n = Turns_Ratio
scaledCurveData = pd.DataFrame.copy(datasheetCurveData)
for x in range(1,colcount):
scaledCurveData[colnames[x]] = datasheetCurveData[colnames[x]]*(Vs/Screen)
for i in range(1,colcount):
iaf[colnames[i]] = {'valueAt': None,'loadLineIntersectionV':0,'loadLineIntersectionI':0}
iaf[colnames[i]]['valueAt'] = interpolate.interp1d(scaledCurveData['PlateVoltage'].tolist(), scaledCurveData[colnames[i]].tolist())
# plot the csv colums versus plate/anode voltage
fig = plt.figure(figsize=(15, 10))
for x in range(1,colcount):
null = plt.plot(scaledCurveData['PlateVoltage'],scaledCurveData[colnames[x]],label='')
plt.grid(linestyle='--', linewidth=0.5)
null = plt.xticks(np.arange(0,VaMAX,20))
null = plt.yticks(np.arange(0,IaMAX,0.01))
# plot power dissipation limit curve
null = plt.plot(PlateVoltages, Pd/PlateVoltages,label='Power Limit',linestyle='--')
null = plt.xlim([0,VaMAX])
null = plt.ylim([0,IaMAX])
null = plt.ylabel("Plate Current in amps")
null = plt.xlabel("Plate Voltage in volts")
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 = 120.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):
x = scaledCurveData['PlateVoltage'][j]
if scaledCurveData['PlateVoltage'][j]*scaledCurveData[colnames[i]][j] > (Pd+1.0):
placeLabel(plt,l,x,scaledCurveData[colnames[i]][j],0)
break
else:
if j == rowcount - 1:
placeLabel(plt,l,x - 15,scaledCurveData[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='r')
gridvoltages = []
platevoltages = []
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
gridvoltages.append(float(colnames[i]))
platevoltages.append(vinter)
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=(-30,120),
textcoords='offset points',
bbox=dict(boxstyle="round", fc="1.0"),
size=12,
arrowprops=dict(arrowstyle="->",
connectionstyle="angle,angleA=0,angleB=110,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=(-120,-10),
textcoords='offset points',
bbox=dict(boxstyle="round", fc="1.0"),
size=12,
arrowprops=dict(arrowstyle="->",
connectionstyle="angle,angleA=0,angleB=50,rad=10"))
# from RCA RC-22 P19
distortion = abs((((icut+isat)/2 - Ia)/(icut-isat))*100)
# figuring out the bias voltage and cathode resistor
row = np.where(scaledCurveData['PlateVoltage'] > Va)[0][0]
# print row, datasheetCurveData['PlateVoltage'][row],datasheetCurveData['PlateVoltage'][row-1]
for col in range(1,colcount):
x = scaledCurveData[colnames[col]][row]
if x < Ia:
x1 = x
y1 = float(colnames[col])
x2 = scaledCurveData[colnames[col-1]][row]
y2 = float(colnames[col-1])
# print x1,y1,x2,y2
m = (y2 - y1)/(x2 - x1)
b = y1 - m*x1
Vc = m*Ia + b
break
Rc = abs(Vc/(Ia+Is))
null = plt.annotate(s="%.1fV@%.1fmA\n%.1fV"%(Va,Ia*1000,Vc),
xy=(Va,Ia),
xycoords='data',
xytext=(-120,-30),
textcoords='offset points',
bbox=dict(boxstyle="round", fc="1.0"),
size=12,
arrowprops=dict(arrowstyle="->",
connectionstyle="angle,angleA=0,angleB=70,rad=10"))
dvlower = Va - vsat
dilower = isat - Ia
dvhigher = vcut - Va
dihigher = Ia - icut
Rs = (Vb-(Vs+(-Vc)))/Is
Prs = Rs*Is*Is
if dvlower < dvhigher:
# closer to saturation
Pout = dvlower*(1/math.sqrt(2))*dilower
title = "~%.2fW@%.1f%% - closer to saturation | %dohms => %dohms\nVb=%.1fV Vp=%.1fV Vs=%.1f/%.1fV Vg=%.1fV | Ip=%.1fmA Is=%.1fmA Ic=%.1fmA\nRc=%dohms Rs=%dohms@%.1fW\nDissipation - Plate=%.1fW Screen=%.1fW"%(
Pout,
distortion,
Rl,
plateImpedance,
Vb,Va,Vs,Vs+(-Vc),Vc,Ia*1000,Is*1000,1000*(Ia+Is),Rc,Rs,Prs,Va*Ia,Vs*Is)
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@%.1f%% - closer to cutoff | %dohms => %dohms\nVb=%.1fV Vp=%.1fV Vs=%.1f/%.1fV Vg=%.1fV | Ip=%.1fmA Is=%.1fmA Ic=%.1fmA\nRc=%dohms Rs=%dohms@%.1fW\nDissipation - Plate=%.1fW Screen=%.1fW"%(
Pout,
distortion,
Rl,
plateImpedance,
Vb,Va,Vs,Vs+(-Vc),Vc,Ia*1000,Is*1000,1000*(Ia+Is),Rc,Rs,Prs,Va*Ia,Vs*Is)
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.legend(loc='upper right')
# https://ipywidgets.readthedocs.io/en/latest/examples/Using%20Interact.html
null = interact(plot,
Plate_I = widgets.FloatSlider(min=0.001,max=0.050,step=0.001,value=Ia,continuous_update=False,readout_format='.3f'),
Screen_V = widgets.FloatSlider(min=50,max=400,step=5.0,value=Screen,continuous_update=False,readout_format='.0f'),
Speaker_R = widgets.FloatSlider(min=2,max=16,step=1,value=4,continuous_update=False,readout_format='.1f'),
Turns_Ratio = widgets.FloatSlider(min=10,max=50,step=1,value=46,continuous_update=False,readout_format='.1f'),
PSVoltage = widgets.FloatSlider(min=10,max=800,step=2,value=PSVoltage,continuous_update=False,readout_format='.0f'),
PSResistance = widgets.FloatSlider(min=10,max=5000,step=100,value=PSResistance,continuous_update=False,readout_format='.0f')
)
di.display_html('''<button onclick="jQuery('.input_area').toggle(); jQuery('.prompt').toggle();">Toggle code alt-t</button>''', raw=True)
In [34]:
def plotAnalysis(Vdrive):
global coeff,platevfit,amplitude,gridvoltages,platevoltages,scaledCurveData,Vc
coeff = np.polyfit(gridvoltages,platevoltages,4)
platevfit = np.poly1d(coeff)
biasvoltage = Vc
amplitude = Vdrive
'''
# Number of samplepoints
N = 16384*2
# sample spacingScreen
T = 1.0 / 500.0
'''
Fs = 1024 # sampling rate
Ts = 1.0/Fs # sampling interval
t = np.arange(0,1,Ts) # time vector
freq = 10.0
gridv = (amplitude/2)*np.sin(2.0*np.pi*freq*t)+biasvoltage
platev = platevfit(gridv)
# platev = np.random.normal(platev,0.001)
n = len(platev)
k = np.arange(n)
T = n/Fs
frq = k/T # two sides frequency range
frq = frq[range(int(n/2))]
yf = np.fft.fft(platev)/(n/2)
yf = yf[range(int(n/2))]
logyf = 20*np.log10(abs(yf))
peakfrq = np.where(logyf > -100)
# print np.where(t > -100)[0]
fig = plt.figure(figsize=(15, 15))
gainplot = fig.add_subplot(311)
x = np.linspace(0,-35)
y = platevfit(x)
null = gainplot.plot(x,y)
gainplot.set_ylabel("Plate Voltage V")
gainplot.set_xlabel("Grid Voltage V")
try:
x1 = np.where(x < (biasvoltage-(amplitude/2)))[0][0]
except IndexError:
x1 = len(x)
try:
x2 = np.where(x < (biasvoltage+(amplitude/2)))[0][0]
except IndexError:
x2 = 0
x21 = np.arange(x2,x1)
gainplot.plot(x[x21],y[x21],color='tab:green',linewidth=6)
null = gainplot.plot(biasvoltage,platevfit(biasvoltage), 'or',label='Op Point',color='r')
platevoltageplot = fig.add_subplot(312)
null = platevoltageplot.plot(t,platev,label='Plate Voltage')
null = platevoltageplot.set_xlim(0,4/freq)
null = platevoltageplot.set_ylim(0,600.0)
null = platevoltageplot.legend(loc='upper left')
null = platevoltageplot.set_ylabel("Plate Voltage V")
null = platevoltageplot.set_xlabel("Time")
gridvoltageplot = platevoltageplot.twinx()
null = gridvoltageplot.plot(t,gridv,color='tab:green',label='Grid Voltage')
null = gridvoltageplot.set_ylim(-35,0)
null = gridvoltageplot.legend(loc='upper right')
null = gridvoltageplot.set_ylabel("Grid Voltage V")
fftplot = fig.add_subplot(313)
null = fftplot.plot(frq, logyf)
#null = fftplot.set_xlim([0,6])
null = fftplot.set_ylim(-60,60)
null = fftplot.set_xlim(0,60)
null = fftplot.grid(axis='y')
fundfrq = peakfrq[0][1]
fundmag = logyf[fundfrq]
null = fftplot.annotate(s="fund",
xy=(fundfrq,fundmag),
xycoords='data',
xytext=(10,-5),
textcoords='offset points',
bbox=dict(boxstyle="round", fc="1.0"),
size=12)
null = fftplot.set_ylabel("Mag dBV")
null = fftplot.set_xlabel("Freq")
for i in range(2,len(peakfrq[0])):
f = peakfrq[0][i]
m = logyf[f]
if i == 2:
mag2 = m
label = "%+.1fdBV"%(m-fundmag)
else:
label = "%+.1fdBV\n2nd %+.1fdBV"%(m-fundmag,m-mag2)
null = fftplot.annotate(s=label,
xy=(f,m),
xycoords='data',
xytext=(10,-5),
textcoords='offset points',
bbox=dict(boxstyle="round", fc="1.0"),
size=12)
null = interact(plotAnalysis,
Vdrive =widgets.FloatSlider(min=1,max=40,step=0.5,value=10,continuous_update=False,readout_format='.1f'))
In [ ]: