All geek questions in one place

How to interpret the output of numpy.fft.rfft2?

Obviously, the rfft2 function simply computes the discrete fft of the input matrix. However, how do I interpret this output index? Given the exit index, what Fourier coefficient am I looking at?
I am especially confused about the size of the release. For an nn-matrix, the result is represented by the matrix n by (n / 2) +1 (for even n). Why does the square matrix end with a non-square Fourier transform?

+4
source share
3 answers

The conclusion numpy.fft.rfft2is simply the left half (plus one column) of the standard two-dimensional FFT, calculated by numpy.fft.fft2. It is not necessary rfft2to provide the correct half of the result, because the FFT of the real array has natural and simple symmetry , and therefore the right half of the full FFT can be obtained from the left half using this symmetry.

Here is an example to illustrate. First, to be easy to reproduce and easy to see, I will set the random state of NumPy and print options:

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)

Create a real input array with 6 rows and 6 columns:

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]])

Now take a look at the full FFT (using fft2, not rfft2):

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

, , axes np.fft.rfft2:

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]])

documentation np.fft.rfftn, NumPy FFT , FFT .

, rfft2_result : , , . [0, 0], [0, 3], [3, 0] [3, 3] , . .

+2

:

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

rfftn . . 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

, rfft, fftn. , rfft , fftn .

, . fft.

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

2d fftn, 1d fft, 2d-:

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)
+2

fft:

: 1- 0 ( ), n/2 + 1 . length-10:

np.fft.fftfreq(10) : array([ 0. , 0.1, 0.2, 0.3, 0.4, -0.5, -0.4, -0.3, -0.2, -0.1])

np.fft.fftshift(cf), cf=np.fft.fft(array), , : array([-0.5, -0.4, -0.3, -0.2, -0.1, 0. , 0.1, 0.2, 0.3, 0.4]) .

. fft2 rfft2 .

+1

Source: https://habr.com/ru/post/1673126/

More articles:


All Articles