Как написать fftshift и ifftshift в R?

Концепция, стоящая за fftshift и ifftshift, довольно проста. Вот рисунок, который я вытащил из MathWorks (создатели MATLAB):

Источник: страница документа MathWorks на fftshift

Представьте, что ваша входная 2D-матрица разбита на квадранты. Квадрант № 1 находится в верхнем левом углу, квадрант № 2 — в верхнем правом, квадрант № 3 — в нижнем правом и квадрант № 4 — в нижнем левом. Для двумерных матриц fftshift по умолчанию меняет местами первый и третий квадранты, а также второй и четвертый квадранты. Вы можете переопределить это поведение, просто выполняя fftshift для одного измерения отдельно. Если вы сделаете это, вы поменяете местами так называемые полупространства. Если вы указали поменять местами строки (т. е. измерение 1), то верхняя половина матрицы заменяется нижней половиной. Если вы указали поменять местами столбцы (т. е. измерение 2), то правая половина заменяется левой половиной:

Источник: страница документа MathWorks на fftshift

Использование fftshift по умолчанию последовательно меняет местами измерения 1 и 2. Если у вас есть матрица четного размера, в которой строки и столбцы четные, то очень однозначно разрезать матрицу на четыре части и поменять местами. Однако если матрица нечетного размера, это зависит от того, какой язык вы просматриваете. Например, в MATLAB и Python numpy место для переключения определяется как (r,c) = ceil(rows/2), ceil(cols/2), где rows и cols — это строки и столбцы матрицы. r и c — это строка и столбец, в которых происходит обмен.

Для ifftshift вы просто отменяете действия, выполненные для fftshift. Таким образом, действие по умолчанию состоит в том, чтобы поменять местами размеры 2, а затем размеры 1. Однако вам необходимо переопределить, где находится центр переключения для матриц нечетных размеров. Вместо ceil вы должны использовать floor, потому что это точно определяет, где были полупространства после выполнения 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, а также fftshift и ifftshift из пакета numpy.fft импортируются, и создается 2D-матрица, аналогичная той, что вы видели в примере, определенном в R. Затем мы последовательно вызываем fftshift, затем fftshift и ifftshift, и мы можем видеть что мы получаем те же результаты, что и в коде R.