1

我正在尝试绘制方波的傅立叶序列，但是对于许多术语，该程序只需要太多时间来计算所有点。我尝试使用多处理模块但没有工作。请帮助我如何使用多处理。我在 fedora linux 上运行并且有一个 AMD FX Octa 内核。谢谢

#!/usr/bin/env python

import pylab,math

#Square Wave approximation by Fourier Series:
#y=4/pi*(sin(x) + sin(3x)/3 + sin(5x)/5 + ... n) terms

n=input("how many terms do you want?")
y=[]
# generate x's to calculate the points. 
x=pylab.arange(0,360,0.5)
#generate de odd numbers to add the sines terms
impar=range(1,n+2,2)
no=4/math.pi #just for normalize the sequence to converge to 1 (not pi/4)
#calculate the points
for ps in x:
    a=0
    for i in impar:
        a=a+(((math.sin(math.radians(i*ps)))/i))
    a=no*a
    y.append(a)

#plot
pylab.figure()
pylab.plot(x,y,"o")
pylab.title(str(n)+str(" terms"))
pylab.grid()
#save a image(just in case)
pylab.savefig(str(n)+str("sqwv.png"))
#show the graph
pylab.show()
4

2 回答 2

1

您可以使用 numpy 的 ufunc 来加速计算：

y = reduce(np.add, (np.sin(np.radians(i*x))/i for i in impar))

这是完整的代码：

import pylab,math
import numpy as np
#Square Wave approximation by Fourier Series:
#y=4/pi*(sin(x) + sin(3x)/3 + sin(5x)/5 + ... n) terms

n = 1000

# generate x's to calculate the points. 
x=pylab.arange(0,360,0.5)
#generate de odd numbers to add the sines terms
impar=range(1,n+2,2)
no=4/math.pi #just for normalize the sequence to converge to 1 (not pi/4)
#calculate the points

y = reduce(np.add, (np.sin(np.radians(i*x))/i for i in impar)) * no

#plot
pylab.figure()
pylab.plot(x,y,"-")
pylab.title(str(n)+str(" terms"))
pylab.grid()
#save a image(just in case)
pylab.savefig(str(n)+str("sqwv.png"))
#show the graph
pylab.show()

我的一台电脑，当 n 为 10000 时大约需要 200 毫秒。

于 2013-06-17T01:43:24.963 回答
1

这个问题是令人尴尬的并行，所以很容易让脚本使用多个 CPU。以@HYRY 的回答为基础：

from multiprocessing.dummy import Pool # use threads

def compute_y(i):
    np.add.reduce(no * np.sin(np.radians(impar * x[i])) / impar, out=y[i])

step = 10**8 // n 
Pool().map(compute_y, (slice(j, j + step) for j in range(0, len(x), step)))

矢量化 numpy 操作释放 GIL，因此该代码应该能够利用多个 CPU。

完整脚本：

#!/usr/bin/env python
"""Square Wave approximation by Fourier Series."""
import math
from multiprocessing.dummy import Pool # use threads
from timeit import default_timer as timer

import matplotlib.pyplot as plt
import numpy as np

n = 100000 # compute first n terms in the series

# generate x's to calculate the points.
x = np.arange(0, 360, .1)

# generate odd numbers to add the sines terms
impar = np.c_[1:n+2:2]
no = 4 / math.pi # just for normalize the sequence to converge to 1 (not pi/4)

# calculate the points
start = timer()
# y = 4/pi*(sin(x) + sin(3x)/3 + sin(5x)/5 + ... n) terms
y = np.empty_like(x)

def compute_y(i):
    t = impar * x[i]
    np.radians(t, out=t)
    np.sin(t, out=t)
    t *= no
    t /= impar
    np.add.reduce(t, out=y[i])

step = 10**8 // n
Pool().map(compute_y, (slice(j, j + step) for j in range(0, len(x), step)))
print("y calc. takes us %.2f seconds" % (timer() - start,))

# plot
plt.plot(x, y, "-")
plt.title("%d terms" % n)
plt.grid()
# save a image(just in case)
plt.savefig("%dsqwv.png" % n)
# show the graph
plt.show()
于 2013-06-17T03:34:25.683 回答