实现欧拉形式的三角插值

Question

目前我正在努力正确实现欧拉的三角插值。为了指导您完成我已经完成的工作和收集的信息，我会将代码与数学定义配对。代码将用 python 编写，并使用 numpy 函数。 请注意，我不会使用快速傅立叶变换。

首先，所有实验都将在以下数据集上进行

(x,y)

：

[('0.0', '0.0'),
 ('0.6283185307179586', '0.6427876096865393'),
 ('1.2566370614359172', '0.984807753012208'),
 ('1.8849555921538759', '0.8660254037844387'),
 ('2.5132741228718345', '0.3420201433256689'),
 ('3.141592653589793', '-0.34202014332566866'),
 ('3.7699111843077517', '-0.8660254037844384'),
 ('4.39822971502571', '-0.9848077530122081'),
 ('5.026548245743669', '-0.6427876096865396'),
 ('5.654866776461628', '-2.4492935982947064e-16')]

我们定义离散傅里叶变换

所以我对这个功能的实现如下

import numpy as np #consider numpy as imported from here on

def F_n(Y):
    n = len(Y)
    Y_hat = []
    for k in range(len(Y)):
        transformed_k = 1/n * sum([y_l * np.exp(-2 * np.pi * 1j* k * l/n) for l, y_l in enumerate(Y) ])
        Y_hat.append(transformed_k)
    return Y_hat

因此，将得到的系数与使用

np.fft.fft

的系数进行比较，看起来系数几乎是正确的。它们的区别仅在于移动了一位数字，即

# F_n(y)
['(-1.33907057366955e-17+0j)',
 '(0.14283712054380923-0.439607454395199j)',
 '(-0.048591754799448425+0.06688081278992913j)',
 '(-0.039133572999081954+0.028432205056635337j)',
 '(-0.036913281031968816+0.01199385205986717j)',
 '(-0.036397023426620205-2.0058074207055733e-17j)',
 '(-0.03691328103196878-0.011993852059867215j)',
 '(-0.03913357299908168-0.028432205056635646j)',
 '(-0.04859175479944824-0.06688081278992904j)',
 '(0.1428371205438091+0.439607454395199j)']
# np.fft.fft(y)
['(-1.1102230246251565e-16+0j)',
 '(1.428371205438092-4.39607454395199j)',
 '(-0.4859175479944836+0.6688081278992911j)',
 '(-0.3913357299908192+0.2843220505663533j)',
 '(-0.36913281031968803+0.11993852059867194j)',
 '(-0.36397023426620184-1.1102230246251565e-16j)',
 '(-0.36913281031968803-0.11993852059867194j)',
 '(-0.3913357299908196-0.2843220505663534j)',
 '(-0.4859175479944836-0.6688081278992911j)',
 '(1.4283712054380922+4.39607454395199j)']

现在我想实现三角插值。这是定义

我用这个实现来实现定理

def trig_interpolation(Y_hat, x_range, depth=1000):
    n = len(Y_hat)

    get_summand = lambda c_j,l,x: c_j*np.exp(2 * np.pi * 1j * l*x)

    y_intp = []
    x_intp = list((i/depth)*x_range for i in range(depth))
    if n%2==0:
        K = n//2 
        for x in x_intp:
            y_intp.append(sum([get_summand(c_j,l,x) for l,c_j in zip(range(-K+1,K+1), Y_hat)]))
    else:
        K = n//2+1
        for x in x_intp:
            y_intp.append(sum([get_summand(c_j,l,x) for l,c_j in zip(range(-K,K+1), Y_hat)]))

    return x_intp, y_intp

x_range = max(x)-min(x)
x_intp, y_intp = trig_interpolation(np.fft.fft(y), x_range)

其中

x_range

表示原始数据点集中

x

值的范围，

depth

表示插值的分辨率，第 4 行中

get_summand

函数表示项 $$ xp(2\pi j x)$$ 使求和过程更容易阅读。

当对 numpy 的 fft 给出的系数运行我的代码时，我得到

plt.plot(x_intp,np.real(y_intp))
plt.plot(x,y, 'o')

虽然这些点与使用 numpy 的 fft 进行插值对齐，但我希望曲线看起来有所不同，表现得更接近实际的正弦曲线。

使用我的实现计算傅里叶系数给出了一条预期错误的曲线。

---

我要求指出我的错误，以便我可以根据我提出的数学特征正确地实现三角插值。

Answer 1

您所做的事情存在许多单独的问题。

最重要的是离散傅立叶变换不会产生您需要的频率顺序。顺序大致如下（参见 https://numpy.org/doc/stable/reference/routines.fft.html ）

freq[0] ....正频率...负频率（从非常正的频率映射）

所以你必须改变你的频率。有一个例程

numpy.fft.fftshift

但它并没有完全按照您想要的编号移动，所以我在这里使用

numpy.roll

来旋转数组。

您需要将 FFT 除以

n

才能获得所需形式的系数。

您的

x_range

是错误的 - 它基于隐含的 x_max，而不是数组的最后一个元素。

K

在一种情况下是错误的。

如果

x/x_range

不在 0 和 1 之间，则您希望在总和中包含

x

。

import numpy as np
import matplotlib.pyplot as plt

n = 10
x = np.linspace( 0.0, 2 * np.pi, n, endpoint=False )
y = np.sin( x )

def trig_interpolation(Y_hat, x_range, depth=1000):
    n = len(Y_hat)

    get_summand = lambda c_j, l, x: c_j * np.exp(2 * np.pi * 1j * l * x / x_range )

    x_intp = [(i / depth) * x_range for i in range(depth)]                     
    y_intp = []
    K = n // 2                                                       # Fix K for the odd/even cases
    if n%2==0:
        Y_hat = np.roll( Y_hat, K - 1 )                              # Rotate
        for x in x_intp:
            y_intp.append(sum([get_summand(c_j,l,x) for l,c_j in zip(range(-K+1,K+1), Y_hat)]))
    else:
        Y_hat = np.roll( Y_hat, K )                                  # Rotate
        for x in x_intp:
            y_intp.append(sum([get_summand(c_j,l,x) for l,c_j in zip(range(-K,K+1), Y_hat)]))
    return x_intp, y_intp


x_range = n * ( x[1] - x[0] )                                        # Fix the range: it ISN'T defined by last element of x
x_intp, y_intp = trig_interpolation( np.fft.fft(y) / n, x_range)     # Divide the FFT by n

plt.plot(x_intp,np.real(y_intp))
plt.plot(x,y, 'o')
plt.show()