R fft 快速离散傅立叶变换 (FFT)

说明

使用快速算法 “Fast Fourier Transform” (FFT) 计算数组的离散傅里叶变换 (DFT)。

用法

fft(z, inverse = FALSE)
mvfft(z, inverse = FALSE)

参数

值

当 z 是向量时，fft 计算并返回的值是 z 中的值序列的非归一化单变量离散傅立叶变换。具体来说，y <- fft(z) 返回

对于 h = 1,\ldots,n 其中 n = length(y) 。如果 inverse 为 TRUE ，则 \exp(-2\pi\ldots) 替换为 \exp(2\pi\ldots) 。

当 z 包含数组时，fft 计算并返回多元(空间)变换。如果 inverse 是 TRUE ，则返回(非标准化)傅里叶逆变换，即，如果 y <- fft(z) ，则 z 是 fft(y, inverse = TRUE) / length(y) 。

相比之下，mvfft 采用实数或复数矩阵作为参数，并返回形状相似的矩阵，但每一列都由其离散傅立叶变换替换。这对于分析vector-valued系列很有用。

当要变换的序列的长度高度复合(即具有许多因子)时，FFT 速度最快。如果不是这种情况，则转换可能需要很长时间来计算，并且会使用大量内存。

例子

x <- 1:4
fft(x)
fft(fft(x), inverse = TRUE)/length(x)

## Slow Discrete Fourier Transform (DFT) - e.g., for checking the formula
fft0 <- function(z, inverse=FALSE) {
  n <- length(z)
  if(n == 0) return(z)
  k <- 0:(n-1)
  ff <- (if(inverse) 1 else -1) * 2*pi * 1i * k/n
  vapply(1:n, function(h) sum(z * exp(ff*(h-1))), complex(1))
}

relD <- function(x,y) 2* abs(x - y) / abs(x + y)
n <- 2^8
z <- complex(n, rnorm(n), rnorm(n))
## relative differences in the order of 4*10^{-14} :
summary(relD(fft(z), fft0(z)))
summary(relD(fft(z, inverse=TRUE), fft0(z, inverse=TRUE)))

来源

在 Singleton (1979) 中使用 Fortran 代码的 C 翻译。

参考

Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988). The New S Language. Wadsworth & Brooks/Cole.

Singleton, R. C. (1979). Mixed Radix Fast Fourier Transforms, in Programs for Digital Signal Processing, IEEE Digital Signal Processing Committee eds. IEEE Press.

Cooley, James W., and Tukey, John W. (1965). An algorithm for the machine calculation of complex Fourier series, Mathematics of Computation, 19(90), 297-301. doi:10.2307/2003354.

也可以看看

convolve、nextn。