我一直在阅读这方面的内容，并且有多种方法可以使用不同大小的 N。我的 Matlab 生锈了，所以它们在 Python 中（N是输入信号的长度x，k是arange(N)= ):$[0, 1, 2, ..., N-1]$

类型 2 DCT 使用 4N FFT 且无移位

信号[a, b, c, d]变成

[0, a, 0, b, 0, c, 0, d, 0, d, 0, c, 0, b, 0, a].

然后取 FFT 得到频谱

[A, B, C, D, 0, -D, -C, -B, -A, -B, -C, -D, 0, D, C, B]

然后扔掉除了第一个之外的所有东西[A, B, C, D]，你就完成了：

u = zeros(4 * N)
u[1:2*N:2] = x
u[2*N+1::2] = x[::-1]

U = fft(u)[:N]
return U.real

使用 2N FFT 镜像的 2 型 DCT (Makhoul)

[a, b, c, d]变成[a, b, c, d, d, c, b, a]. 取它的 FFT 得到 [A, B, C, D, 0, D*, C*, B*]，然后扔掉所有东西，但[A, B, C, D]乘以（半样本移位）得到 DCT：$e^{-j \pi k \over 2 N}$

y = empty(2*N)
y[:N] = x
y[N:] = x[::-1]

Y = fft(y)[:N]

Y *= exp(-1j*pi*k/(2*N))
return Y.real

使用 2N FFT 填充的 2 型 DCT (Makhoul)

[a, b, c, d]变成[a, b, c, d, 0, 0, 0, 0]. 取 FFT 得到 [A, B, C, D, E, D*, C*, B*]，然后扔掉所有东西，但[A, B, C, D]乘以得到 DCT：$2e^{-j \pi k \over 2 N}$

y = zeros(2*N)
y[:N] = x

Y = fft(y)[:N]

Y *= 2 * exp(-1j*pi*k/(2*N))
return Y.real

使用 N FFT (Makhoul) 的 2 型 DCT

信号[a, b, c, d, e, f]变为[a, c, e, f, d, b]，然后用FFT得到[A, B, C, D, C*, B*]。乘以然后取实部：$2 e^{-j \pi k \over 2N}$

v = empty_like(x)
v[:(N-1)//2+1] = x[::2]

if N % 2: # odd length
    v[(N-1)//2+1:] = x[-2::-2]
else: # even length
    v[(N-1)//2+1:] = x[::-2]

V = fft(v)

V *= 2 * exp(-1j*pi*k/(2*N))
return V.real

在我的机器上，这些速度都大致相同，因为生成exp(-1j*pi*k/(2*N))比 FFT 需要更长的时间。:D

In [99]: timeit dct2_4nfft(a)
10 loops, best of 3: 23.6 ms per loop

In [100]: timeit dct2_2nfft_1(a)
10 loops, best of 3: 20.1 ms per loop

In [101]: timeit dct2_2nfft_2(a)
10 loops, best of 3: 20.8 ms per loop

In [102]: timeit dct2_nfft(a)
100 loops, best of 3: 16.4 ms per loop

In [103]: timeit scipy.fftpack.dct(a, 2)
100 loops, best of 3: 3 ms per loop

我不确定您引用的论文使用什么方法，但是我之前已经得出了这个，这就是我最终得到的结果：（假设是输入）$x(n)$

让

$y(n) = \Bigl\lbrace \begin{array}{ll} x(n), & n=0,1,...,N-1\\ x(2N - 1 - n), & n=N,N+1,...,2N-1 \end{array}$

然后 DCT 由下式给出

$C(k) = \mathrm{Re}\Bigl\lbrace e^{-j\pi \frac{k}{2N}} FFT\lbrace y(n)\rbrace\Bigr\rbrace$

因此，您基本上创建了长度序列，其中前半部分是，后半部分是反转。然后只需采用 FFT 并将该向量乘以相位斜坡。最后只取实部，你就有了 DCT。$2N$$y(n)$$x(n)$$x(n)$

这是MATLAB中的代码。

function C = fdct(x)
    N = length(x);
    y = zeros(1,2*N);
    y(1:N) = x;
    y(N+1:2*N) = fliplr(x);
    Y = fft(y);
    k=0:N-1;
    C = real(exp(-j.* pi.*k./(2*N)).*Y(1:N));

编辑：

注意：这里使用的 DCT 公式是：

$C(k) = 2\sum_{n=0}^{N-1} x(n) \cos\left(\frac{\pi k}{2N}(2n+1) \right)$

有几种缩放总和的方法，因此它可能与其他实现不完全匹配。例如，MATLAB 使用：

$C(k) = w(k)\sum_{n=0}^{N-1} x(n) \cos\left(\frac{\pi k}{2N}(2n+1) \right)$

其中和$w(0) = \sqrt\frac{1}{N}$$w(1...N-1) = \sqrt{\frac{2}{N}}$

您可以通过适当缩放输出来解决此问题。

对于真正的科学计算，内存使用量也很重要。因此 N 点 FFT 对我来说更有吸引力。这仅是由于信号的 Hermitian 对称性才可能实现的。这里给出了参考 Makhoul。并且实际上有计算DCT和IDCT或DCT10和DCT01的算法。
http://ieeexplore.ieee.org/abstract/document/1163351/

我为 DCT 和 DST 编写了 python numpy/scipy 函数，如果你想将它与具有可微 FFT 但还没有 DCT 和 DST 的 autodiff 库（如 pytorch）一起使用，这将很有用。

import numpy as np
import scipy.fftpack

def dct_with_fft(x):
    N = x.size
    v = np.empty_like(x)
    v[:(N-1)//2+1] = x[::2]

    if N % 2: # odd length
        v[(N-1)//2+1:] = x[-2::-2]
    else: # even length
        v[(N-1)//2+1:] = x[::-2]

    V = scipy.fftpack.fft(v)

    k = np.arange(N)
    V *= 2 * np.exp(-1j*np.pi*k/(2*N))
    return V.real


def dst_with_fft(x):
    N = x.size
    v = np.empty_like(x)
    v[:(N-1)//2+1] = x[::2]

    if N % 2: # odd length
        v[(N-1)//2+1:] = -x[-2::-2]
    else: # even length
        v[(N-1)//2+1:] = -x[::-2]

    V = scipy.fftpack.fft(1j*v)

    k = np.arange(N)
    V *= 2 * np.exp(-1j*np.pi*k/(2*N))
    return V.imag[::-1]

N = 5
x = np.random.rand(N)

print(dct_with_fft(x))
print(scipy.fft.dct(x))

print(dst_with_fft(x))
print(scipy.fft.dst(x))