8

使用 Rprof 发现 dtt 包中的 dct 是一段运行非常缓慢的 R 代码中的主要违规者。将它换成 stats 包中的 fft（这不是相同的转换，但应该花费相同的时间来计算）我的运行时间显着改善。事实上，我的 Rprof 行中有 2/3 之前是 dct 调用，而在进行切换后，大约 600 行中只有 3 行是 fft 调用。

dtt 包中的 dct 实现不是使用快速离散傅立叶变换完成的吗？如果是这样，是否有一个包有一个？（我知道可以将数据加倍，然后从这些 fft 系数中提取 dct 的系数，但是直接快速 dct 肯定会更好，而且真的应该在某个地方有一个）。

4

3 回答 3

12

似乎没有快速 dct，但在 stats 包中有一个 fft（快速傅立叶变换），所以这里是如何使用 fft 获得快速 dct。

使用它需要您自担风险。我没有对它进行任何认真的检查。我在几个不同大小的向量上检查了它，它给出了与 dtt 包中的函数 dct 相同的结果。如果有人想通过将其与 dct 的输出进行比较来仔细检查我，那么请随时这样做并发布您的结果。

将向量扩展为向量的两倍，如下所示：如果向量是 v=(1,2,3)，则将条目加倍为 w=(1,2,3,3,2,1)。注意排序。如果您的向量是 v=(1,2,4,9)，则将条目加倍为 w=(1,2,4,9,9,4,2,1)

让 N 为您的原始向量的长度（在您将其长度加倍之前）。

那么 .5 * fft(w)/exp(complex(imaginary=pi / 2 / N)*(seq(2*N)-1)) 的前 N ​​个系数应该与计算 dct(v) 一致，除非它应该是在几乎所有情况下都显着加快。

速度考虑。如果你是素数 N，那么计算快速 dct 所需的时间就像为每个素因数做一个慢速 dct 所需的时间。因此，如果 N 为 2^K，它就像对长度为 2 的向量进行 K 个不同的慢速 dct 变换，所以它真的很快。如果 N 是素数（最坏的情况），则根本没有加速。最大的加速是在长度为 2 的幂的向量上。

注意：上面的 R 代码看起来非常不友好，所以让我说一下是怎么回事。以正确的方式将长度加倍后，你得到的 fft 的前 N ​​个系数几乎是正确的。然而，系数需要稍微调整一下。让 P 代表 e^(pi * i / 2 / N)。不理会第一个系数。将第二个系数除以 P，将第三个系数除以 P^2，将第四个系数除以 P^3，等等......然后将所有系数除以 2（包括第一个系数）以与 R 用于 dct 的归一化一致.

这应该与使用 dtt 包中的 dct 函数相同，但在几乎所有情况下都快得多。

于 2012-06-28T18:57:11.407 回答
2

约翰和伊戈尔的回答都是正确的，但我认为他们缺少一些示例代码。

library(dtt)

par(mfrow=c(3, 1), mar=c(2, 3, 1, 0.1), mgp=c(2, 0.8, 0))

set.seed(1)
N <- 60
v <- sin(seq(0, pi*2*4, l=N))/2 + 
     sin(seq(0, pi*2*7, l=N))/3 + 
     sin(seq(0, pi*2*15, l=N))/2 + 
     runif(N, -1, 1)/40 +
     runif(N, -1, 1)/40

plot(v, type="o")

DCT-I：

v.dct1 <- dct(v, variant=1)

w <- c(v, v[(N - 1):2])
w.dct1 <- 0.5 * Re(fft(w)[1:N])


plot(v.dct1, type="l", col="#00000088")
abline(h=0, col="#00000044")
lines(w.dct1, col=2, lty=3, lwd=2)
legend("topright", c("dtt::dct", "fft"), bty="n", col=1:2, lwd=1)
legend("topleft", "DCT-I", bty="n")

DCT-II：

v.dct2 <- dct(v, variant=2)

P <- exp(complex(imaginary=pi / 2 / N)*(seq(2*N)-1)) 

w <- c(v, v[N:1])
w.dct2 <- 0.5 * Re(fft(w)[1:N]/P)

plot(v.dct2, type="l", col="#00000088")
abline(h=0, col="#00000044")
lines(w.dct2, col=2, lty=3, lwd=2)
legend("topright", c("dtt::dct", "fft"), bty="n", col=1:2, lwd=1)
legend("topleft", "DCT-II", bty="n")

在此处输入图像描述

于 2018-08-16T19:07:45.870 回答
2

免责声明：我从未使用过该dtt软件包，也无法将我的结果与其结果进行比较。我的回答是关于 DCT/FFT 的通用答案。

John 的想法是正确的，但他在重复向量方面是两不误的，并且必须通过调整系数（e^(pi * i / 2 / N)因子）来补偿它。通过原始向量的正确扩展，FFT 可以直接产生正确的结果 (*)。引用维基百科：

" DCT-I 与 2N-2 个具有偶数对称性的实数的 DFT 完全等价（总比例因子为 2）。例如，N=5 个实数 abcde 的 DCT-I 完全等价于八个实数 abcdedcb（偶数对称）除以 2 的 DFT。 ”

即，如果我们有一个v = [1, 2, 3, 4, 5]希望对其执行 DCT 的向量，我们应该构造一个新向量w = [1, 2, 3, 4, 5, 4, 3, 2]并对其执行 FFT。请注意， 的第一个和最后一个组件v仅按原始顺序出现一次w！

这是因为偶函数的傅里叶变换（函数关于零对称）纯粹由实数（余弦）系数组成。如果我们像约翰建议的那样w通过包含整个 reversed 来构造向量v，它将围绕 对称-0.5。由于这种微小的偏移，傅立叶变换也会产生虚（正弦）系数。

(*) 这里的方法产生 DCT-I。John 的方法似乎针对的是 DCT-II。

于 2018-04-05T10:33:14.337 回答