Use DCT instead of RFFT/IRFFT in pyfftw_sdp

[NimaSarajpoor]

MASS_V4 uses Discrete-Cosine-Transform(DCT) for sliding-dot-product. See issue #28 for more details. We started to implement this in python and check its performance in PR #30. It might be worth it to try DCT in "pyfftw_sdp". We need DCT (type II) and DCT (type I). According to pyfftw documentation and FFTW documentation (section 4.8.3 DCTs), we can use:

• pyfftw.FFTW(..., direction=FFTW_REDFT00) for DCT(type I)
• pyfftw.FFTW(..., direction='FFTW_REDFT10') for DCT(type II)

[NimaSarajpoor]

According to Scipy documentation:

For a single dimension array x, dct(x, norm='ortho') is equal to MATLAB dct(x)

However, pyFFTW's DCT is equivalent to scipy.fft.dct(x, norm=None), and it does not provide support for the parameter "norm". So we need to modify the output . See the notes in scipy documentation.

DCT Type I

import numpy as np
import pyfftw
import scipy

# input/output size
n = 25  

# SCIPY
a = np.random.rand(n)
ref = scipy.fft.dct(a, type=1, norm='ortho')

# PYFFTW
a_input = pyfftw.empty_aligned(n, dtype=np.float64)
dct_type_I = pyfftw.builders.dct(a_input, type=1)
a_input[:] = a
comp = dct_type_I()

# MODIFY PYFFTW's OUTPUT
comp += (np.sqrt(2) - 1) * (a[0] + ((-1) ** np.arange(n)) * a[-1])
comp[:] = comp / np.sqrt(2 * (n - 1))
comp[0] = comp[0] / np.sqrt(2)
comp[-1] = comp[-1] / np.sqrt(2)

# ASSERTION
np.testing.assert_allclose(comp, ref)

DCT Type 2

import numpy as np
import pyfftw
import scipy

# input/output size
n = 25  

# SCIPY
a = np.random.rand(n)
ref = scipy.fft.dct(a, type=2, norm='ortho')

# PYFFTW
a_input = pyfftw.empty_aligned(n, dtype=np.float64)
dct_type_II = pyfftw.builders.dct(a_input, type=2)
a_input[:] = a
comp = dct_type_II()

# MODIFY PYFFTW's OUTPUT
comp[0] *= 1 / (2 * np.sqrt(n))
comp[1:] *= 1 / np.sqrt(2 * n)

# ASSERTION
np.testing.assert_allclose(comp, ref)