FFTWから各周波数のすべての値を加算する

Question

FFTW（fftw_plan_dft_r2c_3d）を使用して3Dフーリエ変換を実行しました。繰り返し周波数を含むすべての周波数でトランスフォームの（log of）値を合計したいと思います実際に出力配列に格納されていない（私はサイズがNx x Ny x (Nz/2 + 1)であると理解しています）。二重計算をしないでこれを行うにはどうすればよいですか？

Answer 1

大きな質問です。申し訳ありませんが、私の答えは少し長めですが、私は間違いをしないようにしたいと思います。ここでgo-

複素数と複素数の3D FFTの合計ログの大きさは、複数の場合、実数から複素数の3次元FFTの対数の総和に等しくなります前者から失われた後者の（最後の次元の）すべての「スライス」。

• Nzが偶数の場合、それは最初と最後のスライス以外のすべてのスライスをダブルカウントを意味します。
• Nzが奇数の場合は、最初のスライスを除くすべてのスライスをダブルカウントします。

（これは、πのラジアンの角周波数（フェーザー-1に対応）が含まれていますが、奇数長のものはそれを停止します。このパターンを覚えているので、N/4対N = 3のフェーザーを単位円の周りに描画して、奇数かπ-radかを思い出させる。

Numpy/Pythonを使用したアイデアの実験的検証は、実数から複素数のFFTの表記法はFFTWと一致すると信じています。Nx = 10をNy = 20でNz = 8の実配列で生成します。複素数から複雑な3D FFT（Nx複雑な配列でNyをNxに変換する）とを実数から複素数の3D FFT（Nxを（Nz/2 + 1）の複素数配列で生成する）を計算します。 Nzが偶数であるため、最初の&最後のスライス以外のすべてをダブルカウントすると、前者のログ総和の大きさはと同じであることを確認してください。

コード：あなたは奇妙 Nzでこれを再実行した場合

import numpy as np 
import numpy.fft as fft 

Nx = 10 
Ny = 20 
Nz = 8 

x = np.random.randn(Nx, Ny, Nz) 

Xf = fft.fftn(x) 
Xfr = fft.rfftn(x) 

energyProduct1 = np.log10(np.abs(Xf)).sum() 

lastSlice = -1 if Nz % 2 is 0 else None 
energyProduct2 = np.log10(np.abs(np.dstack((Xfr, Xfr[:, :, 1:lastSlice])))).sum() 

print('Difference: %g' % (energyProduct1 - energyProduct2)) 
# Difference: -4.54747e-13 

は、あなたが表示されますその0の機械精度内の複素数 - 複素数と実数 - 複素数遺跡の違いNzもあるのでnp.dstack((Xfr, Xfr[:, :, 1:lastSlice))（dstack用ドキュメント、fft.rfftn）が第3次元、最後から二番目にスライスを最後から二番目のためにその第二とrfftn出力をスタック、そしてあなたは、ダブルカウント0または-πしないこと。

DFTビン。

もちろん

、これを行うための別の方法は、実数 - 複素数配列上の和のログ・大きさを計算することである、ダブルそれは、その後、減算和のログ-大きさを超えます最初のスライス、最後のスライス（Nzが偶数の場合）。

tl; dr実数から複合出力へのログの大きさを合計します。それを倍増させる。この結果から、（3次元の）最初のスライスの合計ログの大きさを差し引きます。 Nzが奇妙な場合は、完了です。 Nzが偶数の場合は、最後のスライスの合計ログの大きさも減算します。