SciPy是一种使用NumPy来做高等数学、信号处理、优化、统计和许多其它科学任务的语言扩展,SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。使用时要首先安装合适版本的scipy模块,win7_32位下载后一步步安装即可,直到安装完成。简单的例子:

from scipy import *
import numpy as np
from pylab import *  
a = zeros(1000)
a[:10]=1
b=fft(a)
plot(abs(b))
show()

信号 PYTHON 槽 scipy信号处理_信号 PYTHON 槽

(1)最小二乘法拟合数据

scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数leastsq。函数原型为:

leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=0.0, factor=100, diag=None, warning=True)


一般我们只要指定前三个参数就可以了。

func 是我们自己定义的一个计算误差的函数,

x0 是计算的初始参数值

args 是指定func的其他参数

下面用一个示例程序来解释他的用法,同样也是使用上面的求一次函数的拟合参数


# -*- coding: utf-8 -*-
import numpy as np
from scipy.optimize import leastsq
import pylab as pl

def func(x,p):
    """数据拟合所用的函数: A*sin(2*pi*k*x + theta)"""
    A, k, theta = p
    return A*np.sin(2*np.pi*k*x+theta)   

def errors(p,y,x):
    """实验数据x, y和拟合函数之间的差,p为拟合需要找到的系数"""
    return y-func(x,p)

x = np.linspace(0,-2*np.pi,100)
A, k, theta = 10,0.34,np.pi/6 # 真实数据的函数参数
y0 = func(x,[A,k,theta]) # 真实数据
y1 = y0+2*np.random.randn(len(x)) # 加入噪声之后的实验数据    

p0 = [7,0.2,0] # 设置第一次猜测的函数拟合参数

#调用leastsq进行数据拟合,errors为计算误差的函数,p0为拟合参数的初始值
#args为需要拟合的实验数据
plsq =leastsq(errors, p0, args=(y1, x))

print("real pares:", [A, k, theta]) 
print("fitting paras", plsq[0]) # 实验数据拟合后的参数

pl.plot(x, y0, label="real data")
pl.plot(x, y1, label="noise data")
pl.plot(x, func(x,plsq[0]),label="fitting data")
pl.legend()
pl.show()

F5运行结果如下:

>>> 
('real pares:', [10, 0.34, 0.5235987755982988])
('fitting paras', array([-10.51740809,   0.34191249,   3.66939538]))

信号 PYTHON 槽 scipy信号处理_信号 PYTHON 槽_02

注:由于正弦函数具有周期性,实际上拟合参数得到的函数和真实参数对应的函数基本是一致的。

(2)非线性方程组求解

optimize库中的fsolve函数可以用来对非线性方程组进行求解。它的基本调用形式如下:


fsolve(func, x0)

func(x)是计算方程组误差的函数,它的参数x是一个矢量,表示方程组的各个未知数的一组可能解,func返回将x代入方程组之后得到的误差;x0为未知数矢量的初始值。如果要对如下方程组进行求解的话:


f1(u1,u2,u3) = 0


f2(u1,u2,u3) = 0


f3(u1,u2,u3) = 0


那么func可以如下定义:


def func(x):
    u1,u2,u3 = x
    return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]

例如,求解如下方程组的解:


x1-2 = 0


x0*x0-*sin(x1*x2) = 0


x1*x2-2 = 0


from scipy.optimize import fsolve
from math import sin,cos

def f(x):
    x0 = float(x[0])
    x1 = float(x[1])
    x2 = float(x[2])
    return [
        x1-2,
        x0*x0-sin(x1*x2),
        x1*x2-2
    ]

result = fsolve(f, [1,1,1]) #设置初始值为[1 1 1]

print (result)
print (f(result))

运行结果如下:


>>> 
[ 0.95357088  2.          1.        ]
[0.0, -6.661338147750939e-16, -4.440892098500626e-16]


参考:

http://www.cnphp6.com/archives/65841

http://old.sebug.net/paper/books/scipydoc/scipy_intro.html