Cooley-Tukey FFT в корпусе R radix-2 DIT - PullRequest
Новая
       15

Cooley-Tukey FFT в корпусе R radix-2 DIT

2 голосов
/ 16 марта 2019

Итак, я пытался (вручную) реализовать алгоритм БПФ Кули-Турции в R (для входов с размером 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 () для сравнения, оба вычисления, когда нормализованы, дают одинаковые значения, но, похоже, не согласны с моим алгоритмом выше.

Если кто-то почувствует заинтересованность и увидит, где я ошибся, помощь будет принята с благодарностью, я схожу с ума в поисках ошибки и начинаю сомневаться, понял ли я когда-либо этот алгоритм БПФ.

ОБНОВЛЕНИЕ: я исправил это, я не уверен на 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)
}

1 Ответ

0 голосов
/ 17 марта 2019

Проблема была со следующей строкой 

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)
}
...