Дискретное косинусное преобразование для дифференцирования реальной симметричной функции

Question

Я хотел бы дифференцировать реальную периодическую функцию на (0,2*pi), которая также симметрична относительно x=pi, используя дискретное преобразование Фурье. Я написал код Python, который делает это с использованием FFT/IFFT, но это не учитывает симметрию функции и поэтому является немного расточительным.

(Общая цель состоит в том, чтобы сделать решение псевдоспектрального потока жидкости, и периодичность и симметрия в одном из направлений должны позволить мне расширить переменные в этом направлении, используя только косинусную часть ряда Фурье)

Я знаю, что для этого мне нужно использовать дискретное косинусное преобразование (DCT), но я не могу понять, что нужно изменить в моем домене (x), векторе волновых чисел (k) и реализации сохранения DCT/IDCT, что должны сделать первые два быть на половину длины.

Любая помощь высоко ценится. "Долгое время читатель, первый постер", поэтому заранее извиняюсь за плохое форматирование!

import sympy as sp
import numpy as np
import matplotlib.pylab as plt
from scipy.fftpack import fft, ifft

# Number of grid points
N = 2**5 

# Test function to check results with (using SymPy)
w = 3.; X = sp.Symbol('x'); Y=sp.cos(w*X)

# Domain of regularly spaced points in [0,2pi)
x=(2*np.pi/N)*np.arange(0,N)

# Calc exact derivatives using SymPy then turn into functions
dY      = Y.diff(X)
d2Y     = dY.diff(X)
d3Y     = d2Y.diff(X)
f       = sp.lambdify(X, Y,'numpy')
df_ex   = sp.lambdify(X, dY, 'numpy')
d2f_ex  = sp.lambdify(X, d2Y, 'numpy')
d3f_ex  = sp.lambdify(X, d3Y, 'numpy')

# Wavenumber vector
k=np.hstack(( np.arange(0,N/2), 0, np.arange(-N/2+1,0) )); 
k2=k**2; k3=k**3;

# Trans. to Fourier domain, diff, then return to phyical space
F   = fft(f(x))
df  = np.real(ifft(1j*k*F))
d2f = np.real(ifft( -k2*F))
d3f = np.real(ifft(-1j*k3*F))

# Plot result
fh=plt.figure(figsize=(8,4)); ah=fh.add_subplot(111)
plt.plot(x,f(x),'b-',x,df_ex(x), 'r-',x,d2f_ex(x),'g-',x,d3f_ex(x),'k-')
plt.plot(x,df,'ro',x,d2f,'go',x,d3f,'ko')
plt.xlim([0,2*np.pi])