FFT代码：
def fft(x):
if x.size ==1:
return x
N = x.size
x0 = my_fft(x[0:N:2])
x1 = my_fft(x[0+1:N:2])
k = numpy.arange(N/2)
e = numpy.exp(-2j*numpy.pi*k/N)
l = x0 + x1 * e
r = x0 - x1 * e
return numpy.hstack([l,r])
DCT代码：
def dct(x):
k = 0
N = x.size
xk = numpy.zeros(N)
for k in range(N):
for n in range(N):
xn = x[n]
xk[k] += xn*numpy.cos(numpy.pi/N*(n+1/2.0)*k)
return xk
FCT试用版：
def my_fct(x):
if x.size ==1:
return x
N = x.size
x0 = my_fct(x[0:N:2]) # have to be set to zero?
x1 = my_fct(x[0+1:N:2])
k = numpy.arange(N/2)
n = # ???
c = numpy.cos(numpy.pi/N*(n+1/2.0)*k)
l = x0 #???
r = x0 #???
return numpy.hstack([l,r])
* J。 Makhoul，“一维和二维快速余弦变换”，IEEE Trans。 co语音信号进程28（1），27-34（1980）。
#1 楼
我一直在阅读有关此内容的信息，有多种方法可以使用不同的大小N来实现。我的Matlab生锈了，所以它们在Python中使用（N是输入信号的长度x，k是arange(N) = $ [0，1 ，2，...，N-1] $）：使用2N FFT且无移位的类型2 DCT
信号
[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]以外的所有内容，然后将其乘以$ e ^ {-j \ pi k \ over 2 N} $（半采样移位）以得到DCT：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]以外的所有内容，然后乘以$ 2e ^ {-j \ pi k \ over 2 N} $来获得DCT：y = zeros(2*N)
y[:N] = x
Y = fft(y)[:N]
Y *= 2 * exp(-1j*pi*k/(2*N))
return Y.real
使用N FFT的类型2 DCT（Makhoul）
信号
[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
评论
$ \ begingroup $
很好的答案，对我的实施有很大帮助！附加说明：如果信号长度为奇数，则最后一种方法“使用N FFT的类型2 DCT”仍然可以正常工作；最后一个元素移到中间元素。我已经验证了这一事实的数学和代码。
$ \ endgroup $
– Nayuki
17年7月4日在23:55
$ \ begingroup $
@Nayuki您正在生成exp（-1j * pi * k /（2 * N））还是该步骤的快捷方式？
$ \ endgroup $
– Endolith
17年7月6日在17:43
$ \ begingroup $
我正在我的代码中生成exp（-1j * pi * k /（2 * N）），因为要使DCT到DFT映射有效，必须进行四分之一样本移位。是什么让你问？
$ \ endgroup $
– Nayuki
17年7月7日在6:36
$ \ begingroup $
嗨，这对于Type III DCT如何计算DCT-II的逆数？
$ \ endgroup $
–杰克H
18年1月25日在20:55
$ \ begingroup $
使用2N FFT镜像（Makhoul）的Type 2 DCT是实现2D的方法。
$ \ endgroup $
–爱德华多·里斯（Eduardo Reis）
20/09/11在21:53
#2 楼
我不确定您所引用的论文使用的是哪种方法，但是我之前已经推导了该方法，这就是我最终得到的结果（假设输入$ x（n）$）let
$ 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 $
因此，您基本上创建了$ 2N $长度序列$ y（n）$，其中上半部分是$ x（n）$，下半部分是$ x（n）$相反。然后只需进行FFT，然后将该矢量乘以一个相位斜坡即可。最后只需要真正的部分，便拥有了DCT。这是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));
编辑： br />注意：使用的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}} $
您可以通过适当缩放输出来解决这个问题。
评论
$ \ begingroup $
y（n）应该是N长度，而不是2N长度。通过从N长度信号中计算N长度DCT而不是从2N信号中进行2N FFT，您将获得4倍的计算速度。 fourier.eng.hmc.edu/e161/lectures/dct/node2.html www-ee.uta.edu/dip/Courses/EE5355/Discrete%20class%201.pdf
$ \ endgroup $
– Endolith
13年8月29日在18:57
#3 楼
对于真正的科学计算，内存使用量也很重要。因此N点FFT对我来说更有吸引力。这仅由于信号的埃尔米特对称性才可能实现。
此处提供了参考Makhoul。并且实际上具有用于计算DCT和IDCT或DCT10和DCT01的算法。http：//ieeexplore.ieee.org/abstract/document/1163351/
评论
您是在问您的DCT代码是否正确吗？谢谢您的意见。我在开头添加了另一句话。我的目标是在FFT的基础上实现FCT。
使用CTT或IFFT在逆DCT（IDCT）计算中给出IDCT。