Numpy中的傅里叶级数。关于之前回答的问题

Question

我正在尝试复制之前一个帖子里的答案：如何在Numpy中计算傅里叶级数？

import numpy as np
import matplotlib.pyplot as plt
import itertools

def func(x):
    if x >= 1.0 or x <= -1.0:
        return 0
    else:
        return (abs(x) - 1.0)
a = 1.0
b = -1.0
N = 128.
time = np.linspace( a, b, N )
y = (np.fromiter(itertools.imap(func, time), 
                  dtype=time.dtype, count=time.shape[0]))
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(time,y) 
period = 2.
def cn(n):
    c = y*np.exp(-1j*2*n*np.pi*time/period)
    return c.sum()/c.size
def f(x, Nh):
    f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/period) for i in range(1,Nh+1)])
    return f.sum()
y2 = np.array([f(t,10).real for t in time])
ax.plot(time, y2)
plt.show()

我得到的结果接近正确答案，但有点偏移。我不太确定我哪里做错了。

Answer 1

你可以用下面这种方式来优化你的代码：

import numpy as np
import matplotlib.pyplot as plt
from optparse import OptionParser

def func(x):
    return np.where(np.abs(x) >= 1, 0., np.abs(x) - 1.0)

def cn(x, y, n, period):
    c = y * np.exp(-1j * 2. * np.pi * n * x / period)
    return c.sum()/c.size

def f(x, y, Nh, period):
    rng = np.arange(.5, Nh+.5)
    coeffs = np.array([cn(x,y,i,period) for i in rng])
    f = np.array([2. * coeffs[i] * np.exp(1j*2*i*np.pi*x/period) for i in rng])
    return f.sum(axis=0)

if __name__=='__main__':

    Version = '0.1'
    usage = "usage: %prog [options]"

    parser = OptionParser(usage = usage,version="%prog "+Version)
    parser.add_option("-a", dest='a', type='float', default=1., help="initial time")
    parser.add_option("-b", dest='b', type='float', default=-1., help="end time")
    parser.add_option("-N", "--Nt", dest='N', type='int', default=128, help="number of time steps")
    parser.add_option("-p", "--period", dest='period', type='float', default=2., help="period [time span]")
    parser.add_option("--Nh", dest='Nh', type='int', default=10, help="number of fourier series terms")

    (options, args) = parser.parse_args()

    for key,value in options.__dict__.iteritems():
        exec key + ' = ' + repr(value)

    time = np.linspace( a, b, N )
    y = func(time)
    period = np.abs(a-b)

    y2 = f(time,y,Nh,period).real

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.plot(time, y)
    ax.plot(time, y2)
    plt.show()

假设你的代码已经保存为"fourier_series.py"这个名字，你可以在普通的终端里尝试输入：

python fourier_series.py -N 512 --Nh 128

或者在ipython控制台里输入：

%run fourier_series.py -N 512 --Nh 128

Answer 2

我觉得你的 DC 项目好像出问题了。我现在不能自己检查，但你确定在 f() 函数里，range(1, Nh+1) 里的 1 是正确的吗？

Answer 3

这个错误似乎和你的黎曼和 计算方法 有关（右边/中间/左边），这是 regularfry 提到的。使用中间方法会得到：

代码：

import numpy as np
import matplotlib.pyplot as plt
import itertools

def func(x):
    if x >= 1.0 or x <= -1.0:
        return 0
    else:
        return (abs(x) - 1.0)

a = 1.0
b = -1.0
N = 128.
time = np.linspace( a, b, N )
y = (np.fromiter(itertools.imap(func, time), 
                  dtype=time.dtype, count=time.shape[0]))

period = 2.
def cn(n):
    c = y*np.exp(-1j*2*n*np.pi*time/period)
    return c.sum()/c.size
def f(x, Nh):
    rng = np.arange(.5, Nh+.5)
    f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/period) for i in rng])
    return f.sum()

y2 = np.array([f(t,10).real for t in time])

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(time, y)
ax.plot(time, y2)
plt.show()

正如Sven在 另一个问题 中提到的，你使用列表推导（还有 imap）而不是数组和 ufuncs，这样做效率比较低（如果你遇到性能问题的话）。