我尝试根据以下公式实现傅立叶级数函数:
...在哪里...
...和...
这是我解决问题的方法:
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?，并且已经实现了这种方法。它的效果很好，但是它使用了幂函数方法，在计算
a_{n}和b_{n}系数的积分时，我想重点介绍三角函数和矩形方法。先感谢您。
更新(已解决)
最后，这是代码的有效示例。但是，我会花更多的时间在上面，因此，如果有任何可以改进的地方，那就可以了。
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()
最佳答案
考虑以一种不同的方式(逐块地)开发代码。如果这样的代码在第一次尝试时会起作用，您应该会感到惊讶。正如@ tom10所说，调试是一种选择。另一个选择是在解释器中逐步逐步编写代码原型(prototype)，甚至使用ipython更好。
在上面，您期望b_1000不为零，因为输入的f(x)是一个正弦曲线，其中包含1000。您还期望所有其他系数都为零，对吗？
然后，您应该只关注b(n, L, accuracy = 1000)函数。看着它，三件事出了问题。这里有一些提示。
dx的乘法在循环内。可以吗i应该是整数吗？真的是整数吗？通过原型(prototype)设计或调试，您将发现此(1/L)或类似表达式时都要小心。如果您使用的是python2.7，则可能做错了。如果不是，请至少在源代码的顶部使用from __future__ import division。如果您不知道我在说什么，请阅读this PEP。 如果您解决了这3点，
b()将起作用。然后以类似的方式考虑a。关于python - 用三角法计算傅里叶级数，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/27708002/