R で離散フーリエ変換を行う方法は以前にも確認したが、DFT には 1/N をつける位置などにいろいろな流儀があってややこしい。
一歩一歩内容を確認しながら検証した。
一次元
N <- 128 x <- seq(from=0,length.out=N, by=1/N) v <- 2 + 3 * sin(x * 2 * 2 * pi) - 2 * cos(x * 8 * 2 * pi) plot(x, v, type='l') rv <- fft(v) / N plot(Mod(rv), type='h') plot(Arg(rv) / pi * 180 * (Mod(rv) > 0.01), type='h') # ある程度の絶対値があるやつだけ rv[1:10] rv[N:(N-9)] # この状態(fft 結果を N で割った状態)で # rv[1] が DC (縮退してるから定数項そのもの) # ここから下は、負の周波数とに等分されるので、上の係数の半分となる # rv[2] は 周波数 1 の 成分 (rv[2] == Conj(rv[128]), rv[3] == Conj(rv[127]), ...) # rv[3] は 周波数 2 の 成分... # 位相については、+sin が -90度、-sin が 90 度。+cos が 0度。-cos が 180度。 # 後半は折り返し(共役) all.equal(rv[2:(N / 2 + 1)], Conj(rv[N:(N / 2 + 1)])) # Nyquist 周波数 N / 2 に該当する成分(1 + N/2 = 65)は、 # 自分自身と重なっているので、ここも実数になる # 実数ということは、Nyquist 周波数の波については、位相情報が失われているということらしい。 rv[N / 2 + 1] # 元に戻す v2 <- fft(rv, inverse=TRUE) all.equal(v, Re(v2)) all.equal(0, sum(abs(Im(v2))))
二次元の場合
N <- 128 x <- seq(from=0,length.out=N, by=1/N) x_2d <- outer(x, x, function(x, y) { sin((4 * x + 3 * y) * 2 * pi) + 1.2 * sin((6 * x + 7 * y) * 2 * pi)}) image(x_2d, asp=1) image(Mod(x_2d), asp=1) # 絶対値なので縞が増える fft1 <- t(apply(x_2d, 1, fft)) / N # 1 の方向に apply したら t が必要 image(Mod(fft1), asp=1) image(Arg(fft1) * (Mod(fft1) > 0.1), asp=1) fft2 <- apply(fft1, 2, fft) / N image(Mod(fft2), asp=1) image(Arg(fft2) * (Mod(fft2) > 0.1), asp=1) ifft1 <- apply(fft2, 2, fft, inverse=T) image(Mod(ifft1), asp=1) image(Arg(ifft1) * (Mod(ifft1) > 0.1), asp=1) ifft2 <- t(apply(ifft1, 1, fft, inverse=T)) image(Mod(ifft2), asp=1) image(Arg(ifft2) * (Mod(ifft2) > 0.1), asp=1) # 確認 image(x_2d, asp=1) image(Re(ifft2), asp=1) all.equal(x_2d, Re(ifft2)) all.equal(0, sum(abs(Im(ifft2)))) # 順序は意味があるのか? 逆の順序で戻してみる ifft1 <- t(apply(fft2, 1, fft, inverse=T)) image(Mod(ifft1), asp=1) image(Arg(ifft1) * (Mod(ifft1) > 0.1), asp=1) ifft2 <- apply(ifft1, 2, fft, inverse=T) persp3d(x, x, Mod(ifft2), col='red') image(Mod(ifft2), asp=1) image(Arg(ifft2) * (Mod(ifft2) > 0.1), asp=1) # 逆順でお同じ結果になることを確認 image(x_2d, asp=1) image(Re(ifft2), asp=1) all.equal(x_2d, Re(ifft2)) all.equal(0, sum(abs(Im(ifft2)))) # 実は一気にやれる fft3 <- fft(x_2d) / N / N image(Mod(fft3), asp=1) image(Arg(fft3) * (Mod(fft3) > 0.1), asp=1) all.equal(fft3, fft2) # mvfft は、column 方向に FFT するのと同じ fftmv <- mvfft(x_2d) / N image(Mod(fftmv), asp=1) fft4 <- apply(x_2d, 2, fft) / N image(Mod(fft4), asp=1) all.equal(fft4, fftmv)
