求傅里叶系数算法
Finding Fourier coefficients algorithm
好的，所以我一直在尝试编写一个 "naive" 方法来计算复数形式的标准傅立叶级数的系数。我想我已经非常接近了，但是有一些奇怪的行为。这可能更像是一道数学题而不是编程题，但我 already asked 在 math.stackexchange 上得到了零答案。这是我的工作代码：
import matplotlib.pyplot as plt
import numpy as np
def coefficients(fn, dx, m, L):
"""
Calculate the complex form fourier series coefficients for the first M
waves.
:param fn: function to sample
:param dx: sampling frequency
:param m: number of waves to compute
:param L: We are solving on the interval [-L, L]
:return: an array containing M Fourier coefficients c_m
"""
N = 2*L / dx
coeffs = np.zeros(m, dtype=np.complex_)
xk = np.arange(-L, L + dx, dx)
# Calculate the coefficients for each wave
for mi in range(m):
coeffs[mi] = 1/N * sum(fn(xk)*np.exp(-1j * mi * np.pi * xk / L))
return coeffs
def fourier_graph(range, L, c_coef, function=None, plot=True, err_plot=False):
"""
Given a range to plot and an array of complex fourier series coefficients,
this function plots the representation.
:param range: the x-axis values to plot
:param c_coef: the complex fourier coefficients, calculated by coefficients()
:param plot: Default True. Plot the fourier representation
:param function: For calculating relative error, provide function definition
:param err_plot: relative error plotted. requires a function to compare solution to
:return: the fourier series values for the given range
"""
# Number of coefficients to sum over
w = len(c_coef)
# Initialize solution array
s = np.zeros(len(range))
for i, ix in enumerate(range):
for iw in np.arange(w):
s[i] += c_coef[iw] * np.exp(1j * iw * np.pi * ix / L)
# If a plot is desired:
if plot:
plt.suptitle("Fourier Series Plot")
plt.xlabel(r"$t$")
plt.ylabel(r"$f(x)$")
plt.plot(range, s, label="Fourier Series")
if err_plot:
plt.plot(range, function(range), label="Actual Solution")
plt.legend()
plt.show()
# If error plot is desired:
if err_plot:
err = abs(function(range) - s) / function(range)
plt.suptitle("Plot of Relative Error")
plt.xlabel("Steps")
plt.ylabel("Relative Error")
plt.plot(range, err)
plt.show()
return s
if __name__ == '__main__':
# Assuming the interval [-l, l] apply discrete fourier transform:
# number of waves to sum
wvs = 50
# step size for calculating c_m coefficients (trap rule)
deltax = .025 * np.pi
# length of interval for Fourier Series is 2*l
l = 2 * np.pi
c_m = coefficients(np.exp, deltax, wvs, l)
# The x range we would like to interpolate function values
x = np.arange(-l, l, .01)
sol = fourier_graph(x, l, c_m, np.exp, err_plot=True)
现在，每个系数乘以2/N。但是，我在我教授打字的笔记中推导了这个总和，其中不包括这个因子 2/N。当我导出 the form myself 时，我得到了一个公式，其中的因子为 1/N，无论我尝试什么技巧都不会抵消。我在 math.stackexchange 询问发生了什么，但没有得到任何答案。
我确实注意到，添加 1/N 后，实际解与傅里叶级数之间的差异大大减小了，但仍然不对。所以我尝试了 2/N 并获得了更好的结果。我真的很想弄清楚这一点，这样我就可以在尝试学习快速傅立叶变换之前为基本傅立叶级数编写一个漂亮、干净的算法。
那我做错了什么？
假设 c_n 由 A_n 给出，如 mathworld
同上c_n = 1/T \int_{-T/2}^{T/2}f(x)e^{-2ipinx/T}dx
我们可以（简单地）计算系数 c_n 分析（这是与梯形积分进行比较的好方法）
k = (1-2in)/2
c_n = 1/(4*pi*k)*(e^{2pik} - e^{-2pik})
所以你的系数很可能被正确计算（两条错误的曲线看起来很相似）
现在请注意，当您重构 f 时，您将系数 c_0 添加到 c_m
但是重建应该发生在 c_{-m} 到 c_m
所以你漏掉了一半的系数。
以下是您修改后的系数函数和理论系数的修正
import matplotlib.pyplot as plt
import numpy as np
def coefficients(fn, dx, m, L):
"""
Calculate the complex form fourier series coefficients for the first M
waves.
:param fn: function to sample
:param dx: sampling frequency
:param m: number of waves to compute
:param L: We are solving on the interval [-L, L]
:return: an array containing M Fourier coefficients c_m
"""
N = 2*L / dx
coeffs = np.zeros(m, dtype=np.complex_)
xk = np.arange(-L, L + dx, dx)
# Calculate the coefficients for each wave
for mi in range(m):
n = mi - m/2
coeffs[mi] = 1/N * sum(fn(xk)*np.exp(-1j * n * np.pi * xk / L))
return coeffs
def fourier_graph(range, L, c_coef, ref, function=None, plot=True, err_plot=False):
"""
Given a range to plot and an array of complex fourier series coefficients,
this function plots the representation.
:param range: the x-axis values to plot
:param c_coef: the complex fourier coefficients, calculated by coefficients()
:param plot: Default True. Plot the fourier representation
:param function: For calculating relative error, provide function definition
:param err_plot: relative error plotted. requires a function to compare solution to
:return: the fourier series values for the given range
"""
# Number of coefficients to sum over
w = len(c_coef)
# Initialize solution array
s = np.zeros(len(range), dtype=complex)
t = np.zeros(len(range), dtype=complex)
for i, ix in enumerate(range):
for iw in np.arange(w):
n = iw - w/2
s[i] += c_coef[iw] * (np.exp(1j * n * ix * 2 * np.pi / L))
t[i] += ref[iw] * (np.exp(1j * n * ix * 2 * np.pi / L))
# If a plot is desired:
if plot:
plt.suptitle("Fourier Series Plot")
plt.xlabel(r"$t$")
plt.ylabel(r"$f(x)$")
plt.plot(range, s, label="Fourier Series")
plt.plot(range, t, label="expected Solution")
plt.legend()
if err_plot:
plt.plot(range, function(range), label="Actual Solution")
plt.legend()
plt.show()
return s
def ref_coefficients(m):
"""
Calculate the complex form fourier series coefficients for the first M
waves.
:param fn: function to sample
:param dx: sampling frequency
:param m: number of waves to compute
:param L: We are solving on the interval [-L, L]
:return: an array containing M Fourier coefficients c_m
"""
coeffs = np.zeros(m, dtype=np.complex_)
# Calculate the coefficients for each wave
for iw in range(m):
n = iw - m/2
k = (1-(1j*n)/2)
coeffs[iw] = 1/(4*np.pi*k)* (np.exp(2*np.pi*k) - np.exp(-2*np.pi*k))
return coeffs
if __name__ == '__main__':
# Assuming the interval [-l, l] apply discrete fourier transform:
# number of waves to sum
wvs = 50
# step size for calculating c_m coefficients (trap rule)
deltax = .025 * np.pi
# length of interval for Fourier Series is 2*l
l = 2 * np.pi
c_m = coefficients(np.exp, deltax, wvs, l)
# The x range we would like to interpolate function values
x = np.arange(-l, l, .01)
ref = ref_coefficients(wvs)
sol = fourier_graph(x, 2*l, c_m, ref, np.exp, err_plot=True)
好的，所以我一直在尝试编写一个 "naive" 方法来计算复数形式的标准傅立叶级数的系数。我想我已经非常接近了，但是有一些奇怪的行为。这可能更像是一道数学题而不是编程题，但我 already asked 在 math.stackexchange 上得到了零答案。这是我的工作代码：
import matplotlib.pyplot as plt
import numpy as np
def coefficients(fn, dx, m, L):
"""
Calculate the complex form fourier series coefficients for the first M
waves.
:param fn: function to sample
:param dx: sampling frequency
:param m: number of waves to compute
:param L: We are solving on the interval [-L, L]
:return: an array containing M Fourier coefficients c_m
"""
N = 2*L / dx
coeffs = np.zeros(m, dtype=np.complex_)
xk = np.arange(-L, L + dx, dx)
# Calculate the coefficients for each wave
for mi in range(m):
coeffs[mi] = 1/N * sum(fn(xk)*np.exp(-1j * mi * np.pi * xk / L))
return coeffs
def fourier_graph(range, L, c_coef, function=None, plot=True, err_plot=False):
"""
Given a range to plot and an array of complex fourier series coefficients,
this function plots the representation.
:param range: the x-axis values to plot
:param c_coef: the complex fourier coefficients, calculated by coefficients()
:param plot: Default True. Plot the fourier representation
:param function: For calculating relative error, provide function definition
:param err_plot: relative error plotted. requires a function to compare solution to
:return: the fourier series values for the given range
"""
# Number of coefficients to sum over
w = len(c_coef)
# Initialize solution array
s = np.zeros(len(range))
for i, ix in enumerate(range):
for iw in np.arange(w):
s[i] += c_coef[iw] * np.exp(1j * iw * np.pi * ix / L)
# If a plot is desired:
if plot:
plt.suptitle("Fourier Series Plot")
plt.xlabel(r"$t$")
plt.ylabel(r"$f(x)$")
plt.plot(range, s, label="Fourier Series")
if err_plot:
plt.plot(range, function(range), label="Actual Solution")
plt.legend()
plt.show()
# If error plot is desired:
if err_plot:
err = abs(function(range) - s) / function(range)
plt.suptitle("Plot of Relative Error")
plt.xlabel("Steps")
plt.ylabel("Relative Error")
plt.plot(range, err)
plt.show()
return s
if __name__ == '__main__':
# Assuming the interval [-l, l] apply discrete fourier transform:
# number of waves to sum
wvs = 50
# step size for calculating c_m coefficients (trap rule)
deltax = .025 * np.pi
# length of interval for Fourier Series is 2*l
l = 2 * np.pi
c_m = coefficients(np.exp, deltax, wvs, l)
# The x range we would like to interpolate function values
x = np.arange(-l, l, .01)
sol = fourier_graph(x, l, c_m, np.exp, err_plot=True)
现在，每个系数乘以2/N。但是，我在我教授打字的笔记中推导了这个总和，其中不包括这个因子 2/N。当我导出 the form myself 时，我得到了一个公式，其中的因子为 1/N，无论我尝试什么技巧都不会抵消。我在 math.stackexchange 询问发生了什么，但没有得到任何答案。
我确实注意到，添加 1/N 后，实际解与傅里叶级数之间的差异大大减小了，但仍然不对。所以我尝试了 2/N 并获得了更好的结果。我真的很想弄清楚这一点，这样我就可以在尝试学习快速傅立叶变换之前为基本傅立叶级数编写一个漂亮、干净的算法。
那我做错了什么？
假设 c_n 由 A_n 给出，如 mathworld
同上c_n = 1/T \int_{-T/2}^{T/2}f(x)e^{-2ipinx/T}dx
我们可以（简单地）计算系数 c_n 分析（这是与梯形积分进行比较的好方法）
k = (1-2in)/2
c_n = 1/(4*pi*k)*(e^{2pik} - e^{-2pik})
所以你的系数很可能被正确计算（两条错误的曲线看起来很相似）
现在请注意，当您重构 f 时，您将系数 c_0 添加到 c_m
但是重建应该发生在 c_{-m} 到 c_m
所以你漏掉了一半的系数。
以下是您修改后的系数函数和理论系数的修正
import matplotlib.pyplot as plt
import numpy as np
def coefficients(fn, dx, m, L):
"""
Calculate the complex form fourier series coefficients for the first M
waves.
:param fn: function to sample
:param dx: sampling frequency
:param m: number of waves to compute
:param L: We are solving on the interval [-L, L]
:return: an array containing M Fourier coefficients c_m
"""
N = 2*L / dx
coeffs = np.zeros(m, dtype=np.complex_)
xk = np.arange(-L, L + dx, dx)
# Calculate the coefficients for each wave
for mi in range(m):
n = mi - m/2
coeffs[mi] = 1/N * sum(fn(xk)*np.exp(-1j * n * np.pi * xk / L))
return coeffs
def fourier_graph(range, L, c_coef, ref, function=None, plot=True, err_plot=False):
"""
Given a range to plot and an array of complex fourier series coefficients,
this function plots the representation.
:param range: the x-axis values to plot
:param c_coef: the complex fourier coefficients, calculated by coefficients()
:param plot: Default True. Plot the fourier representation
:param function: For calculating relative error, provide function definition
:param err_plot: relative error plotted. requires a function to compare solution to
:return: the fourier series values for the given range
"""
# Number of coefficients to sum over
w = len(c_coef)
# Initialize solution array
s = np.zeros(len(range), dtype=complex)
t = np.zeros(len(range), dtype=complex)
for i, ix in enumerate(range):
for iw in np.arange(w):
n = iw - w/2
s[i] += c_coef[iw] * (np.exp(1j * n * ix * 2 * np.pi / L))
t[i] += ref[iw] * (np.exp(1j * n * ix * 2 * np.pi / L))
# If a plot is desired:
if plot:
plt.suptitle("Fourier Series Plot")
plt.xlabel(r"$t$")
plt.ylabel(r"$f(x)$")
plt.plot(range, s, label="Fourier Series")
plt.plot(range, t, label="expected Solution")
plt.legend()
if err_plot:
plt.plot(range, function(range), label="Actual Solution")
plt.legend()
plt.show()
return s
def ref_coefficients(m):
"""
Calculate the complex form fourier series coefficients for the first M
waves.
:param fn: function to sample
:param dx: sampling frequency
:param m: number of waves to compute
:param L: We are solving on the interval [-L, L]
:return: an array containing M Fourier coefficients c_m
"""
coeffs = np.zeros(m, dtype=np.complex_)
# Calculate the coefficients for each wave
for iw in range(m):
n = iw - m/2
k = (1-(1j*n)/2)
coeffs[iw] = 1/(4*np.pi*k)* (np.exp(2*np.pi*k) - np.exp(-2*np.pi*k))
return coeffs
if __name__ == '__main__':
# Assuming the interval [-l, l] apply discrete fourier transform:
# number of waves to sum
wvs = 50
# step size for calculating c_m coefficients (trap rule)
deltax = .025 * np.pi
# length of interval for Fourier Series is 2*l
l = 2 * np.pi
c_m = coefficients(np.exp, deltax, wvs, l)
# The x range we would like to interpolate function values
x = np.arange(-l, l, .01)
ref = ref_coefficients(wvs)
sol = fourier_graph(x, 2*l, c_m, ref, np.exp, err_plot=True)