0

我正在尝试使用for循环在 c++ 中实现 fftshift 函数（来自 MATLAB），这真的很耗时。这是我的代码：

const int a = 3;
    const int b = 4;
    const int c = 5;
    int i, j, k;
    int aa = a / 2;
    int bb = b / 2;
    int cc = c / 2;

    double ***te, ***tempa;
    te = new double **[a];
    tempa = new double **[a];
    for (i = 0; i < a; i++)
    {
        te[i] = new double *[b];
        tempa[i] = new double *[b];
        for (j = 0; j < b; j++)
        {
            te[i][j] = new double [c];
            tempa[i][j] = new double [c];
            for (k = 0; k < c; k++)
            {
                te[i][j][k] = i + j+k;
            }

        }
    }
    /*for the row*/
    if (c % 2 == 1)
    {
        for (i = 0; i < a; i++)
        {
            for (j = 0; j < b; j++)
            {
                for (k = 0; k < cc; k++)
                {
                    tempa[i][j][k] = te[i][j][k + cc + 1];
                    tempa[i][j][k + cc] = te[i][j][k];
                    tempa[i][j][c - 1] = te[i][j][cc];
                }
            }
        }
    }
    else
    {
        for (i = 0; i < a; i++)
        {
            for (j = 0; j < b; j++)
            {
                for (k = 0; k < cc; k++)
                {
                    tempa[i][j][k] = te[i][j][k + cc];
                    tempa[i][j][k + cc] = te[i][j][k];
                }
            }
        }
    }

    for (i = 0; i < a; i++)
    {
        for (j = 0; j < b; j++)
        {
            for (k = 0; k < c; k++)
            {
                te[i][j][k] = tempa[i][j][k];
            }
        }
    }

    /*for the column*/
    if (b % 2 == 1)
    {
        for (i = 0; i < a; i++)
        {
            for (j = 0; j < bb; j++)
            {
                for (k = 0; k < c; k++)
                {
                    tempa[i][j][k] = te[i][j + bb + 1][k];
                    tempa[i][j + bb][k] = te[i][j][k];
                    tempa[i][b - 1][k] = te[i][bb][k];
                }
            }
        }
    }
    else
    {
        for (i = 0; i < a; i++)
        {
            for (j = 0; j < bb; j++)
            {
                for (k = 0; k < c; k++)
                {
                    tempa[i][j][k] = te[i][j + bb][k];
                    tempa[i][j + bb][k] = te[i][j][k];
                }
            }
        }
    }

    for (i = 0; i < a; i++)
    {
        for (j = 0; j < b; j++)
        {
            for (k = 0; k < c; k++)
            {
                te[i][j][k] = tempa[i][j][k];
            }
        }
    }

    /*for the third dimension*/
    if (a % 2 == 1)
    {

        for ( i = 0; i < aa; i++)
        {
            for (j = 0; j < b; j++)
            {
                for ( k = 0; k < c; k++)
                {
                    tempa[i][j][k] = te[i + aa + 1][j][k];
                    tempa[i + aa][j][k] = te[i][j][k];
                    tempa[a - 1][j][k] = te[aa][j][k];
                }
            }
        }
    }
    else
    {
        for (i = 0; i < aa; i++)
        {
            for ( j = 0; j < b; j++)
            {
                for ( k = 0; k < c; k++)
                {
                    tempa[i][j][k] = te[i + aa][j][k];
                    tempa[i + aa][j][k] = te[i][j][k];

                }
            }
        }
    }

    for (i = 0; i < a; i++)
    {
        for (j = 0; j < b; j++)
        {
            for (k = 0; k < c; k++)
            {
                cout << te[i][j][k] << ' ';
            }
            cout << endl;
        }
        cout << "\n";
    }
    cout << "and then" << endl;
    for (i = 0; i < a; i++)
    {
        for (j = 0; j < b; j++)
        {
            for (k = 0; k < c; k++)
            {
                cout << tempa[i][j][k] << ' ';
            }
            cout << endl;
        }
        cout << "\n";
    }

现在我想重写它memmove以提高运行效率。对于第三维，我使用：

    memmove(tempa, te + aa, sizeof(double)*(a - aa));
    memmove(tempa + aa+1, te, sizeof(double)* aa);

此代码适用于 1d 和 2d 数组，但不适用于 3d 数组。另外，我不知道如何移动列和行元素memmove。任何人都可以帮助我解决所有这些问题吗？非常感谢！！

现在我修改了代码如下：

double ***te, ***tempa1,***tempa2, ***tempa3;

    te = new double **[a];
    tempa1 = new double **[a];
    tempa2 = new double **[a];
    tempa3 = new double **[a];
    for (i = 0; i < a; i++)
    {
        te[i] = new double *[b];
        tempa1[i] = new double *[b];
        tempa2[i] = new double *[b];
        tempa3[i] = new double *[b];
        for (j = 0; j < b; j++)
        {
            te[i][j] = new double [c];
            tempa1[i][j] = new double [c];
            tempa2[i][j] = new double [c];
            tempa3[i][j] = new double [c];
            for (k = 0; k < c; k++)
            {
                te[i][j][k] = i + j+k;
            }

        }
    }

    /*for the third dimension*/
    memmove(tempa1, te + (a-aa), sizeof(double**)*aa);
    memmove(tempa1 + aa, te, sizeof(double**)* (a-aa));
    //memmove(te, tempa, sizeof(double)*a);
    /*for the row*/
    for (i = 0; i < a; i++)
    {
        memmove(tempa2[i], tempa1[i] + (b - bb), sizeof(double*)*bb);
        memmove(tempa2[i] + bb, tempa1[i], sizeof(double*)*(b - bb));
    }

    /*for the column*/
    for (j = 0; i < a; i++)
    {
        for (k = 0; j < b; j++)
        {
            memmove(tempa3[i][j], tempa2[i][j] + (c - cc), sizeof(double)*cc);
            memmove(tempa3[i][j] + cc, tempa2[i][j], sizeof(double)*(c-cc));
        }
    }

但问题是我定义了太多新的动态数组，而且 tempa3 的结果也不正确。有人可以提出一些建议吗？

4

1 回答 1

0

我相信你想要这样的东西：

memmove(tempa, te + (a - aa), sizeof(double**) * aa);
memmove(tempa + aa, te, sizeof(double**) * (a - aa));

或者

memmove(tempa, te + aa, sizeof(double**) * (a - aa));
memmove(tempa + (a - aa), te, sizeof(double**) * aa);

取决于您是否要交换前半部分“向上或向下舍入”（我假设您希望将其向上舍入，这是第一个版本）。

不过，我不太喜欢您的代码设计：

首先，尽可能避免动态分配和std::vector使用std::array。你可能会争辩说它会阻止你安全地使用memmove而不是swap第一个维度（嗯，它应该可以工作，但我不是 100% 确定它不是实现定义的）但我认为这不会改善那么多效率。此外，如果你想要一个 N 维数组，我通常更喜欢避免“链接指针”（尽管使用你的算法，你实际上可以使用这个结构，所以它不是那么糟糕）。
例如，如果您坚持使用 动态分配数组new，则可以使用类似的方法来减少内存使用量（尽管差异可能可以忽略不计；它也可能稍微快一点，但同样，可能可以忽略不计）：

#include <cstddef>
#include <iostream>

typedef std::size_t index_t;

constexpr index_t width = 3;
constexpr index_t height = 4;
constexpr index_t depth = 5;

// the cells (i, j, k) and (i, j, k+1) are adjacent in memory
// the rows (i, j, _) and (i, j+1, _) are adjacent in memory
// the "slices" (i, _, _) and (i+1, _, _) are adjacent in memory
constexpr index_t cell_index(index_t i, index_t j, index_t k) {
  return (i * height + j) * depth + k;
}

int main() {
  int* array = new int[width * height * depth]();
  for( index_t i = 0 ; i < width ; ++i )
    for( index_t j = 0 ; j < height ; ++j )
      for( index_t k = 0 ; k < depth ; ++k ) {
        // do something on the cell (i, j, k)
        array[cell_index(i, j, k)] = i + j + k;
        std::cout << array[cell_index(i, j, k)] << ' ';
      }
  std::cout << '\n';
  // alternatively you can do this:
  //*
  for( index_t index = 0 ; index < width * height * depth ; ++index) {
    index_t i = index / (height * depth);
    index_t j = (index / depth) % height;
    index_t k = index % depth;
    array[index] = i + j + k;
    std::cout << array[index] << ' ';
  }
  std::cout << '\n';
  //*/
  delete[] array;
}

不同之处在于内存中的组织。在这里，您有一个 60*sizeof(int) 字节的大块（通常是 240 或 480 字节），而使用您的方法，您将拥有： - 1 个 3*sizeof(int**) 字节块 - 3 个 4*sizeof 块(int*) 字节 - 12 个 5*sizeof(int) 字节块（在 64 位架构上多 120 个字节，每个单元访问的两个额外间接，以及用于分配/释放所有内存的更多代码）当然，你可以' t 做 array[i][j][k] 了，但仍然......

与vectors 相同（您可以制作 anstd::vector<std::vector<std::vector<int>>>或 a std::vector<int>）

还有太多的代码重复：您的算法基本上将表的两半交换了三次（每个维度一次），但是您将相同的东西重写了 3 次，但有一些差异。

还有太多的内存分配/复制（您的算法有效，并且可以通过简单地交换指针来交换整个行/切片来利用指针数组的结构，在这种特定情况下，您可以利用此数据结构来避免使用您的算法进行复制...但你没有）

您应该选择更明确的变量名称，这会有所帮助。例如使用width, height,depth而不是a, b, c。

例如，这是一个带有向量的实现（虽然我不知道 matlab 的 fftshift 函数，但根据你的代码和这个页面，我认为它基本上是“交换角落”）：

（另外，使用 -std=c++11 编译）

#include <cstddef>
#include <iostream>
#include <vector>
#include <algorithm>

typedef std::size_t index_t;
typedef double element_t;
typedef std::vector<element_t> row_t;
typedef std::vector<row_t> slice_t;
typedef std::vector<slice_t> array_3d_t;

// for one dimension
// you might overload this for a std::vector<double>& and use memmove
// as you originally wanted to do here
template<class T>
void fftshift_dimension(std::vector<T>& row)
{
  using std::swap;
  const index_t size = row.size();
  if(size <= 1)
    return;
  const index_t halved_size = size / 2;
  // swap the two halves
  for(index_t i = 0, j = size - halved_size ; i < halved_size ; ++i, ++j)
    swap(row[i], row[j]);
  // if the size is odd, rotate the right part
  if(size % 2)
  {
    swap(row[halved_size], row[size - 1]);
    const index_t n = size - 2;
    for(index_t i = halved_size ; i < n ; ++i)
      swap(row[i], row[i + 1]);
  }
}

// base case
template<class T>
void fftshift(std::vector<T>& array) {
  fftshift_dimension(array);
}

// reduce the problem for a dimension N+1 to a dimension N
template<class T>
void fftshift(std::vector<std::vector<T>>& array) {
  fftshift_dimension(array);
  for(auto& slice : array)
    fftshift(slice);
}

// overloads operator<< to print a 3-dimensional array
std::ostream& operator<<(std::ostream& output, const array_3d_t& input) {
  const index_t width = input.size();
  for(index_t i = 0; i < width ; i++)
  {
    const index_t height = input[i].size();
    for(index_t j = 0; j < height ; j++)
    {
      const index_t depth = input[i][j].size();
      for(index_t k = 0; k < depth; k++)
        output << input[i][j][k] << ' ';
      output << '\n';
    }
    output << '\n';
  }
  return output;
}

int main()
{
  constexpr index_t width = 3;
  constexpr index_t height = 4;
  constexpr index_t depth = 5;

  array_3d_t input(width, slice_t(height, row_t(depth)));

  // initialization
  for(index_t i = 0 ; i < width ; ++i)
    for(index_t j = 0 ; j < height ; ++j)
      for(index_t k = 0 ; k < depth ; ++k)
        input[i][j][k] = i + j + k;
  std::cout << input;

  // in place fftshift
  fftshift(input);

  std::cout << "and then" << '\n' << input;
}

活生生的例子

您可以通过避免多次交换同一个单元格和/或使用 memmove 来制定一个稍微更有效的算法，但我认为它已经足够快，可以用于许多用途（在我的机器上，对于 1000x1000x100 表，fftshift 大约需要 130 毫秒）。

于 2015-08-13T18:58:51.493 回答