python - 如何在numpy中向量化傅立叶级数部分和

Question

给定周期为 T 的函数的傅里叶级数系数​​ a[n] 和 b[n](分别对应余弦和正弦) > 和 t 一个等间隔的间隔下面的代码将计算间隔 t 中所有点的部分和 (a,b,t 都是 numpy 数组)。明确了 len(t) <> len(a).

yn=ones(len(t))*a[0]
for n in range(1,len(a)):
    yn=yn+(a[n]*cos(2*pi*n*t/T)-b[n]*sin(2*pi*n*t/T))

我的问题是:这个 for 循环可以向量化吗？

Answer 1

这是一种使用 broadcasting 创建余弦/正弦输入的 2D 数组版本的矢量化方法:2*pi*n*t/T 然后使用 matrix-multiplication 与 np.dot 用于 sum-reduction -

r = np.arange(1,len(a))
S = 2*np.pi*r[:,None]*t/T
cS = np.cos(S)
sS = np.sin(S)
out = a[1:].dot(cS) - b[1:].dot(sS) + a[0]

进一步提升性能

为了进一步提升，我们可以利用 numexpr module 来计算那些三角步 -

import numexpr as ne
cS = ne.evaluate('cos(S)')
sS = ne.evaluate('sin(S)')

运行时测试 -

方法-

def original_app(t,a,b,T):
    yn=np.ones(len(t))*a[0]
    for n in range(1,len(a)):
        yn=yn+(a[n]*np.cos(2*np.pi*n*t/T)-b[n]*np.sin(2*np.pi*n*t/T))
    return yn

def vectorized_app(t,a,b,T):    
    r = np.arange(1,len(a))
    S = (2*np.pi/T)*r[:,None]*t
    cS = np.cos(S)
    sS = np.sin(S)
    return a[1:].dot(cS) - b[1:].dot(sS) + a[0]

def vectorized_app_v2(t,a,b,T):    
    r = np.arange(1,len(a))
    S = (2*np.pi/T)*r[:,None]*t
    cS = ne.evaluate('cos(S)')
    sS = ne.evaluate('sin(S)')
    return a[1:].dot(cS) - b[1:].dot(sS) + a[0]

此外，还包括来自@Paul Panzer 的帖子的函数 PP。

时间 -

In [22]: # Setup inputs
    ...: n = 10000
    ...: t = np.random.randint(0,9,(n))
    ...: a = np.random.randint(0,9,(n))
    ...: b = np.random.randint(0,9,(n))
    ...: T = 3.45
    ...: 

In [23]: print np.allclose(original_app(t,a,b,T), vectorized_app(t,a,b,T))
    ...: print np.allclose(original_app(t,a,b,T), vectorized_app_v2(t,a,b,T))
    ...: print np.allclose(original_app(t,a,b,T), PP(t,a,b,T))
    ...: 
True
True
True

In [25]: %timeit original_app(t,a,b,T)
    ...: %timeit vectorized_app(t,a,b,T)
    ...: %timeit vectorized_app_v2(t,a,b,T)
    ...: %timeit PP(t,a,b,T)
    ...: 
1 loops, best of 3: 6.49 s per loop
1 loops, best of 3: 6.24 s per loop
1 loops, best of 3: 1.54 s per loop
1 loops, best of 3: 1.96 s per loop