fftshift/ifftshift C/C++ source code

Question

Does anyone know if there is any free and open source library that has implemented these two functions the way they are defined in matlab?

Thanks

Answer 1

FFTHIFT / IFFTSHIFT is a fancy way of doing CIRCSHIFT. You can verify that FFTSHIFT can be rewritten as CIRCSHIFT as following. You can define macros in C/C++ to punt FFTSHIFT to 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) 

Similar equivalents can be found for IFFTSHIFT.

Circular shift can be implemented very simply with the following code (Can be improved with parallel versions ofcourse).

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];
    }
  }
}

And then

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

This was done a bit impromptu. Bear with me if there are any formatting / syntactical problems.

Answer 2

Possible this code may help. It perform fftshift/ifftshift only for 1D array within one buffer. Algorithm of forward and backward fftshift for even number of elements is fully identical.

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;
    }
}

UPDATED: Also FFT library (including fftshift operations) for arbitrary points number could be found in Optolithium (under the OptolithiumC/libs/fourier)

Answer 3

Normally, centering the FFT is done with v(k)=v(k)*(-1)**k in the time domain. Shifting in the frequency domain is a poor substitute, for mathematical reasons and for computational efficiency. See pp 27 of: http://show.docjava.com/pub/document/jot/v8n6.pdf

I am not sure why Matlab documentation does it the way they do, they give no technical reference.

Answer 4

Or you can do it yourself by typing type fftshift and recoding that in C++. It's not that complicated of Matlab code.

Edit: I've noticed that this answer has been down-voted a few times recently and commented on in a negative way. I recall a time when type fftshift was more revealing than the current implementation, but I could be wrong. If I could delete the answer, I would as it seems no longer relevant.

Here is a version (courtesy of Octave) that implements it without circshift.

Answer 5

I tested the code provided here and made an example project to test them. For 1D code one can simply use 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);
}

I prefer using std::rotate over the code from Alexei due to its simplicity.

For 2D it gets more complicated. For even numbers it is basically a flip left right and flip upside down. For odd it is the circshift algorithm:

//    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

Here I implemented the circshift code with an interface using only one array for in and output. For even numbers only a single array is required, for odd numbers a second array is temporarily created and copied back to the input array. This causes a performance decrease because of the additional time for copying the array.

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

You could also use arrayfire's shift function as replacement for Matlab's circshift and re-implement the rest of the code. This could be useful if you are interested in any of the other features of AF anyway (such as portability to GPU by simply changing a linker flag).

However if all your code is meant to be run on the CPU and is quite sophisticated or you don't want to use any other data format (AF requires af::arrays) stick with one of the other options.

I ended up changing to AF because I would have had to re-implement fftshift as an OpenCL kernel otherwise back in the time.

Answer 7

Notice: There are better answers provided, I just keep this here for a while for... I do not know what.

Try this:

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
};

For matrices with even dimensions you can do it in-place, just passing the same pointer into in and out parameters.

Also note that for 1D arrays fftshift is just std::rotate.

Answer 8

It will give equivalent result to ifftshift in matlab

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

You can use kissfft. It's reasonable fast, extremely simple to use, and free. Arranging the output like you want it requires only to:

a) shift by (-dim_x/2, -dim_y/2, ...), with periodic boundary conditions

b) FFT or IFFT

c) shift back by (dim_x/2, dim_y/2, ...) , with periodic boundary conditions

d) scale ? (according to your needs IFFT*FFT will scale the function by dim_x*dim_y*... by default)

Answer 10

Octave uses fftw to implement (i)fftshift.