我有一个C++代码,它对给定的meshgrid
执行fftshift
。
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;
这给出以下移位矩阵:
{{{-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.}}}
而对于Python,相同的fftshift
()会导致:
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]]]
如何使c++ fftshift的输出矩阵与scipy的输出矩阵完全相等?
我尝试使用xt::roll
、xt::transpose + xt::swap
和手动循环移位组合,但都不起作用。
更新:尝试使用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);
}
然而,由于某种原因,仅获得与scipy.fft.fftshift相同的正确矩阵,其大小= 5或大小= 125。我不知道为什么?
更新2:根据@chris的回答,我添加了滚动手动换档。它似乎复制了scipy的fftshift,但似乎相当慢。
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());
}
}
1条答案
按热度按时间wtlkbnrh1#
xtensor
fftshift
在the source code中有注解:部分模仿np.fftshift(仅1D阵列)
所以,它不会为你的2D案例工作。
roll
应该可以。我找不到任何关于xtensor的有用文档,但是这个函数的声明给出了很好的提示:auto roll(E&& e,std::ptrdiff_t shift,std::ptrdiff_t axis);
因此,您需要依次应用
roll
每个轴。我不知道我是哪个方向的转变发生,你需要实验一下。移位距离应为size / 2
,对于fftshift
,在一个方向,对于ifftshift
,在另一个方向。请注意,
fftshift
将原点从左上角移到中间,ifftshift
将原点从中间移到角落。对于偶数大小的数组,这是完全相同的事情,但对于奇数大小的数组,如OP中的数组,则不是这样。OP有一个数组,原点位于非标准点(不是中间,不是角落),应用
fftshift
恰好将原点移到角落。但这并不是使用此函数的正常方式。要定义原点在中间的坐标系(对于偶数和奇数大小的数组),请执行以下操作:
现在,应用
ifftshift
将原点移动到角点,而在此结果上应用fftshift
将原点移动回原来的位置。