计算与三角函数的方法傅立叶级数

Question

我试着按照下面的公式来实现傅立叶级数功能：

...其中...

...和...

这里是我的问题解决办法：

import numpy as np 
import pylab as py 

# Define "x" range. 
x = np.linspace(0, 10, 1000) 

# Define "T", i.e functions' period. 
T = 2 
L = T/2 

# "f(x)" function definition. 
def f(x): 
    return np.sin(np.pi * 1000 * x) 

# "a" coefficient calculation. 
def a(n, L, accuracy = 1000): 
    a, b = -L, L 
    dx = (b - a)/accuracy 
    integration = 0 
    for i in np.linspace(a, b, accuracy): 
     x = a + i * dx 
     integration += f(x) * np.cos((n * np.pi * x)/L) 
    integration *= dx 
    return (1/L) * integration 

# "b" coefficient calculation. 
def b(n, L, accuracy = 1000): 
    a, b = -L, L 
    dx = (b - a)/accuracy 
    integration = 0 
    for i in np.linspace(a, b, accuracy): 
     x = a + i * dx 
     integration += f(x) * np.sin((n * np.pi * x)/L) 
    integration *= dx 
    return (1/L) * integration 

# Fourier series. 
def Sf(x, L, n = 10): 
    a0 = a(0, L) 
    sum = 0 
    for i in np.arange(1, n + 1): 
     sum += ((a(i, L) * np.cos(n * np.pi * x)) + (b(i, L) * np.sin(n * np.pi * x))) 
    return (a0/2) + sum  

# x axis. 
py.plot(x, np.zeros(np.size(x)), color = 'black') 

# y axis. 
py.plot(np.zeros(np.size(x)), x, color = 'black') 

# Original signal. 
py.plot(x, f(x), linewidth = 1.5, label = 'Signal') 

# Approximation signal (Fourier series coefficients). 
py.plot(x, Sf(x, L), color = 'red', linewidth = 1.5, label = 'Fourier series') 

# Specify x and y axes limits. 
py.xlim([0, 10]) 
py.ylim([-2, 2]) 

py.legend(loc = 'upper right', fontsize = '10') 

py.show() 

...这里是绘制后的结果我得到什么：

我读过How to calculate a Fourier series in Numpy?，我已经实施了这种方法。它工作的很好，但它使用expotential方法，我想集中在三角函数和矩形方法的情况下计算a_{n}和b_{n}系数的integraions。

预先感谢您。

UPDATE（解决）

最后，这里是代码的工作示例。但是，我会花更多时间在它上面，所以如果有什么可以改进的地方，那就完成了。

from __future__ import division 
import numpy as np 
import pylab as py 

# Define "x" range. 
x = np.linspace(0, 10, 1000) 

# Define "T", i.e functions' period. 
T = 2 
L = T/2 

# "f(x)" function definition. 
def f(x): 
    return np.sin((np.pi) * x) + np.sin((2 * np.pi) * x) + np.sin((5 * np.pi) * x) 

# "a" coefficient calculation. 
def a(n, L, accuracy = 1000): 
    a, b = -L, L 
    dx = (b - a)/accuracy 
    integration = 0 
    for x in np.linspace(a, b, accuracy): 
     integration += f(x) * np.cos((n * np.pi * x)/L) 
    integration *= dx 
    return (1/L) * integration 

# "b" coefficient calculation. 
def b(n, L, accuracy = 1000): 
    a, b = -L, L 
    dx = (b - a)/accuracy 
    integration = 0 
    for x in np.linspace(a, b, accuracy): 
     integration += f(x) * np.sin((n * np.pi * x)/L) 
    integration *= dx 
    return (1/L) * integration 

# Fourier series. 
def Sf(x, L, n = 10): 
    a0 = a(0, L) 
    sum = np.zeros(np.size(x)) 
    for i in np.arange(1, n + 1): 
     sum += ((a(i, L) * np.cos((i * np.pi * x)/L)) + (b(i, L) * np.sin((i * np.pi * x)/L))) 
    return (a0/2) + sum 

# x axis. 
py.plot(x, np.zeros(np.size(x)), color = 'black') 

# y axis. 
py.plot(np.zeros(np.size(x)), x, color = 'black') 

# Original signal. 
py.plot(x, f(x), linewidth = 1.5, label = 'Signal') 

# Approximation signal (Fourier series coefficients). 
py.plot(x, Sf(x, L), '.', color = 'red', linewidth = 1.5, label = 'Fourier series') 

# Specify x and y axes limits. 
py.xlim([0, 5]) 
py.ylim([-2.2, 2.2]) 

py.legend(loc = 'upper right', fontsize = '10') 

py.show() 

Answer 1

考虑块开发代码以不同的方式，块。如果像这样的代码在第一次尝试时就会有效，你应该感到惊讶。 @ tom10说，调试是一种选择。另一种选择是在解释器中逐步快速构建代码原型，ipython甚至更好。

以上，您期望b_1000非零，因为输入f(x)是一个正弦曲线，其中有1000。你也期待所有其他系数都是零的权利？

然后您应该只关注功能b(n, L, accuracy = 1000)。看着它，有三件事情出错了。这里有一些提示。

• dx的乘法在循环内。对吗？
• 在循环中，i应该是一个整数吗？它真的是一个整数吗？通过原型或调试，你会发现这
• 每当你写(1/L)或类似的表达时要小心。如果你使用python2.7，你可能会犯错。如果没有，至少在源代码顶部使用from __future__ import division。如果你不知道我在说什么，请阅读this PEP。

如果你解决了这3点，b()将工作。然后以类似的方式考虑a。