fftshift/ifftshift C/C++ 源代码 [关闭]
Posted
技术标签:
【中文标题】fftshift/ifftshift C/C++ 源代码 [关闭]【英文标题】：fftshift/ifftshift C/C++ source code [closed] 【发布时间】：2011-05-06 17:50:34 【问题描述】：有谁知道是否有任何免费的开源库按照 matlab 中定义的方式实现了这两个功能？
谢谢
【问题讨论】：
我投错票了。尽管this question 相似，但您的问题并不完全重复，因为您要求免费/开放库，而其他人只是想将 MATLAB 的实现移植到 C。 一些方便的shift C++函数可以在FFTW++中找到：github.com/dealias/fftwpp/blob/master/fftw%2B%2B.h 除上述之外，NFFT库在C语言中提供了移位函数：github.com/NFFT/nfft/blob/… 【参考方案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))
这有点即兴。如果有任何格式/语法问题，请多多包涵。
【讨论】：
【参考方案2】：通常，FFT 的中心化是通过 v(k)=v(k)*(-1)**k in 时域。在频域中移动是一个糟糕的替代品，因为 数学原因和计算效率。 见第 27 页： http://show.docjava.com/pub/document/jot/v8n6.pdf
我不确定为什么 Matlab 文档会按照他们的方式进行操作， 他们没有提供技术参考。
【讨论】：
+1。就代码而言，对于行向量x，fftshift(fft(x)) 在fft(x .* (-1).^(0 : numel(x) - 1)) 的机器精度范围内。 ifftshift（以及二维情况）也有类似的结果。如果正在计算 FFT 的输入，有时很容易在 FFT 之前插入符号交替，这比在 FFT 之后移动内存要快得多。
这是否也适用于元素数量为奇数的情况？
不，它不适用于奇数个元素。对于奇数个元素，您需要将输入乘以exp(i*pi*numel(x)/(numel(x)+1)*(0 : numel(x) - 1))。但是当我这样做时，我会遇到数值精度问题。请注意，即使对于实际输入，结果也很复杂，因此您将无法就地执行此操作。对于奇数大小的信号，最好在傅里叶域中移动数据。我不明白为什么它是一个糟糕的替代品，当然我看不到数学上的原因！【参考方案3】：
此代码可能会有所帮助。它仅对一个缓冲区内的一维数组执行 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/libs/fourier 下）
【讨论】：
【参考方案4】：或者您可以自己输入type fftshift 并用C++ 重新编码。 Matlab 代码没那么复杂。
编辑：我注意到这个答案最近几次被否决，并以负面的方式发表评论。我记得有一次type fftshift 比当前的实现更具启发性，但我可能是错的。如果我可以删除答案，我会因为它似乎不再相关。
Here 是一个没有实现它的版本（由 Octave 提供）
circshift.
【讨论】：
fftshift 被实现为对circshift 的调用，而circshift 是一个内置函数，所以你看不到它是如何编码的。这是一个无用的答案。【参考方案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，它变得更加复杂。对于偶数，它基本上是左右翻转和上下翻转。奇怪的是 circshift 算法：
// 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
在这里，我使用一个接口实现了 circshift 代码，该接口仅使用一个数组用于输入和输出。对于偶数，只需要一个数组，对于奇数，临时创建第二个数组并将其复制回输入数组。由于复制数组需要额外的时间，这会导致性能下降。
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]);
【讨论】：
【参考方案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
;
对于具有偶数维度的矩阵，您可以就地执行此操作，只需将相同的指针传递给 in 和 out 参数。
另请注意，对于一维数组，fftshift 只是 std::rotate。
【讨论】：
您能否提供此代码的示例。 IT 对一维或二维数组的计算和处理内容并不十分明显。 这仅适用于二维数组。它所做的是在知道正确的四分之一的源数组和目标数组之间进行半行复制。我自己实际上不再喜欢这段代码。我建议从 Pavan Yalamanchili 回答开始，只是不要定义函数可以执行的宏。 该示例的问题在于，它需要一个输入和输出数组。然而，对于大型矩阵，这需要分配我想避免的额外内存。我喜欢的是它是更通用的解决方案。 对于偶数维度，您可以将相同的矩阵用于输入和输出。但实际上，这里有更好的答案。例如，喜欢在真实空间领域做（虽然我不喜欢那里的推理，但这个想法本身是值得的）。【参考方案7】：您还可以使用 arrayfire 的 shift 函数代替 Matlab 的 circshift 并重新实现其余代码。如果您对 AF 的任何其他功能（例如通过简单地更改链接器标志可移植到 GPU）感兴趣，这可能会很有用。
但是，如果您的所有代码都打算在 CPU 上运行并且非常复杂，或者您不想使用任何其他数据格式（AF 需要 af::arrays），请坚持使用其他选项之一。
我最终改用 AF，因为我不得不将 fftshift 重新实现为 OpenCL 内核，否则回到过去。
【讨论】：
【参考方案8】：Octaveusesfftw 实现 (i)fftshift。
【讨论】：
fftw 没有任何我认为没有的转换功能的 API 无需进入频域再返回即可执行这些简单的移位操作。【参考方案9】：它会给出与 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;
【讨论】：
【参考方案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*... 缩放函数）
【讨论】：
以上是关于fftshift/ifftshift C/C++ 源代码 [关闭]的主要内容，如果未能解决你的问题，请参考以下文章