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。