从numpy手动恢复原始函数

Question

我对一个函数执行了numpy.fft.rfft运算，以获得傅里叶系数。由于the docs似乎不包含所使用的精确公式，我一直假设在我的教科书中找到一个公式：

S(x) = a_0/2 + SUM(real(a_n) * cos(nx) + imag(a_n) * sin(nx))

其中imag(a_n)是傅里叶系数的n_th元素的虚部。你知道吗

为了将其转换为python speak，我实现了以下功能：

def fourier(freqs, X):
    # input the fourier frequencies from np.fft.rfft, and arbitrary X
    const_term = np.repeat(np.real(freqs[0])/2, X.shape[0]).reshape(-1,1)
    # this is the "n" part of the inside of the trig terms
    trig_terms = np.tile(np.arange(1,len(freqs)), (X.shape[0],1))
    sin_terms = np.imag(freqs[1:])*np.sin(np.einsum('i,ij->ij', X, trig_terms))
    cos_terms = np.real(freqs[1:])*np.cos(np.einsum('i,ij->ij', X, trig_terms))
    return np.concatenate((const_term, sin_terms, cos_terms), axis=1)

这应该给我一个[X.shape[0], 2*freqs.shape[0] - 1]数组，在i,j项中包含X的i_th元素，在傅里叶分解的j_th项处计算（其中j_th项是奇数j的sin项）。你知道吗

通过在傅里叶项的轴上求和这个数组，我应该得到在X中的i_th项处计算的函数：

import numpy as np
import matplotlib.pyplot as plt

X = np.linspace(-1,1,50)
y = X*(X-0.8)*(X+1)

reconstructed_y = np.sum(
    fourier(
        np.fft.rfft(y), 
        X
    ),
    axis = 1
)

plt.plot(X,y)
plt.plot(X, reconstructed_y, c='r')
plt.show()

无论如何，红线基本上应该在蓝线之上。在我关于numpy.fft.rfft返回什么的假设中，或者在我的具体实现中，出现了一些错误，但是我很难追踪到这个bug。有人能解释一下我做错了什么吗？你知道吗