r - R radix-2 DIT 案例中的 Cooley-Tukey FFT

Question

所以我一直在尝试(手动)在 R 中实现 Cooley-Turkey FFT 算法(对于大小为 N=n^2 的输入)。我试过:

myfft <- function(s){
  N <- length(s)
  if (N != 1){
    s[1:(N/2)] <- myfft(s[(1:(N/2))*2-1])
    s[(N/2+1):N] <- myfft(s[(1:(N/2))*2])

    for (k in 1:(N/2)){
      t <- s[k]
      s[k] <- t + exp(-1i*2*pi*(k-1)/N) * s[k+N/2]
      s[k+N/2] <- t - exp(-1i*2*pi*(k-1)/N) * s[k+N/2]
    }
  }
  s
}

这可以编译，但对于 n>1，N=2^n，它不会计算出正确的值。我实现了一个 DFT 函数并使用 fft() 函数进行比较，两者都计算，当归一化时，给出相同的值，但似乎不同意我上面的算法。

如果有人感兴趣并看到我哪里出错了，将不胜感激，我会疯狂地寻找错误并开始质疑，如果我曾经理解过这个 FFT 算法。

更新:我修复了它，我不是 100% 确定问题到底出在哪里，但这是有效的实现:

myfft <- function(s){
  N <- length(s)
  if (N != 1){
    t <- s
    t[1:(N/2)] <- myfft(s[(1:(N/2))*2-1]) # 1 3 5 7 ... 
    t[(N/2+1):N] <- myfft(s[(1:(N/2))*2]) # 2 4 6 8 ... 

    s[1:(N/2)] <- t[1:(N/2)] + exp(-1i*2*pi*(0:(N/2-1))/N) * t[(N/2+1):N]
    s[(N/2+1):N] <- t[1:(N/2)] - exp(-1i*2*pi*(0:(N/2-1))/N) * t[(N/2+1):N]

  }
  return(s)
}

Answer 1

问题出在下面这行

s[1:(N/2)] <- myfft(s[(1:(N/2))*2-1])

它覆盖了后续行所需的部分未转换值:

s[(N/2+1):N] <- myfft(s[(1:(N/2))*2])

例如，当N=4时，第二次调用myfft使用s[2]和s[4] ，但是第一次调用 myfft 的赋值写入了 s[1] 和 s[2](因此覆盖了s[2] 中需要原始值)。

您复制整个数组的解决方案可防止这种覆盖。

常用的替代解决方案是分别复制偶数和奇数部分:

myfft <- function(s){
  N <- length(s)
  if (N != 1){
    odd <- s[(1:(N/2))*2-1]
    even <- s[(1:(N/2))*2]
    s[1:(N/2)] <- myfft(odd)
    s[(N/2+1):N] <- myfft(even)

    s[1:(N/2)] <- t[1:(N/2)] + exp(-1i*2*pi*(0:(N/2-1))/N) * t[(N/2+1):N]
    s[(N/2+1):N] <- t[1:(N/2)] - exp(-1i*2*pi*(0:(N/2-1))/N) * t[(N/2+1):N]
  }
  return(s)
}