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[0];
    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

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

> 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
>     quadrants and the second and fourth quadrants.  For N-D
>     arrays, FFTSHIFT(X) swaps "half-spaces" of X along each
>     dimension.
> --
>
> Thanks in advance,
> 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

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
    quadrants and the second and fourth quadrants.  For N-D
    arrays, FFTSHIFT(X) swaps "half-spaces" of X along each
    dimension.
--

Thanks in advance,
Graham