fftshift/ifftshift C/C++源代码

Question

Closed. This question is seeking recommendations for software libraries, tutorials, tools, books, or other off-site resources. It does not meet Stack Overflow guidelines. It is not currently accepting answers.

---

我们不允许提出寻求软件库、教程、工具、书籍或其他场外资源推荐的问题。您可以编辑问题，以便用事实和引文来回答。

7 年前已关闭。

Answer 1

FFTHIFT / IFFTSHIFT 是一种执行 CIRCSHIFT 的奇特方法。
您可以验证 FFTSHIFT 是否可以重写为 CIRCSHIFT，如下所示。
您可以在 C/C++ 中定义宏，将 FFTSHIFT 转为 CIRCSHIFT。

A = rand(m, n);
mm = floor(m / 2);
nn = floor(n / 2);
% All three of the following should provide zeros.
circshift(A,[mm, nn]) - fftshift(A)
circshift(A,[mm,  0]) - fftshift(A, 1)
circshift(A,[ 0, nn]) - fftshift(A, 2) 

IFFTSHIFT 可以找到类似的等效项。

循环移位可以通过以下代码非常简单地实现（当然可以使用并行版本进行改进）。

template<class ty>
void circshift(ty *out, const ty *in, int xdim, int ydim, int xshift, int yshift)
{
  for (int i = 0; i < xdim; i++) {
    int ii = (i + xshift) % xdim;
    for (int j = 0; j < ydim; j++) {
      int jj = (j + yshift) % ydim;
      out[ii * ydim + jj] = in[i * ydim + j];
    }
  }
}

然后

#define fftshift(out, in, x, y) circshift(out, in, x, y, (x/2), (y/2))
#define ifftshift(out, in, x, y) circshift(out, in, x, y, ((x+1)/2), ((y+1)/2))

这是临时完成的。如果有任何格式/语法问题，请耐心等待。

Answer 2

这段代码可能会有所帮助。它仅对一个缓冲区内的一维数组执行 fftshift/ifftshift。偶数个元素的前向和后向 fftshift 算法完全相同。

void swap(complex *v1, complex *v2)
{
    complex tmp = *v1;
    *v1 = *v2;
    *v2 = tmp;
}

void fftshift(complex *data, int count)
{
    int k = 0;
    int c = (int) floor((float)count/2);
    // For odd and for even numbers of element use different algorithm
    if (count % 2 == 0)
    {
        for (k = 0; k < c; k++)
            swap(&data[k], &data[k+c]);
    }
    else
    {
        complex tmp = data[0];
        for (k = 0; k < c; k++)
        {
            data[k] = data[c + k + 1];
            data[c + k + 1] = data[k + 1];
        }
        data[c] = tmp;
    }
}

void ifftshift(complex *data, int count)
{
    int k = 0;
    int c = (int) floor((float)count/2);
    if (count % 2 == 0)
    {
        for (k = 0; k < c; k++)
            swap(&data[k], &data[k+c]);
    }
    else
    {
        complex tmp = data[count - 1];
        for (k = c-1; k >= 0; k--)
        {
            data[c + k + 1] = data[k];
            data[k] = data[c + k];
        }
        data[c] = tmp;
    }
}

更新：
另外，任意点数的 FFT 库（包括 fftshift 操作）可以在 Optolithium 中找到（在 OptolithiumC/库/傅立叶）

Answer 3

通常，FFT 的居中通过 v(k)=v(k)*(-1)**k 完成
时域。频域移动是一个很差的替代品，对于
数学原因和计算效率。
参见第 27 页：
http://show.docjava.com/pub/document/jot/v8n6。 pdf

我不确定为什么 Matlab 文档会这样做，
他们没有提供技术参考。

Answer 4

或者您可以通过键入 type fftshift 并用 C++ 重新编码来自己完成此操作。 Matlab 代码并没有那么复杂。

编辑：我注意到这个答案最近被否决了几次，并以负面的方式发表了评论。我记得有一段时间 type fftshift 比当前的实现更具启发性，但我可能是错的。如果我可以删除答案，我会的，因为它似乎不再相关。

这里是一个版本（由 Octave 提供），无需
圆周移位。

Answer 5

我测试了此处提供的代码并制作了一个示例项目来测试它们。对于一维代码，可以简单地使用 std::rotate

template <typename _Real>
static inline
void rotshift(complex<_Real> * complexVector, const size_t count)
{
    int center = (int) floor((float)count/2);
    if (count % 2 != 0) {
        center++;
    }
    // odd: 012 34 changes to 34 012
    std::rotate(complexVector,complexVector + center,complexVector + count);
}

template <typename _Real>
static inline
void irotshift(complex<_Real> * complexVector, const size_t count)
{
    int center = (int) floor((float)count/2);
    // odd: 01 234 changes to 234 01
    std::rotate(complexVector,complexVector +center,complexVector + count);
}

我更喜欢使用 std::rotate 而不是 Alexei 的代码，因为它很简单。

对于 2D 来说，情况变得更加复杂。对于偶数，它基本上是左右翻转和上下翻转。奇怪的是循环移位算法：

//    A =
//    1     2     3
//    4     5     6
//    7     8     9


//    fftshift2D(A)
//    9  |   7     8
//    --------------
//    3  |   1     2
//    6  |   4     5


//    ifftshift2D(A)
//    5     6   |  4
//    8     9   |  7
//    --------------
//    2     3   |  1

这里我通过一个接口实现了循环移位代码，仅使用一个数组进行输入和输出。对于偶数，仅需要一个数组，对于奇数，临时创建第二个数组并将其复制回输入数组。由于复制数组需要额外的时间，这会导致性能下降。

template<class _Real>
static inline
void fftshift2D(complex<_Real> *data, size_t xdim, size_t ydim)
{
    size_t xshift = xdim / 2;
    size_t yshift = ydim / 2;
    if ((xdim*ydim) % 2 != 0) {
        // temp output array
        std::vector<complex<_Real> > out;
        out.resize(xdim * ydim);
        for (size_t x = 0; x < xdim; x++) {
            size_t outX = (x + xshift) % xdim;
            for (size_t y = 0; y < ydim; y++) {
                size_t outY = (y + yshift) % ydim;
                // row-major order
                out[outX + xdim * outY] = data[x + xdim * y];
            }
        }
        // copy out back to data
        copy(out.begin(), out.end(), &data[0]);
        }
    else {
        // in and output array are the same,
        // values are exchanged using swap
        for (size_t x = 0; x < xdim; x++) {
            size_t outX = (x + xshift) % xdim;
            for (size_t y = 0; y < yshift; y++) {
                size_t outY = (y + yshift) % ydim;
                // row-major order
                swap(data[outX + xdim * outY], data[x + xdim * y]);
            }
        }
    }
}

template<class _Real>
static inline
void ifftshift2D(complex<_Real> *data, size_t xdim, size_t ydim)
{
    size_t xshift = xdim / 2;
    if (xdim % 2 != 0) {
        xshift++;
    }
    size_t yshift = ydim / 2;
    if (ydim % 2 != 0) {
        yshift++;
    }
    if ((xdim*ydim) % 2 != 0) {
        // temp output array
        std::vector<complex<_Real> > out;
        out.resize(xdim * ydim);
        for (size_t x = 0; x < xdim; x++) {
            size_t outX = (x + xshift) % xdim;
            for (size_t y = 0; y < ydim; y++) {
                size_t outY = (y + yshift) % ydim;
                // row-major order
                out[outX + xdim * outY] = data[x + xdim * y];
            }
        }
        // copy out back to data
        copy(out.begin(), out.end(), &data[0]);
        }
    else {
        // in and output array are the same,
        // values are exchanged using swap
        for (size_t x = 0; x < xdim; x++) {
            size_t outX = (x + xshift) % xdim;
            for (size_t y = 0; y < yshift; y++) {
                size_t outY = (y + yshift) % ydim;
                // row-major order
                swap(data[outX + xdim * outY], data[x + xdim * y]);
            }
        }
    }
}

Answer 6

注意：提供了更好的答案，我只是将其保留在这里一段时间......我不知道什么。

尝试一下：

template<class T> void ifftShift(T *out, const T* in, size_t nx, size_t ny)
{
    const size_t hlen1 = (ny+1)/2;
    const size_t hlen2 = ny/2;
    const size_t shft1 = ((nx+1)/2)*ny + hlen1;
    const size_t shft2 = (nx/2)*ny + hlen2;

    const T* src = in;
    for(T* tgt = out; tgt < out + shft1 - hlen1; tgt += ny, src += ny) { // (nx+1)/2 times
        copy(src, src+hlen1, tgt + shft2);          //1->4
        copy(src+hlen1, src+ny, tgt+shft2-hlen2); } //2->3
    src = in;
    for(T* tgt = out; tgt < out + shft2 - hlen2; tgt += ny, src += ny ){ // nx/2 times
        copy(src+shft1, src+shft1+hlen2, tgt);         //4->1
        copy(src+shft1-hlen1, src+shft1, tgt+hlen2); } //3->2
};

对于偶数维度的矩阵，您可以就地执行此操作，只需将相同的指针传递给输入和输出参数即可。

另请注意，对于一维数组，fftshift 只是 std::rotate。

Answer 7

您还可以使用 arrayfire 的 shift 函数替代 Matlab 的 circshift 并重新实现其余代码。如果您对 AF 的任何其他功能感兴趣（例如通过简单更改链接器标志即可移植到 GPU），这可能会很有用。

但是，如果您的所有代码都旨在在 CPU 上运行并且非常复杂，或者您不想使用任何其他数据格式（AF 需要 af::arrays），请坚持使用其他选项之一。

我最终更改为 AF，因为否则我必须将 fftshift 重新实现为 OpenCL 内核。

Answer 8

它会给出与 matlab 中的 ifftshift 等效的结果

ifftshift(vector< vector <double> > Hlow,int RowLineSpace, int ColumnLineSpace)
{
   int pivotRow=floor(RowLineSpace/2);
   int pivotCol=floor(ColumnLineSpace/2);

   for(int i=pivotRow;i<RowLineSpace;i++){
      for(int j=0;j<ColumnLineSpace;j++){
        double temp=Hlow.at(i).at(j);
        second.push_back(temp);
      }
    ifftShiftRow.push_back(second);
    second.clear();
}

for(int i=0;i<pivotRow;i++){
        for(int j=0;j<ColumnLineSpace;j++){
            double temp=Hlow.at(i).at(j);
            first.push_back(temp);
        }
        ifftShiftRow.push_back(first);
        first.clear();
    }


 double** arr = new double*[RowLineSpace];
  for(int i = 0; i < RowLineSpace; ++i)
      arr[i] = new double[ColumnLineSpace];
  int i1=0,j1=0;
  for(int j=pivotCol;j<ColumnLineSpace;j++){
        for(int i=0;i<RowLineSpace;i++){
            double temp2=ifftShiftRow.at(i).at(j);
            arr[i1][j1]=temp2;
            i1++;
        }
        j1++;
        i1=0;
    }
 for(int j=0;j<pivotCol;j++){
    for(int i=0;i<RowLineSpace;i++){
        double temp1=ifftShiftRow.at(i).at(j);
            arr[i1][j1]=temp1;
                i1++;

            }

           j1++;
            i1=0;
        }

for(int i=0;i<RowLineSpace;i++){
    for(int j=0;j<ColumnLineSpace;j++){
        double value=arr[i][j];
        temp.push_back(value);

    }
    ifftShiftLow.push_back(temp);
    temp.clear();
}
    return ifftShiftLow;


 }

Answer 9

八度 使用 fftw 实施(i)快速移位。

Answer 10

您可以使用kissfft。它速度相当快，使用起来非常简单，而且免费。按照您想要的方式排列输出只需要：

a) 移位 (-dim_x/2, -dim_y/2, ...)，具有周期性边界条件

b) FFT 或 IFFT

c) 向后移位 (dim_x/2, dim_y/2, ...) ，具有周期性边界条件

d) 尺度 ? （根据您的需要 IFFT*FFT 默认情况下会按 dim_x*dim_y*... 缩放函数）