# https://docs.scipy.org/doc/scipy-0.17.0/reference/generated/scipy.signal.welch.html
from scipy import signal
import matplotlib.pyplot as plt
np.random.seed(1234)
fs = 10e3
N = 1e5
amp = 2*np.sqrt(2)
freq = 1234.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
x = amp*np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
f, Pxx_den = signal.welch(x, fs, nperseg=1024)
plt.semilogy(f, Pxx_den)
plt.ylim([0.5e-3, 1])
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD [V**2/Hz]')
plt.show()
f, Pxx_spec = signal.welch(x, fs, 'flattop', 1024, scaling='spectrum')
plt.figure()
plt.semilogy(f, np.sqrt(Pxx_spec))
plt.xlabel('frequency [Hz]')
plt.ylabel('Linear spectrum [V RMS]')
plt.show()
from math import cos
from scipy import optimize
def f(x):
d = 140
l = 156
a, r = x.tolist()
return [
cos(a) - 1 + (d*d)/(2*r*r),
l - r * a
]
result = optimize.fsolve(f,[2, 2])
print result
print f(result)
from scipy import optimize
def f(x):
return x**2 + 10 * np.sin(x)
x = np.arange(-10, 10, 0.1)
plt.plot(x, f(x))
plt.show()
# 梯度下降法求最小值
optimize.fmin_bfgs(f, 0)
# optimize.fmin_bfgs(f, 3, disp=0) 初始值=3,找到局部最小值
# 给定网格,暴力搜索
grid = (-10, 10, 0.1)
xmin_global = optimize.brute(f, (grid,))
# 模拟退火:scipy.optimize.anneal()
# 限制搜索区间
optimize.fminbound(f, 0, 10)
# minmize:
# -7x1+7x2-2x3-x4-6x5
# s.t.:
# 3x1-x2+x3-2x4=-3
# 2x1+x2+x4+x5=4
# -x1+3x2-3x4+x6=12
# xi>=0
import numpy as np
import matplotlib.pyplot as mpl
from scipy import optimize
c=np.array([-7,7,-2,-1,-6,0])
a=np.array([[3,-1,1,-2,0,0],[2,1,0,1,1,0],[-1,3,0,-3,0,1]])
b=np.array([-3,4,12])
res=optimize.linprog(c,A_eq=a,b_eq=b,bounds=((0,None),(0,None),(0,None),(0,None),(0,None),(0,None)))
print (res.x)
print (res.fun)
# maxmize:
# 2x1+3x2-5x3
# s.t.:
# x1+x2+x3=7
# 2x1-5x2+x3>=10
# x1+3x2+x3<=12
# xi>=0
from scipy import optimize as op
import numpy as np
c=np.array([2,3,-5])
A_ub=np.array([[-2,5,-1],[1,3,1]])
B_ub=np.array([-10,12])
A_eq=np.array([[1,1,1]])
B_eq=np.array([7])
x1=(0,7)
x2=(0,7)
x3=(0,7)
res=op.linprog(-c,A_ub,B_ub,A_eq,B_eq,bounds=(x1,x2,x3))
print(res)