fftshift/ifftshift C/C++ source code [closed]

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.

Normally, centering the FFT is done with v(k)=v(k)*(-1)**k inthe time domain. Shifting in the frequency domain is a poor substitute, formathematical 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.

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)