如何在 R 中编写 fftshift 和 ifftshift？ [关闭]

Question

在numpy中，我们有以下功能：

import numpy
from numpy.fft import fft2, ifft2, fftshift, ifftshift

我想在 R 中重写这些函数。R 中的 fft 与 python 中的 fft 或 fft2 一样。同样对于ifft2，我们必须做fft(,inverse=T)

现在我想知道如何在 R 中高效地重写 fftshift 和 ifftshift 函数（用于矩阵）。

Answer 1

fftshift 和 ifftshift 背后的概念非常简单。这是我从 MathWorks（MATLAB 的创建者）中提取的一个图：

^{Source: MathWorks doc page on fftshift}

假设您的输入 2D 矩阵被分成多个象限。第 1 象限在左上角，第 2 象限在右上角，第 3 象限在右下角，第 4 象限在左下角。对于二维矩阵，fftshift 默认交换第一和第三象限以及第二和第四象限。您可以覆盖此行为，您可以简单地在一个维度上单独执行fftshift。如果你这样做，你正在交换所谓的半空格。如果您指定沿行交换（即维度 1），则矩阵的上半部分与下半部分交换。如果您指定沿列交换（即维度 2），则右半部分与左半部分交换：

^{Source: MathWorks doc page on fftshift}

默认情况下，使用fftshift 将依次交换维度 1 和维度 2。如果您有一个偶数大小的矩阵，其中行和列是偶数，那么将矩阵切割成四个部分并进行交换是非常明确的。但是，如果矩阵是 odd 大小的，则取决于您正在查看的语言。例如，在 MATLAB 和 Python numpy 中，执行切换的位置定义为(r,c) = ceil(rows/2), ceil(cols/2)，其中rows 和cols 是矩阵的行和列。 r 和 c 是发生交换的行和列。

对于ifftshift，您只需反转在fftshift 上执行的操作。因此，默认操作是交换维度 2 和维度 1。但是，对于奇数维度的矩阵，您必须重新定义切换中心的位置。您必须使用floor 而不是ceil，因为这精确地确定了在执行fftshift 之后 半空间的位置，并且您现在正在撤消对原始矩阵所做的操作。因此，新的切换中心是(r,c) = floor(rows/2), floor(cols/2)。除此之外，单个维度之间的交换逻辑是一样的——只是切换的中心已经改变了。

因此，这就是我想出的。这些函数采用 R 中定义的矩阵以及第二个可选参数来确定要交换的维度。省略它会执行我刚才谈到的默认象限交换：

fftshift <- function(input_matrix, dim = -1) {

    rows <- dim(input_matrix)[1]    
    cols <- dim(input_matrix)[2]    

    swap_up_down <- function(input_matrix) {
        rows_half <- ceiling(rows/2)
        return(rbind(input_matrix[((rows_half+1):rows), (1:cols)], input_matrix[(1:rows_half), (1:cols)]))
    }

    swap_left_right <- function(input_matrix) {
        cols_half <- ceiling(cols/2)
        return(cbind(input_matrix[1:rows, ((cols_half+1):cols)], input_matrix[1:rows, 1:cols_half]))
    }

    if (dim == -1) {
        input_matrix <- swap_up_down(input_matrix)
        return(swap_left_right(input_matrix))
    }
    else if (dim == 1) {
        return(swap_up_down(input_matrix))
    }
    else if (dim == 2) {
        return(swap_left_right(input_matrix))
    }
    else {
        stop("Invalid dimension parameter")
    }
}

ifftshift <- function(input_matrix, dim = -1) {

    rows <- dim(input_matrix)[1]    
    cols <- dim(input_matrix)[2]    

    swap_up_down <- function(input_matrix) {
        rows_half <- floor(rows/2)
        return(rbind(input_matrix[((rows_half+1):rows), (1:cols)], input_matrix[(1:rows_half), (1:cols)]))
    }

    swap_left_right <- function(input_matrix) {
        cols_half <- floor(cols/2)
        return(cbind(input_matrix[1:rows, ((cols_half+1):cols)], input_matrix[1:rows, 1:cols_half]))
    }

    if (dim == -1) {
        input_matrix <- swap_left_right(input_matrix)
        return(swap_up_down(input_matrix))
    }
    else if (dim == 1) {
        return(swap_up_down(input_matrix))
    }
    else if (dim == 2) {
        return(swap_left_right(input_matrix))
    }
    else {
        stop("Invalid dimension parameter")
    }
}

在每个函数中，我定义了交换左半部分和右半部分以及上半部分和下半部分的函数。要交换左右两半，只需确定您用于执行交换的列并使用索引通过将这两半连接在一起来创建一个新矩阵，其中前半部分是右半部分，后半部分是左半部分.在交换上半部和下半部时您会做同样的事情，但要找到正确的行来执行交换。

第二个参数dim可以选择设置为1或2来交换对应的维度。请注意，fftshift 和 ifftshift 之间的唯一区别是定义的交换中心以及默认交换两个维度时的操作顺序。其余代码相同。

作为演示的一种方式，假设我已经声明了一个 5 x 5 的数字矩阵，如下所示：

> O <- matrix(1:25, 5, 5)
> O
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    6   11   16   21
[2,]    2    7   12   17   22
[3,]    3    8   13   18   23
[4,]    4    9   14   19   24
[5,]    5   10   15   20   25

请注意，矩阵的大小在两个维度上都是奇数。执行fftshift 会给我们：

> P <- fftshift(O)
> P
     [,1] [,2] [,3] [,4] [,5]
[1,]   19   24    4    9   14
[2,]   20   25    5   10   15
[3,]   16   21    1    6   11
[4,]   17   22    2    7   12
[5,]   18   23    3    8   13

可以看到对应的维度已经交换了。用ifftshift 反转这个应该给我们原来的矩阵，它确实：

> Q <- ifftshift(P)
> Q
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    6   11   16   21
[2,]    2    7   12   17   22
[3,]    3    8   13   18   23
[4,]    4    9   14   19   24
[5,]    5   10   15   20   25

当然，您可以覆盖它并指定要交换的维度，假设您只想交换行：

> fftshift(O, 1)
     [,1] [,2] [,3] [,4] [,5]
[1,]    4    9   14   19   24
[2,]    5   10   15   20   25
[3,]    1    6   11   16   21
[4,]    2    7   12   17   22
[5,]    3    8   13   18   23

> ifftshift(fftshift(O, 1), 1)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    6   11   16   21
[2,]    2    7   12   17   22
[3,]    3    8   13   18   23
[4,]    4    9   14   19   24
[5,]    5   10   15   20   25

请注意，当您对一维交换的矩阵执行 ifftshift 时，您必须使用与使用 fftshift 时用于交换的相同维度，以确保获得原矩阵回来了。

我们也可以对第二个维度做同样的事情：

> ifftshift(O, 2)
     [,1] [,2] [,3] [,4] [,5]
[1,]   11   16   21    1    6
[2,]   12   17   22    2    7
[3,]   13   18   23    3    8
[4,]   14   19   24    4    9
[5,]   15   20   25    5   10

> ifftshift(fftshift(O, 2), 2)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    6   11   16   21
[2,]    2    7   12   17   22
[3,]    3    8   13   18   23
[4,]    4    9   14   19   24
[5,]    5   10   15   20   25

---

作为一个有趣的说明，我们可以验证numpy 与我们上面讨论的相同。这是一个 IPython 会话：

In [16]: import numpy as np

In [17]: from numpy.fft import fftshift, ifftshift

In [18]: O = np.reshape(np.arange(1,26), (5,5)).T

In [19]: O
Out[19]:
array([[ 1,  6, 11, 16, 21],
       [ 2,  7, 12, 17, 22],
       [ 3,  8, 13, 18, 23],
       [ 4,  9, 14, 19, 24],
       [ 5, 10, 15, 20, 25]])

In [20]: fftshift(O)
Out[20]:
array([[19, 24,  4,  9, 14],
       [20, 25,  5, 10, 15],
       [16, 21,  1,  6, 11],
       [17, 22,  2,  7, 12],
       [18, 23,  3,  8, 13]])

In [21]: ifftshift(fftshift(O))
Out[21]:
array([[ 1,  6, 11, 16, 21],
       [ 2,  7, 12, 17, 22],
       [ 3,  8, 13, 18, 23],
       [ 4,  9, 14, 19, 24],
       [ 5, 10, 15, 20, 25]])

numpy 以及numpy.fft 包中的fftshift 和ifftshift 被导入并创建一个与您在R 中定义的示例中看到的相同的二维矩阵。然后我们调用@ 987654365@，然后是fftshift 和ifftshift，我们可以看到我们得到的结果与 R 代码中看到的相同。