scipy fftshift and fftw fftshift gives different values (shifted matrix are not the same)

Question

I have a C++ code which performs fftshift for a given meshgrid.

    const std::size_t size = 5;
    auto ar = xt::meshgrid(xt::arange<double>(0, size), xt::arange<double>(0, size));
    int translate = (size + 1) / 2;
    
    xt::xarray<double> x = std::get<0>(ar) - translate;
    xt::xarray<double> y = std::get<1>(ar) - translate;
    xt::xarray<double> xy_ = xt::stack(xt::xtuple(x, y));
    auto p = xt::fftw::fftshift(xy_);
    std::cout << p << std::endl;

Which gives the following shifted matrix:

{{{-3., -2., -1.,  0.,  1.},
  {-3., -2., -1.,  0.,  1.},
  {-3., -2., -1.,  0.,  1.},
  {-3., -2., -1.,  0.,  1.},
  {-3., -2., -1.,  0.,  1.}},
 {{-1., -1., -1., -1., -1.},
  { 0.,  0.,  0.,  0.,  0.},
  { 1.,  1.,  1.,  1.,  1.},
  {-3., -3., -3., -3., -3.},
  {-2., -2., -2., -2., -2.}}}

whereas for python the same fftshift ()results in:

np.mgrid[:size, :size] - int( (size + 1)/2 )
fftshifted_mat = scipy.fftpack.fftshift(mat)
print(fftshifted_mat)

[[[ 0  1 -3 -2 -1]
  [ 0  1 -3 -2 -1]
  [ 0  1 -3 -2 -1]
  [ 0  1 -3 -2 -1]
  [ 0  1 -3 -2 -1]]

 [[ 0  0  0  0  0]
  [ 1  1  1  1  1]
  [-3 -3 -3 -3 -3]
  [-2 -2 -2 -2 -2]
  [-1 -1 -1 -1 -1]]]

How can I make the c++ fftshift output matrix exactly equal to scipy's output matrix?

I tried using xt::roll, xt::transpose + xt::swap, and manual circular shift combinations, but none of them worked.

Update: Tried using roll

    for (std::size_t axis = 0; axis < xy_.shape().size(); ++axis){
        std::size_t dim_size = xy_.shape()[axis];
        std::size_t shift = (dim_size - 1) / 2;
        xy_ = xt::roll(xy_, shift, axis);
    }

However, for some reason only getting right matrix that is same as the scipy.fft.fftshift one with size = 5 or size = 125. I am not sure why?

Update 2: As per @chris' answer I added manual shift with roll. It seems to replicate scipy's fftshift but seems to be quite slow.

template <typename T>
void fftshift_roll(xt::xarray<T>& array)
{
    std::size_t ndims = array.dimension();
    std::vector<std::ptrdiff_t> shift_indices(ndims);

    for (std::size_t i = 0; i < ndims; ++i) {
        std::ptrdiff_t shift = static_cast<std::ptrdiff_t>(array.shape(i)) / 2;
        shift_indices[i] = shift;
    }

    for (std::size_t i = 0; i < ndims; ++i) {
        auto rolled = xt::roll(array, shift_indices[i], i);
        array = xt::view(rolled, xt::all(), xt::all()); 
    }
}

Maybe copy back into `array` only after the loop? Something like `xt::xarray rolled = array; for() { rolled = xt::roll(rolled, ...); } array = ...;`? I haven't used xtensor before, so I don't know if `array = xt::view(rolled, xt::all(), xt::all())` is the best way to copy data over, but it seems weird syntax to me. — Cris Luengo, Jun 23 '23 at 22:25
xtensor has some terrible documentation... not much seems to be explained. — Cris Luengo, Jun 23 '23 at 22:34

Cris Luengo · Answer 1 · 2023-06-23T15:31:25.673

1

xtensor fftshift has a comment in the source code:

partly mimic np.fftshift (only 1D arrays)

So, it won’t work for your 2D case.

roll should do the trick. I can’t find any useful documentation for xtensor, but the function declaration for this one gives good hints:

auto roll(E&& e, std::ptrdiff_t shift, std::ptrdiff_t axis);

So you need to apply roll your each axis in turn. I don’t know I’m which direction the shift happens, you’ll need to experiment a bit. The shift distance should be size / 2, in one direction for fftshift, in the other direction for ifftshift.

Do note that fftshift shifts the origin from the top-left corner to the middle, and ifftshift shifts it from the middle to the corner. For even-sized arrays this is exactly the same thing, but for odd-sized ones like the one in OP it is not.

OP has an array with the origin in a non-standard spot (not the middle, not the corner), and applying fftshift just happens to shift it to the corner. But this is not the normal way of using this function.

To define a coordinate system with the origin in the middle (for both even and odd-sized arrays) do

np.mgrid[:size, :size] - size // 2

Now applying ifftshift will shift the origin to the corner, and applying fftshift on this result will move the origin back to where it was before.

edited Jun 23 '23 at 15:31

answered Jun 23 '23 at 15:00

Cris Luengo

55,762
10
62
120

thanks for the answer. Alright. I already tried roll (please see the edited question above) but that doesn't work. Could you suggest what exactly I could play around with while experimenting with `roll`? I didn't get the second part. You mean that I don't need to use fftshift at all? You last code snippet is already there `np.mgrid[:size, :size] - int( (size + 1)/2 )` – Jarwin Jun 23 '23 at 15:17
@Jarwin I've edited the answer to include the right distance to shift. I didn't understand your sentence "only getting right that corrosponds to scipy one with size = 5 or size = 125." ---- The second part is just to point out that the way you define your coordinate system is not standard, and does not match what `fftshift`/`ifftshift` do. `(size + 1)/2` is not the center of the array according to these functions, it's must be `size/2` (integer division). – Cris Luengo Jun 23 '23 at 15:35
updated the question. I meant the right matrix. And right size/2 I get that now. I'll try the roll directions as with your answer. But it seems slower than scipy already. On a different note, are you aware of any other lib that performs fftshift similar to scipy? – Jarwin Jun 23 '23 at 15:45
@Jarwin I implemented the (`i`)`fftshift` operation in C++ following [this answer](https://stackoverflow.com/a/19752002/7328782), which makes it a very efficient in-place operation. They're the `ShiftCornerToCenter()` and `ShiftCenterToCorner()` functions at the top of [this file](https://github.com/DIPlib/diplib/blob/master/src/transform/fourier.cpp). Of course you can use the DIPlib library directly, it's `dip::FourierTransform()` has an option for whether the origin must be in the middle or in the corner, so there's no need to separately call a shift function. – Cris Luengo Jun 23 '23 at 16:44
ah alright. I'll have to look into it. dip also seems quite straightforward to use. However, since I am using the fftshift for GRF for a FEM library - it would be prudent to use something lightweight / headeronly ":D – Jarwin Jun 23 '23 at 22:25
I updated my question with an attempt. Could you please suggest again how I could implement `ShiftCornerToCenter()` in that? – Jarwin Jun 23 '23 at 22:27
@Jarwin That function runs on a single row or column of data. You'd have to loop over all rows, then over all columns. I don't know enough about the xtensor internals to quickly show you how to do this, sorry. – Cris Luengo Jun 23 '23 at 22:35

scipy fftshift and fftw fftshift gives different values (shifted matrix are not the same)

1 Answers1