# fftshift equivilent in C

Started by August 23, 2003
```Does anyone have an equivilent to the matlab/octave/sci.py 'fftshift'
routine that is written in C?  I am using FFTW directly to compute
rank 1 and 2 FFTs but need to 'shift' the zero frequency component.

--
FFTSHIFT
Shift zero-frequency component to center of spectrum.
For vectors, FFTSHIFT(X) swaps the left and right halves of
X.  For matrices, FFTSHIFT(X) swaps the first and third
arrays, FFTSHIFT(X) swaps "half-spaces" of X along each
dimension.
--

Graham
```
```> Does anyone have an equivilent to the matlab/octave/sci.py 'fftshift'
> routine that is written in C?  I am using FFTW directly to compute
> rank 1 and 2 FFTs but need to 'shift' the zero frequency component.
>
> --
> FFTSHIFT
>     Shift zero-frequency component to center of spectrum.
>     For vectors, FFTSHIFT(X) swaps the left and right halves of
>     X.  For matrices, FFTSHIFT(X) swaps the first and third
>     arrays, FFTSHIFT(X) swaps "half-spaces" of X along each
>     dimension.
> --
>
> Graham

Graham,

to reorder in 1-D:

int n;         // FFT output vector length
int n2;
int i;
complex x[n];  // you might typedef your own complex type
complex tmp;

n2 = n / 2;    // half of vector length

for (i = 0; i < n2; i++)
{
tmp     = x[i];
x[i]    = x[i+n2];
x[i+n2] = tmp;
}

to reorder in 2-D:

int m, n;      // FFT row and column dimensions might be different
int m2, n2;
int i, k;
complex x[m][n];
complex tmp13, tmp24;

m2 = m / 2;    // half of row dimension
n2 = n / 2;    // half of column dimension

// interchange entries in 4 quadrants, 1 <--> 3 and 2 <--> 4

for (i = 0; i < m2; i++)
{
for (k = 0; k < n2; k++)
{
tmp13         = x[i][k];
x[i][k]       = x[i+m2][k+n2];
x[i+m2][k+n2] = tmp13;

tmp24         = x[i+m2][k];
x[i+m2][k]    = x[i][k+n2];
x[i][k+n2]    = tmp24;
}
}

Hope this helps,

Andrew
```
```On 23 Aug 2003, Andrew Nesterov wrote:

> int n;         // FFT output vector length
> int n2;
> int i;
> complex x[n];  // you might typedef your own complex type
> complex tmp;
>
> n2 = n / 2;    // half of vector length
>
> for (i = 0; i < n2; i++)
> {
>      tmp     = x[i];
>      x[i]    = x[i+n2];
>      x[i+n2] = tmp;
> }

This is close to what I came up with but it only works for even n.

Matlab gives:

>> fftshift([1,2,3,4,5,6,7,8])
ans =     5     6     7     8     1     2     3     4

>> fftshift([1,2,3,4,5,6,7,8,9])
ans =     6     7     8     9     1     2     3     4     5

Above gives:

n = 8
n2 = 4
in = 1.00   2.00   3.00   4.00   5.00   6.00   7.00   8.00
out = 5.00   6.00   7.00   8.00   1.00   2.00   3.00   4.00

n = 9
n2 = 4
in = 1.00   2.00   3.00   4.00   5.00   6.00   7.00   8.00   9.00
out = 5.00   6.00   7.00   8.00   1.00   2.00   3.00   4.00   9.00

Any ideas on how to get the same answer as matlab for both even and odd n?

Thanks,
Graham
```
```Graham,

It can be done out of place rather simply, but note a new array, z[n]
was introduced. It also works for even n:

for (i = 0; i < n; i++)
{
//   find circulant index
k    = (i+n2) % n;
//   or, without modulo call
k    = (i+n2);
if (k >= n) k = k - n;

//   move
z[k] = x[i];
}

For this particular case, an odd n and shift by [n/2], it might be
done in place, but more complicated. As it is seen, the only problem
is the overlap of i-th and (i+n2)-th entries:

0  1  2  3  4  5  6  7  8
5  6  7  8  0  1  2  3  4
^
0-th entry is moved into 4-th location, but 4-th location instead of
been moved to 0-th position, as in an even length case, is moved into
the last location, etc. This happens due to curcular addressing mode:
e.g. n  = 9, n2 = 4 and 0+4 mod 9 = 4, 4+4 mod 9 = 8, but 5+4 mod 9 = 0
and the 5th entry eventually goes into the first location.

In an even length case, e.g. n = 8, the 4-th location is moved into
4+4 mod 8 = 0th location, thus the circulant shift might be accomplished
by swap operations.

The code works for even n as well (again, n2 = n/2):

for (i = 0; i < n2; i++) // swap what can be swapped
{
tmp       = x[i];
x[i]      = x[i+n2];
x[i+n2]   = tmp;
}
if (n & 1)               // odd n, shift the rest
{
tmp    = x[n-1];
x[n-1] = x;
for (i = 1; i < n2; i++)
{
x[i-1] = x[i];
}
x[n2-1] = tmp;
}

It might be done in place more simply, with 4 load/store ops less,
but works only for odd n. For even n, (ip1+n2) is n when i becomes
n2-1, thus load/store is being done out of array bound on the last
iteration. For odd n, (ip1+n2) is n-1 on the last iteration:

for (i = 0; i < n2; i++)
{
ip1       = i + 1
tmp       = x[i+n2];
x[i+n2]   = x[i];
x[i]      = x[ip1+n2];
x[ip1+n2] = tmp;
}

Regards,

Andrew

> > int n;         // FFT output vector length
> > int n2;
> > int i;
> > complex x[n];  // you might typedef your own complex type
> > complex tmp;
> >
> > n2 = n / 2;    // half of vector length
> >
> > for (i = 0; i < n2; i++)
> > {
> >      tmp     = x[i];
> >      x[i]    = x[i+n2];
> >      x[i+n2] = tmp;
> > }
>
> This is close to what I came up with but it only works for even n.
>
> Matlab gives:
>
> >> fftshift([1,2,3,4,5,6,7,8])
> ans =     5     6     7     8     1     2     3     4
>
> >> fftshift([1,2,3,4,5,6,7,8,9])
> ans =     6     7     8     9     1     2     3     4     5
>
> Above gives:
>
> n = 8
> n2 = 4
> in = 1.00   2.00   3.00   4.00   5.00   6.00   7.00   8.00
> out = 5.00   6.00   7.00   8.00   1.00   2.00   3.00   4.00
>
> n = 9
> n2 = 4
> in  = 1.00   2.00   3.00   4.00   5.00   6.00   7.00   8.00   9.00
> out = 5.00   6.00   7.00   8.00   1.00   2.00   3.00   4.00   9.00
>
> Any ideas on how to get the same answer as matlab for both even and odd n?
>
> Thanks,
> Graham
```