Python 在信号处理中的优势

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

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

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