DSP

Python 在信号处理中的优势

2019-07-13 11:18发布

休息了几天回来了

前言

本篇是对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