如何利用scipy.fftpack (DST，DCT)计算光谱导数

Question

scipy.fftpack软件包提供了大量与离散傅里叶变换相关的例程。我只需要使用DST (DCT)计算函数的一阶和二阶导数。但是，该包包含diff例程，它使用FFT返回kth派生函数。

有人知道如何使用DST获得第一和第二衍生产品吗？这是我的草稿：

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import dst, idst, dct, idct

L   = 10
N   = 100
a   = 0.4
x   = np.linspace(0,L,N)
# function
u   = np.sin(2*np.pi*x/L)*np.exp(-a*x)

# exact 1st derivative
du  = np.exp(-a*x)*(-a*np.sin(2*np.pi*x/L) + np.cos(2*np.pi*x/L)*2*np.pi/L)

# get 1st derivative
dufft = idst(-dct(u))

plt.figure()
plt.plot(x,u)
plt.plot(x,du)
plt.figure()
plt.plot(x,dufft)
plt.show()

​

​

Answer 1

dct(u)将Ⅱ型离散余弦变换应用于信号u。DCT类型III与DCT II相反。在缩放时，它写道：

​

​

因此，周期T甚至在零附近的信号u被写成周期T/j余弦的加权和，在坐标x=(2i+1)T/4N处采样ui值。因此，该信号可描述为：

​

​

它的衍生物是：

​

​

一旦在同一点取样，它看起来几乎像第三类DST：

​

​

然而，为了与dst III的惯例相匹配，需要改变和指数。

​

​

下面是一个使用DCTⅡ和IDSTⅡ来计算偶数信号的奇数导数的示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import dst, idst, dct, idct


N=100
x=np.linspace(0,100,N)
u=np.linspace(0,100,N)

hatu=dct(u,type=2)
#multiply by frequency, minus sign
for i in range(N):
    hatu[i]=-(i)*hatu[i]
#shift to left
hatu[0:N-1]=hatu[1:N]
hatu[N-1]=0.

#dst type III, or equivalently IDST type II.
dotu=idst(hatu,type=2)

dotu=dotu/(2*N)


plt.figure()
plt.plot(x,u,label='the function')
plt.plot(x,dotu,label='its derivative')
plt.legend()
#plt.figure()
#plt.plot(x,dufft)
plt.show()

输入信号是一个斜坡，向三角形周期信号展开(偶数-偶对称).因此，计算得到的导数几乎是平坦的，并向方波展开。

​

​

三角信号中的高频被化成N以下的频率，计算导数放大高频。因此，估计的导数被一些高频率的振荡所淹没。

当输入信号受到明显的白噪声干扰时，DFT/DST/DCT导数具有巨大的高频噪声。对于更高的衍生品来说，情况就更糟了，因为它再次放大了高频。为了避免这种不必要的行为，导数，就像斜坡滤波器，可以组合成一个低通滤波器，以减少噪音。它是一种广泛应用的滤波反投影算法在层析重建中的应用技术。参见不同过滤器的测试这里。

Answer 2

上述程序基本上是正确的。只想说几句：

1. 使用Gauss网格θ(J)=(nn+0.5)*pi/nn和映射:例如r=L(theta/2.0)对径向半无限格。
2. 导数遵循链规则，所以映射dr/dtheta的导数必须附加。
3. 这一点很重要:所有这些只有在计算在Chebyshev多项式中可扩展的函数的导数时才能工作。例如，如果函数被非整数幂控制，或者在其域中不规则，则误差会显著增大。否则，获得第一导数和(非常接近！)的机器精度是相对简单的！第二种衍生产品。