import matplotlib.pyplot as plt
import numpy
import scipy.integrate
t = numpy.arange(0,4,0.001)
# duty cycle on phase A and B
Da = 0.70
Db = 0.40
def extendrange(ra,rb):
if ra is None:
return rb
elif rb is None:
return ra
else:
return (min(ra[0],rb[0]),max(ra[1],rb[1]))
def createLimits(margin, *args):
r = None
for x in args:
r = extendrange(r, (numpy.min(x),numpy.max(x)))
rmargin = (r[1]-r[0])*margin/2.0
return (r[0]-rmargin,r[1]+rmargin)
def showripple(centeralign=False):
# voltage waveforms on phases A and B
if centeralign:
sawtooth = abs(2*(t % 1) - 1)
Va = sawtooth < Da
Vb = sawtooth < Db
else:
ramp = t % 1
Va = ramp < Da
Vb = ramp < Db
Vab = Va - Vb
def ripple(x,t):
T = t[-1]-t[0]
meanval = numpy.mean(x)
# cumtrapz produces a vector of length N-1
# so we need to add one element back in
return numpy.append([0],scipy.integrate.cumtrapz(x - meanval,t))
Iab = ripple(Vab, t)
# plot results
margin = 0.1
fig = plt.figure(figsize=(8, 6), dpi=80)
ax = fig.add_subplot(3,1,1)
y = [Va*0.8, Vb*0.8+1]
ax.plot(t,y[0],t,y[1])
ax.set_yticks([0.4,1.4])
ax.set_yticklabels(['A','B'])
ax.set_ylim(createLimits(margin,y[0],y[1]))
ax.set_ylabel('Phase duty cycles')
ax = fig.add_subplot(3,1,2)
ax.plot(t,Vab)
ax.set_ylim(createLimits(margin,Vab))
ax.set_ylabel('Load voltage')
ax = fig.add_subplot(3,1,3)
ax.plot(t,Iab)
ax.set_ylim(createLimits(margin,Iab))
ax.set_ylabel('Ripple current')
savefile = 'pwm-%s-1.png' % ('center' if centeralign else 'edge')
fig.savefig(savefile, dpi=fig.dpi)
showripple(centeralign=False)
showripple(centeralign=True)
plt.show()
或者比较两个2级RC滤波器,一个具有相同RC并且一个具有第二级阻抗的滤波器增加10以减少负载(注意:下面的示意图不是用Python画的,而是在CircuitLab中手动画的)。
相应的Python代码片段为:
import matplotlib.pyplot as plt
import numpy
import itertools
# array version of the zip() function
def azip(*args):
iters = [iter(arg) for arg in args]
for i in itertools.count():
yield tuple([it.next() for it in iters])
# special case for 2 args
def azip2(a1,a2):
it1 = iter(a1)
it2 = iter(a2)
for i in itertools.count():
yield (it1.next(), it2.next())
def rcfilt(t,Vin,R,C):
N = len(C)
Vc = [0]*N
tprev = None
for (tj,Vj) in azip2(t,Vin):
if tprev is not None:
I = [(Vj-Vc[0])/R[0]] + [(Vc[k-1]-Vc[k])/R[k] for k in range(1,N)] + [0]
dt = tj - tprev
for k in range(N):
Vc[k] += (I[k]-I[k+1])/C[k]*dt
tprev = tj
yield numpy.array(Vc)
# 0-100 microseconds
t = numpy.arange(0,100,0.1)*1e-6
tus = t*1e6
Vin = (tus >= 10) * 1.0
# R1 = 1kohm, C1 = 10nF
# R2 = 10kohm, C2 = 1nF
R = [1000, 10000]
C = [10e-9, 1e-9]
Vc_a = numpy.array(list(rcfilt(t,Vin,R,C)))
R = [1000, 1000]
C = [10e-9, 10e-9]
Vc_b = numpy.array(list(rcfilt(t,Vin,R,C)))
fig = plt.figure(figsize=[8,6], dpi=80)
ylabels = ['Vc_a', 'Vc_b']
for (k,Vc) in enumerate([Vc_a,Vc_b]):
ax = fig.add_subplot(3,1,k+1)
ax.plot(tus,Vin,tus,Vc)
ax.legend(['Vin','Vc1','Vc2'])
ax.set_ylabel(ylabels[k])
ax.grid('on')
ax = fig.add_subplot(3,1,3)
ax.plot(tus,Vc_a[:,-1],tus,Vc_b[:,-1])
ax.legend(['Vc2_a','Vc2_b'])
ax.set_ylabel('Vc2')
ax.grid('on')
fig.suptitle('2-pole RC filters: Vc_a = 1K,10nF,10K,1nF; Vc_b = 1K,10nF,1K,10nF')
fig.savefig('rcfilt1.png',dpi=fig.dpi)
plt.show()
或者使用Python的sympy符号代数包来计算分段线性代码段的均方值:
from sympy import *
x0,x1,y0,y1,m,h = symbols('x0 x1 y0 y1 m h')
simplify(integrate((m*(x-x0)+y0)**2,(x,x0,x0+h)).subs(m,(y1-y0)/h))
你甚至可以在 SymPy Live server上自己亲自试试: