c++ - 如何将递归代码更改为迭代形式

Question

我在下面实现了一个简单的 FFT(忽略最终缩放):

typedef complex<double> base;
vector<base> w;
int FFTN = 1024;
void fft(vector<base> &fa){
    int n = fa.size();
    if (n==1) return;
    int half = (n>>1);
    vector<base> odd(half),even(half);
    for(int i=0,j = 0;i<n;i+=2,j++) {
        even[j] = fa[i];
        odd[j] = fa[i+1];    
    }    
    fft(odd);
    fft(even);    
    int fact = FFTN/n;    
    for (int i=0;i<half;i++){        
        fa[i] = even[i] + odd[i] * w[i * fact];
        fa[i + half] = even[i] - odd[i] * w[i * fact];
    }
}

效果很好。但我一直坚持将其转换为迭代形式。到目前为止我已经尝试过:

int n = fa.size();
int fact = (FFTN>>1);
int half = 1;
while(half<n){
    for(int i=0;i<n/half;i+=2){
        base even = fa[i], odd = fa[i+1];
        fa[i] = even + odd * w[i*fact];
        fa[i+half] = even - odd*w[i*fact];                            
    }
    for(int j=0;j<n/half;j++) 
        fa[j] = fa[j+half];
    fact >>= 1;
    half <<= 1;
}

有人可以帮我解决转换方面的技巧吗？

Answer 1

我要做的第一件事就是让你的函数“更加递归”。

void fft(base* fa, size_t stride, size_t n) {
  if (n==1) return;
  int half = (n>>1);
  fft(fa+stride, stride*2, half); // odd
  fft(fa, stride*2, half); // even
  int fact = FFTN/n;  
  for (int i=0;i<half;i++){    
    fa[i] = fa[stride*2*i] + fa[stride*2+i+stride] * w[i * fact];
    fa[i + half] = fa[stride*2*i] - fa[stride*2+i+stride] * w[i * fact];
  }
}
void fft(std::vector<base>& fa){ fft(fa.data(), 1, fa.size()); }

现在我们在缓冲区内就地进行fft。

由于我们现在有两个不同的 fft 实现，因此我们可以对它们进行相互测试。此时构建一些单元测试，以便可以针对已知的“良好”(或至少稳定)行为集来测试进一步的更改。

接下来，我们可以检查元素在原始 vector 中的组合顺序。检查长度为 4 的缓冲区。

a   b   c   d

我们递归执行奇数和偶数

a[e]  b[o]  a[e]  d[o]

然后递归执行奇数和偶数

a[ee]  b[oe]  a[eo]   d[oo]

这些集合的大小为 1。将它们单独保留，然后我们将奇数和偶数组合起来。

现在我们看 8。两次递归后，元素“拥有”于:

0[ee]   1[oe]   2[eo]   3[oo]   4[ee]   5[oe]   6[eo]  7[oo]

3之后:

0[eee]   1[oee]   2[eoe]   3[ooe]   4[eeo]   5[oeo]   6[eoo]  7[ooo]

如果我们反转这些标签，并调用 e 0 和 o 1，我们会得到:

0[000]   1[001]   2[010]   3[011]   4[100]   5[101]   6[110]   7[111]

这是二进制计数。第一位被丢弃，现在相等的元素在倒数第二个递归调用中组合。

然后丢弃前两位，并合并最后一位匹配的元素。

我们可以不查看位，而是查看每个组合的开头和步幅长度。

第一个组合的步幅长度等于数组长度(每个元素 1 个)。

第二个是长度/2。第三个是长度/4。

这一直持续到步幅长度为 1。

要组合的子数组数量等于步幅长度。

所以

for(size_t stride = n; stride = stride/2; stride!=0) {
  for (size_t offset = 0; offset != stride; ++offset) {
    fft_process( array+offset, stride, n/stride );
  }
}

其中 fft_process 基于:

  int fact = FFTN/n;  
  for (int i=0;i<half;i++){    
    fa[i] = fa[stride*2*i] + fa[stride*2+i+stride] * w[i * fact];
    fa[i + half] = fa[stride*2*i] - fa[stride*2+i+stride] * w[i * fact];
  }

也许是这样的:

void fft_process( base* fa, size_t stride, size_t n ) {
  int fact = FFTN/n; // equals stride I think!  Assuming outermost n is 1024.
  for (int i=0;i<half;i++){    
    fa[i] = fa[stride*2*i] + fa[stride*2+i+stride] * w[i * fact];
    fa[i + half] = fa[stride*2*i] - fa[stride*2+i+stride] * w[i * fact];
  }
}

这些都没有经过测试，但它给出了如何执行此操作的分步示例。您需要在此迭代版本上释放之前编写的单元测试(以测试 fft 的两个早期版本)。