[翻译] Python 在信号处理中的优势

休息了几天回来了

前言

本篇是对Pylab的小试牛刀,也是对许多其他主题的过渡——包括《编码速度估计的长时间等待的后果》。

在工作中,我们使用 MATLAB 作为数据分析和可视化软件。但是在我的组里它仅仅是以共享平台方式来使用。并且我讨厌必须要共享。:-)所以我开始看看另外的选择。
Scilab,Octave,Sage...所有都有点脆弱,并且似乎没有我想要的功能特点和丰富性。之后我发现了 Pylab 。

Pylab 是 Python 环境的科学计算,包含了以下的包:

  1. matplotlib:图形和数据可视化;
  2. numpy:基本的数值分析(向量,矩阵,针对这些运算的科学函数);
  3. scipy:科学和工程应用。

等等。这是一篇嵌入式系统的博文,对吗?!Python 不能运行于资源-有限的嵌入式系统,并且实际上,Python 是符合我的三个标准之一的。。。可惜你不能运行在资源有限的嵌入式系统:

  1. 你具有命令行的操作系统
  2. 你可以运行 Python
  3. 有编译器运行在你的操作系统中,所以你不必需要交叉-编译

所以如果你正在使用 Python,你不会真正做嵌入式系统的开发,不过这个没事,因为你需要拓展你的技能跨度。不要做一个只会一技之长的人而且只会用 C,为你喜欢的处理器选择集成开发!

不管怎样,有许多时间需要我停止编程而需要对我的某些想法理论化。后来Python给了很大的帮助。

我们真的需要臭恶的 MATLAB 吗?

我们需要清楚的是本篇针对的是工程师(尤其是嵌入式系统的工程师),他们的信号处理,数据分析和可视化工作是作为他们工作的次要部分而言的。

对于全职且一直做信号处理或控制系统设计的工程师,MATLAB 是合适的选择。
如果您的公司有能力支付每周 40 小时的费用,他们也可以负担得起MATLAB的费用。

如果对费用不关系,我喜欢使用 MATLAB,并且我会尽可能的拿到所有的工具箱。

我不会在这里深入阐述信号处理或控制系统算法(z-变换,FFTs,根轨迹图,Nichols 图等等)。我会一步步的对使用 Python 和 Pylab 进行介绍。Pylab 的基本使用纯粹是激发你们的兴趣。

应用例子

假设你需要理解具有有感负荷的H-bridge的波纹电流,在边缘对齐和中心对齐的脉冲宽度调制。
这里有一些波纹电流图,是用一些Python脚本语言产生的。

边缘对齐PWM(脉冲宽度调制)

file

中心对齐PWM(脉冲宽度调制)

file

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中手动画的)。
file
file
相应的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上自己亲自试试:
file

安装

Python 核心的安装是非常简单的;OSX 系统用户可以直接安装 Python,但是不管你是什么操作系统,在 Python 官网 python.org 有编译好的二进制安装文件。

如果你想安装 scipy / numpy / matplotlib 库而不依赖于安装了正确的编译环境,那么这事情会变得有点棘手。

scipy.org 网站上列出了一些很好的解决方案;我原以为我也会分享自己的经验。大师我没有使用 Linux 所以请查看 scipy.org 页面的解决方案。

WINDOWS系统

我使用了三种免费的 PyLab 预打包版本

  1. Enthought Canopy
  2. PortablePython
  3. PythonXY

PortablePython 具有最可靠的安装/运行时。
PythonXY 具有最大的功能集(以及最大的安装占用空间)。
Enthought Canopy 不错;。Enthought 提供免费的版本试用,如果你想要更多的库,可以购买非免费版本 - 比如他们早期的发行版,EPD。这些版本从命令行运行起来有点容易,但我不知道如何稳定地跑起来。
还有 Anaconda,我一开始在 Mac OSX 系统上用过,但没在 Windows 上用过。

MAC OSX系统

我在家里的 Mac 上运行 Snow Leopard(OSX 10.6)。我还没有为 PyLab 找到一个很好的解决方案,但我正在努力。

PyLab 最简单的免费安装似乎是来自 Continuum Analytics 的 Anaconda。安装很简单,很快就能工作...除了我运行为这篇文章编写的脚本(脚本确实正常工作)时有一些关于内存分配的警告。当我去运行我常规的 Python 安装时,我的matplotlib 安装搞砸了。哎呀,希望这些问题能得到理顺。Anaconda 貌似很有前景。

Mac 上常用的免费软件进程是使用像 fink 或 MacPorts 这样的包管理器。 MacPorts 进程(sudo port install blahblahblah ...来自命令终端)有点脆弱。如果你的设置有问题,那么整个过程就会停止,并且带有一个神秘的信息。

Enthought Canopy 也有 OSX 和 Linux 版本,但我还没有试过。

也可以使用已存在的预编译二进制文件用于各种包。虽然 Python 是预安装在 Mac 上的,但请确保您的 Python 版本与您要安装的库兼容。我也建议安装一个最新版本的 Python 。至少,这些是你需要的:

  1. Python
  2. matplotlib
  3. scipy
  4. numpy

原文地址:https://www.dsprelated.com/showarticle/359/adventures-in-signal-processing-with-python-matlab-we-don-t-need-no-stinkin-matlab

追求梦想,做最好的自己