numpy.fft.rfft2の出力をどのように解釈すればよいですか？

Question

明らかに、rfft2関数は単純に入力行列の離散fftを計算します。しかし、どのように私は出力の特定のインデックスを解釈するのですか？出力のインデックスが与えられていると、フーリエ係数はどのように見えますか？
私は特に出力のサイズが混乱しています。 n×n行列の場合、出力はn×（n/2）+1行列（偶数nの場合）のように見えます。なぜ正方行列は非正方フーリエ変換で終わるのですか？ドキュメントから

Answer 1

numpy.fft.rfft2の出力がありますnumpy.fft.fft2によって計算されるように、標準的な2次元FFTの左半分（プラス1列）を単純に計算します。実際の配列のFFTはnatural and simple symmetryであるため、rfft2は結果の右半分を供給する必要はなく、完全半分のFFTの右半分はその対称性を使用して左半分から導き出すことができます。

例を示します。まず、私はnumpyののランダムな状態と印刷オプションを設定します、それが再現しやすく、見やすいようにするには：

In [1]: import numpy as np 

In [2]: np.set_printoptions(precision=3, suppress=True, linewidth=128) 

In [3]: random = np.random.RandomState(seed=15206) 

のは、6行6列で、実際の入力配列を作成してみましょう：

In [4]: x = random.randn(6, 6) 

In [5]: x 
Out[5]: 
array([[ 1.577, 0.426, 0.322, -0.891, -0.793, 0.017], 
     [ 0.238, 0.603, -0.094, -0.087, -0.936, -1.139], 
     [-0.583, 0.394, 0.323, -1.384, 1.255, 0.457], 
     [-0.186, 0.687, -0.815, -0.54 , 0.762, -0.674], 
     [-1.604, -0.557, 1.933, -1.122, -0.516, -1.51 ], 
     [-1.683, -0.006, -1.648, -0.016, 1.145, 0.809]]) 

（ rfft2を、 fft2を使用していない）

今フルFFTを見てみましょう：ここでの対称性があることを

In [6]: fft2_result = np.fft.fft2(x) 

In [7]: fft2_result 
Out[7]: 
array([[ -5.834+0.j , 1.084-2.33j , -6.504-3.884j, 3.228-0.j , -6.504+3.884j, 1.084+2.33j ], 
     [ 1.475-3.311j, 1.865-3.699j, 2.777-0.095j, -2.570-1.152j, 4.705-3.373j, 4.555-3.657j], 
     [ 2.758+3.339j, -3.512+0.398j, 5.824-4.045j, 1.149-3.705j, 0.661-2.127j, 12.368+1.464j], 
     [ 1.326-0.j , 1.191-4.479j, -3.263+6.19j , 8.939-0.j , -3.263-6.19j , 1.191+4.479j], 
     [ 2.758-3.339j, 12.368-1.464j, 0.661+2.127j, 1.149+3.705j, 5.824+4.045j, -3.512-0.398j], 
     [ 1.475+3.311j, 4.555+3.657j, 4.705+3.373j, -2.570+1.152j, 2.777+0.095j, 1.865+3.699j]]) 

お知らせ：すべての指標iためとjと0 <= i < 6と0 <= j < 6,fft2_result[i, j]は、fft_result[-i, -j]の複素共役である。たとえば、次のように我々ははそれが左半分に由来することができるので、出力の右半分を含めるを必要としないことを意味し

In [8]: fft2_result[2, 4] 
Out[8]: (0.66075993512998199-2.127249005984857j) 

In [9]: fft2_result[-2, -4].conj() 
Out[9]: (0.66075993512998199-2.127249005984857j) 

。フルFFTの左半分を計算するだけで、メモリを節約できます。そして、それはrfft2はまったく同じものです：

In [10]: rfft2_result = np.fft.rfft2(x) 

In [11]: rfft2_result 
Out[11]: 
array([[ -5.834+0.j , 1.084-2.33j , -6.504-3.884j, 3.228+0.j ], 
     [ 1.475-3.311j, 1.865-3.699j, 2.777-0.095j, -2.570-1.152j], 
     [ 2.758+3.339j, -3.512+0.398j, 5.824-4.045j, 1.149-3.705j], 
     [ 1.326-0.j , 1.191-4.479j, -3.263+6.19j , 8.939-0.j ], 
     [ 2.758-3.339j, 12.368-1.464j, 0.661+2.127j, 1.149+3.705j], 
     [ 1.475+3.311j, 4.555+3.657j, 4.705+3.373j, -2.570+1.152j]]) 

お知らせrfft2_resultその試合fft2_result[:, :4]、少なくとも数値誤差まで：

In [12]: np.allclose(rfft2_result, fft2_result[:, :4]) 
Out[12]: True 

我々はまた、左ではなく、出力の上半分を維持することを選択することができますnp.fft.rfft2にaxes引数を使用して、半分、：

として

In [13]: np.fft.rfft2(x, axes=[1, 0]) 
Out[13]: 
array([[ -5.834+0.j , 1.084-2.33j , -6.504-3.884j, 3.228-0.j , -6.504+3.884j, 1.084+2.33j ], 
     [ 1.475-3.311j, 1.865-3.699j, 2.777-0.095j, -2.570-1.152j, 4.705-3.373j, 4.555-3.657j], 
     [ 2.758+3.339j, -3.512+0.398j, 5.824-4.045j, 1.149-3.705j, 0.661-2.127j, 12.368+1.464j], 
     [ 1.326+0.j , 1.191-4.479j, -3.263+6.19j , 8.939-0.j , -3.263-6.19j , 1.191+4.479j]]) 

210の場合np.fft.rfftnの場合、NumPyは指定された最後の軸で実数FFTを実行し、他の軸では複素数FFTを実行します。

もちろん、rfft2_resultにはまだ冗長性があります。最初の列の下半分と最後の列の下半分を破棄し、前と同じ対称性を使用して再構成することができます。 [0, 0]、[0, 3]、[3, 0]、[3, 3]のエントリはすべて実数になるので、虚数部は破棄できます。しかしそれは、あまり便利でない配列表現を私たちに残すでしょう。 DOCによれば

：

Answer 2

：

https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfft2.html

This is really just rfftn with different default behavior. For more details see rfftn.

numpy.fft.rfft2(a, s=None, axes=(-2, -1), norm=None)[source] 

対

numpy.fft.rfftn(a, s=None, axes=None, norm=None)[source] 

https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfftn.html#numpy.fft.rfftn

Notes

The transform for real input is performed over the last transformation axis, as by rfft, then the transform over the remaining axes is performed as by fftn. The order of the output is as for rfft for the final transformation axis, and as for fftn for the remaining transformation axes.

See fft for details, definitions and conventions used.

私は2D fftnで働いていないが、私は1D FFTの解釈のためにこれを作成し、あなたの2D出力の解釈にいくつかの洞察力を与える可能性がある

https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html#numpy.fft.fft

：

import math 
import numpy as np 

PERIOD   = 30 
SIFT   = 2 # integer from 1 to PERIOD/2 

def fourier_series(array, period, sift): 

    # Create an array of length data period; then reverse its order 
    signal   = (array[-period:])[::-1] 

    # Extract amplitude and phase in (a + bi) complex form 
    complex_fft  = np.fft.fft(signal) 

    ''' Calculate amplitude, phase, frequency, and velocity ''' 
    # Define empty lists for later use 
    amplitude  = [] 
    phase   = [] 
    frequency  = [] 
    velocity  = [] 

    # Extract real and imaginary coefficients from complex scipy output 
    for n in range(period, 0, -1): 
     amplitude.append(complex_fft.real[-n]) 
     phase.append(complex_fft.imag[-n]) 

    # The final equation will need to be divided by period 
    # I do it here so that it is calculated once saving cycles  
    amplitude = [(x/period) for x in amplitude]  

    # Extract the carrier 
    carrier = max(amplitude) 

    # The frequency is a helper function of fft 
    # It only has access to the length of the data set 
    frequency.append(np.fft.fftfreq(signal.size, 1)) 

    # Convert frequency array to list 
    frequency  = frequency[-1] 

    # Velocity is 2*pi*frequency; I do this here once to save cpu time  
    velocity  = [x*2*math.pi for x in frequency] 

    ''' Calculate the Full Spectrum Sinusoid ''' 
    # Recombine ALL elements in the form An*sin(2*pi(Fn) + Pn) for full spectrum 
    full_spectrum = 0 
    for m in range(1, period+1): 
     full_spectrum += amplitude[-m]*(1+math.sin(velocity[-m] + phase[-m])) 

    ''' Calculate the Filtered Sinusoid ''' 
    # Normalize user sift input as an integer 
    sift = int(sift) 

    # If sift is more than half of the period, return full spectrum 
    if sift >= period/2: 
     filtered_transform = full_spectrum 

    # If sift is 0 or 1, return the carrier  
    else: 
     filtered_transform = carrier 

     # For every whole number of sift over 1, but less than 0.5*period: 
     # Add an 2 elements to the sinusoid respective of 
     # a negative and positive frequency pair 
     if sift > 1: 
      for m in range(1, sift): 
       p = period - m 
       filtered_transform += amplitude[-m]*(1+math.sin(velocity[-m] + phase[-m])) 
       filtered_transform += amplitude[-p]*(1+math.sin(velocity[-p] + phase[-p])) 

    ''' Print Elements and Return FFT''' 
    if 1: 
     print('**********************************') 
     print('Carrier: %.3f' % amplitude[-period]) 
     print(['%.2f' % x for x in amplitude])  
     print(['%.2f' % x for x in velocity]) 
     print(['%.2f' % x for x in phase]) 

    return filtered_transform, carrier, full_spectrum 

stochastic = # Your 1d input array 
y, y_carrier, y_full = fourier_series(stochastic, PERIOD, SIFT) 

Answer 3

もfft出力における係数の順序を注意デフォルトで第1要素は、0周波数成分（有効和または配列の平均値）の係数であり、そして出発します2番目から私たちは昇順にpostive周波数のためのcoefcientsを持ち、n/2 + 1から降順で負の周波数のために始まります。

np.fft.fftfreq(10) 出力である：出力は、それがこの周波数の順序に対応するようにシフトされるcf=np.fft.fft(array) array([ 0. , 0.1, 0.2, 0.3, 0.4, -0.5, -0.4, -0.3, -0.2, -0.1])

使用np.fft.fftshift(cf)：長さ10のアレイの周波数のビューを有すること を array([-0.5, -0.4, -0.3, -0.2, -0.1, 0. , 0.1, 0.2, 0.3, 0.4]) これはプロットの意味合いを示しています。

2Dの場合は同じです。 fft2とrfft2の違いは、他の人によって説明されているとおりです。