How to calculate spectral derivative using scipy.fftpack (DST, DCT)

Question

scipy.fftpack package provides a large number of routines related to Discrete Fourier transforms. I need to calculate 1st and 2nd derivative of a function using DST (DCT) only. However the package contains diff routine which returns kth derivate using FFT

Does anyone know how to get 1st and 2nd derivatives using DST for example? Here is my draft:

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) applies the type II discrete cosine transform to the signal u. The DCT type III is the inverse of the DCT II. Appart for a scaling, it writes:

Hence, the signal u of period T even around zero is written as a weighted sum of cosines of period T/j. Values u[i] are sampled at coordinates x=(2i+1)T/4N. Hence, the signal could be described as:

Its derivative is:

Once sampled at the same point, it looks almost like a type III DST:

Nevertheless shift of sum index is required to match the conventions of scipy dst III:

Here goes a sample making use of DCT II and IDST II to compute the odd-odd derivative of an even-even signal:

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()

The input signal is a ramp, which unfolds to a triangular periodic signal (even-even symetry). Hence, the computed derivative is almost flat, and it unfolds to a square wave.

The high frequencies in the triangular signal are aliased to frequencies below N. Computing the derivative amplifies the high frequencies. As a result, the estimated derivative is plaged by some high freqency oscilations.

If the input signal is corrupted by a significant white noise, the DFT/DST/DCT derivative features a huge high frequency noise. It is even worse for higher derivatives, as it amplifies high frequencies again. To avoid such an unwanted behavior, the derivative, just like the ramp filter, can be combined to a low-pass filter to attenuate that noise. It is a widely used technic as the filtered backprojection algorithm is applied for tomographic reconstruction. See tests of different filters here.

Answer 2

The procedure shown above is essentially right. Just a few remarks:

1. Use Gauss grid theta(j) = (nn-j+0.5)*pi/nn, and a mapping: for example r = L cot(theta/2.0) for a radial semi-infinite lattice.
2. The derivative follows the chain rule, so the the derivative of your mapping dr/dtheta must be appended.
3. This is important: all of this only works when computing derivatives of functions expandable in Chebyshev polynomials. The error gets significantly larger if the function, for example, is dominated by non-integer powers, or is irregular in its domain. Otherwise, it is relatively simple to get machine precision for first derivatives and (very close!) for second derivatives.