FFT 数据 C语言 C++

快速傅里叶变换毛毛虫法

亚历山大·谢苗诺夫

5.00/5 (13投票s)

2017 年 7 月 1 日

CPOL

10分钟阅读

22265

基于预编译工具的数据驱动方法开发最快的 FFT 实现

引言

随着处理需求的增加，在涉及DSP计算的应用程序中，性能正成为一个瓶颈。在大多数情况下，使用专用硬件并非一个可行的选择，因此问题需要一个软件解决方案。为了在将实现成本降至最低的同时实现更高性能，可以使用数据驱动的编程方法。在这种编程模型中，不同的数据集根据程序逻辑的通用形式控制程序的流程。通常有人声称，数据驱动方法过于狭窄，不适用于广泛使用的算法。然而，可以使用预编译工具或编译器扩展来克服这些担忧。本文描述了基于预编译工具构建两个数据驱动FFT算法的方法，并推广了所发现的“毛毛虫方法”作为广泛使用的蝶形方法的替代方案。

背景

通过数字傅里叶变换将时域信号转换为频域信号

X_k = _n=0^N-1∑ x_n* ω_N^kn, k= 0..N-1, ω_N= e^-2πi/N, i = ²√-1 ;

在此方案中，两个域都包含一个由N个复点组成的信号。这些复点x中的每一个都由两部分组成：实值和虚值。复数旋转因子ω表示一个“旋转向量”，它映射到复平面中的单位圆上。

实际上，DSP应用程序基于计算复杂度较低的快速傅里叶变换 (FFT) 算法。计算FFT的第一步是将当前的N点信号 X=(x₀, x₁_… x_n-1) 分成两部分 A=(x₀, x₂_… x_n-2) 和 B=(x₁, x₃_… x_n-1)，每部分包含 N/2 个点

A_k = _n=0^N/2-1∑ x_{2 n}* ω_N^k(2n) ;

B_k= _n=0^N/2-1∑ x_{2 n+1} * ω_N^k(2n+1) = _n=0^N/2-1∑ x_{2 n+1} * ω_N^k* ω_N^k(2n) = ω_N^k* _n=0^N/2-1∑ x_{2 n+1}* ω_N^k(2n);

因此，X_k = A_k + ω_N^k * B_k;

将偶数 N/2 点的 DFT 与奇数 N/2 点的 DFT 结合起来，生成一个 N 点的 DFT。在下一步中，N/2 点 DFT A_k 和 B_k 以相同的方式使用 N/4 点 DFT 计算。这种推导也适用于逆 FFT，这不在本文的讨论范围之内。

狭义上，标准乘累加操作（MAC）是 a = a + b * c，其中“a”是累加器。广义上，MAC可以表示为 c = a + w * b，这是FFT的基本操作。FFT的点数必须是2的幂。直接DFT实现需要N个点的N²次操作。对于FFT，总计算量与 **N*log₂(N)** 成比例。

使用 C++ 代码

根据上述公式直接实现FFT并不困难。FFT递归算法的伪代码可以在[1]中找到。为了验证这个概念，该算法被转换为现代C++代码，如下所示：

#include <complex>
#include <vector>
#include <algorithm>
#include <iostream>
#include <math.h>

#define M_PI 3.14159265358979323846
using namespace std;
typedef complex<double> w_type; 

static vector<w_type> fft(const vector<w_type> &In) 
{ 
  int i = 0, wi = 0;
  int n = In.size();
  vector<w_type> A(n / 2), B(n / 2), Out(n);
  if (n == 1) {
      return vector<w_type>(1, In[0]);
  }
  i = 0;
  copy_if( In.begin(), In.end(), A.begin(), [&i] (w_type e) {
           return !(i++ % 2);
  } );
  copy_if( In.begin(), In.end(), B.begin(), [&i] (w_type e) {
           return (i++ % 2);
  } );

  vector<w_type> At = fft(A);
  vector<w_type> Bt = fft(B);

  transform(At.begin(), At.end(), Bt.begin(), Out.begin(), [&wi, &n]
  (w_type& a, w_type& b) { 
       return  a + b * exp(w_type(0, 2 * M_PI * wi++ / n)); 
  });
  transform(At.begin(), At.end(), Bt.begin(), Out.begin() + n / 2, [&wi, &n]
   (w_type& a, w_type& b) { 
      return  a + b * exp(w_type(0, 2 * M_PI * wi++ / n)); 
  });
  return Out;
}
void main(int argc, char* argv[])
{
    int ln = (int)floor( log(argc - 1.0) / log(2.0) ); 
    vector<w_type> In(1 << ln);
    std::transform(argv + 1, argv + argc, In.begin(),[&](const char* arg) {
        return w_type(atof(arg), 0);
    });
    vector<w_type> Out = fft(In);
    for (vector<w_type>::iterator itr = Out.begin(); itr != Out.end(); itr++) {
        cout  << *itr << endl;
    }
}

在主程序中，输入程序参数的数量被截断为最接近2的幂，这定义了变换大小和输入向量的大小。标准C++复数模板库用于这些计算。STL复数向量用于处理。输入向量复数值的实部由输入程序参数组成。虚部用零填充，表示时域信号。

第一级递归将n点信号`In`分成两个信号`A`和`B`，每个信号包含n/2点。下一级将数据分成4个n/4点信号。当只剩下1点信号时，递归停止。变换后的`At`和`Bt`向量用于组成输出向量`Out`。

最后，程序从输出向量打印频域信号的实部和虚部。

不同风格的实现

C++ 实现的抽象级别过高，因此可能不适用于低成本平台。这是提供 C 语言实现的原因。然而，问题在于在可移植性和效率之间取得正确的平衡。软件工程师需要根据应用程序的要求应用不同的编程风格。代码效率的平台优化可以适用于嵌入式应用程序。应用的编程风格级别越高，所需的编程工作就越少。一般来说，代码效率与此呈反比关系。

MAC是FFT的基本操作，在复数表示中可以表示为 **c = a + w * b**。以下示例演示了不同编程风格在将复数MAC操作实现为`w_mac`函数时的适用性。

表 1. MAC 函数示例
原型	void w_mac(w_type* cc, w_type a, w_type w, w_type b)
高级	cc = a + b exp(w_type(0, 2 * M_PI * w / n)
经典	cc->r = a.r + w.r * b.r - w.i * b.i; cc->i = a.i + w.r * b.i + w.i * b.r;
嵌入式（已为平台 MAC 优化准备）	register double reg; reg = mac(a.r, w.r, b.r); cc->r = mac(reg, -w.i, b.i); reg = mac(a.i, w.r, b.i ); cc->i = mac(reg, w.i * b.r );

FFT大小的以2为底的对数是FFT递归的次数。下面提供了几种用于此类计算的n2ln实用例程的样式。

表 1. Log₂ 函数示例
原型	int n2ln( int n )
高级	return (int)floor( log(n) / log(2.0) );
经典在最坏情况下需要 15 个周期	while (n >>= 1) ln++; return ln;
嵌入式最坏情况等于平均情况每个结果需要 3 次汇编器移位、4 次分支和 1 次移动可以实现为宏，以便在编译时计算常量值	return n/256 ?n/4048 ?n/16348 ?n/32768?15:14 :n/8096?13:12 :n/1024 ?n/2048?11:10 :n/512?9:8 :n/16 ?n/64 ?n/128?7:6 :n/32?5:4 :n/4 ?n/8?3:2 :n/2?1:0;

经典 C 语言实现

此 C 语言实现在算法上与上面的 C++ 实现类似。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define M_PI 3.14159265358979323846

typedef struct { double r; double i; } w_type;

int n2ln( int n );
void w_mac( w_type* cc, w_type a, w_type w, w_type b );

static void fft0( w_type InOut[], int n ) 
{ 
    int i;
    w_type w, *A, *B;

    if (n == 1) return; 
    A = malloc( sizeof(w_type) * n / 2 );
    B = malloc( sizeof(w_type) * n / 2 );
    for (i = 0; i < n / 2; i++) {
        A[i] = InOut[i * 2];
        B[i] = InOut[i * 2 + 1];
    }
    fft0( A, n / 2 );
    fft0( B, n / 2 );
    for (i = 0; i < n; i++) {
        w.r = cos(2 * M_PI * i / n);
        w.i = sin(2 * M_PI * i / n);
        w_mac( &InOut[i], A[i % (n / 2)], w, B[i % (n / 2)] );
    }
    free( A );
    free( B );
}
void main( int argc, char * argv[] )
{
    int i;
    int ln = n2ln(argc - 1);
    w_type* InOut = malloc( (1 << ln) * sizeof(w_type) ); 
    for (i = 0; i < (1 << ln); i++) {
        InOut[i].r = atof(argv[i+1]); 
        InOut[i].i = 0;
    }
    fft0( InOut, 1 << ln );
    for(i = 0; i < (1 << ln); i++) {
        printf("%.4f %.4f\n", InOut[i].r, InOut[i].i);
    }
}

经典 C 语言示例引入了几处算法修改

复数计算手动完成
由于显式内存分配，实际输入缓冲区大小被传递给递归过程
C++ 实现的输入和输出向量被组合在一起，并作为复数值的单个 `InOut` 缓冲区共享。

数据驱动的修改

以下示例基于经典 C 语言实现，旨在演示如何构建一个可创建嵌入式风格程序的预编译工具。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define M_PI 3.14159265358979323846

#define LN_FFT 4    /* number of samples is 1 << LN_FFT */

typedef struct { double r; double i; } w_type;

int n2ln( int n ); 
void w_mac( w_type* cc, w_type a, w_type w, w_type b );

static struct tMac {
    w_type *c, *a, *b, w;
} Mac[LN_FFT * (1 << LN_FFT)];
static int nMac;
static w_type DataTrace[LN_FFT + 1][1 << LN_FFT];
static int BusyDataTrace[LN_FFT + 1];

static void calculate_macs()
{
    int i;
    for (i = 0; i < nMac; i++) { 
        w_mac(Mac[i].c, *Mac[i].a, Mac[i].w, *Mac[i].b);
    }
}

static void record_mac( w_type** cc, w_type* a, w_type w, w_type *b, int n )
{
    int ln = n2ln(n);
    int i = BusyDataTrace[ln]++;

    *cc = &DataTrace[ln][i]; 
    Mac[nMac].c = &DataTrace[ln][i];
    Mac[nMac].w = w;
    Mac[nMac].a = a;
    Mac[nMac].b = b;
    nMac++;
}

static void fft0( w_type* InOut[], int n ) 
{ 
    int i;
    w_type w, **A, **B;

    if (n == 1) return;  
    A = malloc( sizeof(w_type*) * n / 2 );
    B = malloc( sizeof(w_type*) * n / 2 );
    for (i = 0; i < n / 2; i++) {
        A[i] = InOut[i * 2];
        B[i] = InOut[i * 2 + 1];
    }
    fft0( &A[0], n / 2 );
    fft0( &B[0], n / 2 );
    for (i = 0; i < n; i++) {
        w.r = cos(2 * M_PI * i / n);
        w.i = sin(2 * M_PI * i / n);
        record_mac( &InOut[i], A[i % (n / 2)], w, B[i % (n / 2)], n );
    }
    free(A); 
    free(B);
}

void main( int argc, char* argv[] )
{
    int i;
    w_type** InOut = malloc( sizeof(w_type*) * (1 << LN_FFT) );
    for (i = 0; i < (1 << LN_FFT); i++) {
        InOut[i] = &DataTrace[0][i];
        DataTrace[0][i].r = atof( argv[i % (argc - 1) + 1] );  
        DataTrace[0][i].i = 0;
    }
    fft0( InOut, 1 << LN_FFT );
    calculate_macs();
    for(i = 0; i < (1 << LN_FFT); i++) {
        printf("%.4f %.4f\n", DataTrace[LN_FFT][i].r, DataTrace[LN_FFT][i].i);
    }
    free(InOut);
}

这个FFT递归算法已经重写，采用了数据驱动的方法进行FFT计算。递归过程已经重新基于支持复杂类型指针的缓冲区，而不是直接使用复杂类型。实际计算被记录到 `Mac` 数组中，以便稍后计算。二维数组 `DataTrace` 用于跟踪所有计算。当FFT的所有递归完成后，用户需要调用 `calculate_macs`。所有重新编码的复数mac都按照它们被重新编码的顺序进行计算。 `calculate_macs` 过程只有一个循环，因此它尽可能快。该算法具有理论FFT复杂度n*log₂(n)，现在这也与内存使用相关。这意味着 `Mac` 和 `DataTrace` 数组都具有n*log₂(n)个元素。这也可以适用于效率至上的现代计算系统，但低成本系统将需要额外的修改。

基于生成代码的表格实现

对于跨平台系统，数据驱动方法转变为静态配置方法。这意味着上述示例中的 Mac 数组被写入为 C 语言初始化的结构，作为新程序构建和运行。Mac 数组元素的转置有两个自由度

相同递归级别的 `Mac` 数组元素可以安全地互换
Mac 数组元素的指针可以更改以利用其他 DataTrace 元素。当然，需要注意引用更改位置的其他指针。

这使得所需的优化目标得以实现。上一个示例中的 Mac 数组利用了来自 `DataTrace` 数组的 N*log₂(N) 个 RAM 元素。之前的程序已修改，将 RAM 使用优化到 **N+2** 个元素，并生成以下代码：

#include <stdlib.h>
#include <stdio.h>

#define LN_FFT  4  /* power of 2 is number of fft points */
/*   */
#define W_0_02   1.000000000000000 /* angle 0.00 dg */ 
#define W_1_04   0.000000000000000 /* angle 90.00 dg */ 
#define W_1_08   0.707106781186548 /* angle 45.00 dg */ 
#define W_1_16   0.923879532511287 /* angle 22.50 dg */ 
#define W_3_16   0.382683432365090 /* angle 67.50 dg */ 

typedef struct { double r; double i; } w_type;
static const struct fft_tbl {
    double wr, wi;
    unsigned char c, a, b;
} tbl[] = {
{ W_0_02,+W_1_04,16, 0, 8}, {-W_0_02,+W_1_04,17, 0, 8},
{ W_0_02,+W_1_04, 0, 4,12}, {-W_0_02,+W_1_04, 8, 4,12},
{ W_0_02,+W_1_04, 4, 2,10}, {-W_0_02,+W_1_04,12, 2,10},
{ W_0_02,+W_1_04, 2, 6,14}, {-W_0_02,+W_1_04,10, 6,14},
{ W_0_02,+W_1_04, 6, 1, 9}, {-W_0_02,+W_1_04,14, 1, 9},
{ W_0_02,+W_1_04, 1, 5,13}, {-W_0_02,+W_1_04, 9, 5,13},
{ W_0_02,+W_1_04, 5, 3,11}, {-W_0_02,+W_1_04,13, 3,11},
{ W_0_02,+W_1_04, 3, 7,15}, {-W_0_02,+W_1_04,11, 7,15},
{ W_0_02,+W_1_04, 7,16, 0}, {-W_0_02,+W_1_04,15,16, 0},
{ W_1_04,+W_0_02,16,17, 8}, {-W_1_04,-W_0_02, 0,17, 8},
{ W_0_02,+W_1_04,17, 4, 2}, {-W_0_02,+W_1_04, 8, 4, 2},
{ W_1_04,+W_0_02, 4,12,10}, {-W_1_04,-W_0_02, 2,12,10},
{ W_0_02,+W_1_04,12, 6, 1}, {-W_0_02,+W_1_04,10, 6, 1},
{ W_1_04,+W_0_02, 6,14, 9}, {-W_1_04,-W_0_02, 1,14, 9},
{ W_0_02,+W_1_04,14, 5, 3}, {-W_0_02,+W_1_04, 9, 5, 3},
{ W_1_04,+W_0_02, 5,13,11}, {-W_1_04,-W_0_02, 3,13,11},
{ W_0_02,+W_1_04,13, 7,17}, {-W_0_02,+W_1_04,11, 7,17},
{ W_1_08,+W_1_08, 7,16, 4}, {-W_1_08,-W_1_08,17,16, 4},
{ W_1_04,+W_0_02,16,15, 8}, {-W_1_04,-W_0_02, 4,15, 8},
{-W_1_08,+W_1_08,15, 0, 2}, { W_1_08,-W_1_08, 8, 0, 2},
{ W_0_02,+W_1_04, 0,12,14}, {-W_0_02,+W_1_04, 2,12,14},
{ W_1_08,+W_1_08,12, 6, 5}, {-W_1_08,-W_1_08,14, 6, 5},
{ W_1_04,+W_0_02, 6,10, 9}, {-W_1_04,-W_0_02, 5,10, 9},
{-W_1_08,+W_1_08,10, 1, 3}, { W_1_08,-W_1_08, 9, 1, 3},
{ W_0_02,+W_1_04, 1,13, 0}, {-W_0_02,+W_1_04, 3,13, 0},
{ W_1_16,+W_3_16,13, 7,12}, {-W_1_16,-W_3_16, 0, 7,12},
{ W_1_08,+W_1_08, 7,16, 6}, {-W_1_08,-W_1_08,12,16, 6},
{ W_3_16,+W_1_16,16,15,10}, {-W_3_16,-W_1_16, 6,15,10},
{ W_1_04,+W_0_02,15,11, 2}, {-W_1_04,-W_0_02,10,11, 2},
{-W_3_16,+W_1_16,11,17,14}, { W_3_16,-W_1_16, 2,17,14},
{-W_1_08,+W_1_08,17, 4, 5}, { W_1_08,-W_1_08,14, 4, 5},
{-W_1_16,+W_3_16, 4, 8, 9}, { W_1_16,-W_3_16, 5, 8, 9},
};

static const unsigned char OutOrder[]={
1,13,7,16,15,11,17,4,3,0,12,6,10,2,14,5,};

static struct { double r; double i; } XY[(1 << LN_FFT) + 2];

void fft_table()
{
    int i;
    register const struct fft_tbl* t;

    for (i = 0, t = tbl; i < sizeof(tbl) / sizeof(tbl[0]); i++, t++) {
        XY[t->c].r = XY[t->a].r + t->wr * XY[t->b].r - t->wi * XY[t->b].i;
        XY[t->c].i = XY[t->a].i + t->wr * XY[t->b].i + t->wi * XY[t->b].r;
    }
}
void main(int argc, char* argv[])
{
    int i;
    for (i = 0; i < (1 << LN_FFT); i++) {
        XY[i].r = atof( argv[i % (argc - 1) + 1] );  
        XY[i].i = 0;
    }
    fft_table();
    for(i = 0; i < (1 << LN_FFT); i++) {
        printf("%.4f %.4f\n", XY[OutOrder[i]].r, XY[OutOrder[i]].i);
    }
}

这是一个16点FFT的例子。`tbl` 数组的一个元素包含一个旋转因子的复数值和3个偏移量，以基于XY内存数组[1]的3个RAM元素提供一个复数MAC操作。该代码速度快，因为它只使用一个“`for`”循环。

一个FFT复数运算由4个浮点运算组成，因此可以生成一个用于浮点或定点运算的数据驱动表，以提供优化能力，实现最高效率。这类矩阵或算法形式的优化在参考文献[2]中广为人知，但现在可以在所提出的框架基础上进行描述。

基于特定 FFT 特征的优化[2]
- 输入时域信号由N个复数点信号表示，其虚部为零
- 输出频域信号具有复制方案。点N/2+1的实部等于点N/2-1的实部。这一直应用到点N-1的实部，它与点1的实部相同。点N/2+1的虚部是点N/2-1虚部的负值。样本0和N/2没有匹配点。因此，点N/2+1到N-1是冗余的，输出频域信号可以由点0到N/2充分表示。

这些优化可以达到 N*log₂(N - 1/2 - 1/2) = **N*log₂(N)/2** 的性能。

基于退化旋转因子的优化
- 旋转因子为0、1和-1的运算可以简化为简单的加法；
- 旋转因子为0.7071/* 45.00度角 */ 的运算（实部和虚部相同）可以节省一个MAC操作。
输入和输出表示的优化
- 上面的例子使用 `OutOrder` 数组来恢复频域信号样本的标准FFT顺序。事实上，这种顺序很少用于后续处理。该工具可以生成具有所需自定义顺序的 `OutOrder` 数组。
- 通常，在进行 FFT 之前，时域信号可以进行加窗和/或归一化。此操作可以引入到 FFT 表方法中，以便由该工具支持。

[1] 本文描述了数据驱动方法的通用原理，并未考虑任何常量表优化。对于大多数架构而言，一个包含5个元素的结构化数组不如5个单独的数组有效。

[2] 逆FFT的优化有所不同。

基于生成代码的毛毛虫实现

FFT表方法利用一个包含 N*log₂(N) 个元素的表。另一种优化方法是组合在某些字段中相似的表元素。以下示例演示了根据旋转因子对复数 MAC 进行分组，以将它们生成为代码语句

#include <stdlib.h>
#include <stdio.h>

#define LN_FFT  5  /* power of 2 is number of fft points */
/*   */
#define W_0_02   1.000000000000000 /* angle 0.00 dg */ 
#define W_0_04   0.000000000000000 /* angle 90.00 dg */ 
#define W_0_08   0.707106781186547 /* angle 45.00 dg */ 
#define W_0_16   0.923879532511287 /* angle 22.50 dg */ 
#define W_1_16   0.382683432365090 /* angle 67.50 dg */ 
#define W_0_32   0.980785280403230 /* angle 11.25 dg */ 
#define W_1_32   0.831469612302545 /* angle 33.75 dg */ 
#define W_2_32   0.555570233019602 /* angle 56.25 dg */ 
#define W_3_32   0.195090322016128 /* angle 78.75 dg */ 

typedef struct { double r; double i; } w_type;

#define X2Y(a) (a + (1 << (LN_FFT - 1))) 
#define XMAC(c, a, wr, wi) \
    c->r = a->r + wr * X2Y(a)->r - wi * X2Y(a)->i; \
    c->i = a->i + wr * X2Y(a)->i + wi * X2Y(a)->r;

static w_type XY[2][(1 << LN_FFT)];
static const w_type* pOut = LN_FFT % 2 ? &XY[1][0] : &XY[0][0];
static const unsigned char OutOrder[]={31,15,23,14,27,13,22,12,29,11,21,
10,26,9,20,8,30,7,19,6,25,5,18,4,28,3,17,2,24,1,16,0,};

void fft_caterpillar()
{
    int i, j, lim;
    register  w_type *pc, *pa; /* pb = a + (1 << (LN_FFT - 1)) */
    for (i = 1; i <= LN_FFT; i++) {
        pc = i % 2 ? &XY[1][0] : &XY[0][0];
        pa = i % 2 ? &XY[0][0] : &XY[1][0];
        lim = 1 << (LN_FFT - i);
        for (j = 0; j < lim; j++) {
            switch (i) {
            case 5:
                XMAC(pc, pa,  W_0_32, -W_3_32); pc++; pa += 1;
                XMAC(pc, pa,  W_1_32, -W_2_32); pc++; pa += 1;
                XMAC(pc, pa,  W_2_32, -W_1_32); pc++; pa += 1;
                XMAC(pc, pa,  W_3_32, -W_0_32); pc++; pa += 1;
                XMAC(pc, pa, -W_3_32, -W_0_32); pc++; pa += 1;
                XMAC(pc, pa, -W_2_32, -W_1_32); pc++; pa += 1;
                XMAC(pc, pa, -W_1_32, -W_2_32); pc++; pa += 1;
                XMAC(pc, pa, -W_0_32, -W_3_32); pc++; pa += 1;
                pa -= 8;
                XMAC(pc, pa, -W_0_32, +W_3_32); pc++; pa += 1;
                XMAC(pc, pa, -W_1_32, +W_2_32); pc++; pa += 1;
                XMAC(pc, pa, -W_2_32, +W_1_32); pc++; pa += 1;
                XMAC(pc, pa, -W_3_32, +W_0_32); pc++; pa += 1;
                XMAC(pc, pa,  W_3_32, +W_0_32); pc++; pa += 1;
                XMAC(pc, pa,  W_2_32, +W_1_32); pc++; pa += 1;
                XMAC(pc, pa,  W_1_32, +W_2_32); pc++; pa += 1;
                XMAC(pc, pa,  W_0_32, +W_3_32); pc++; pa += 1;
            case 4:
                XMAC(pc, pa,  W_0_16, -W_1_16); pc++; pa += 1;
                XMAC(pc, pa,  W_1_16, -W_0_16); pc++; pa += 1;
                XMAC(pc, pa, -W_1_16, -W_0_16); pc++; pa += 1;
                XMAC(pc, pa, -W_0_16, -W_1_16); pc++; pa += 1;
                pa -= 4;
                XMAC(pc, pa, -W_0_16, +W_1_16); pc++; pa += 1;
                XMAC(pc, pa, -W_1_16, +W_0_16); pc++; pa += 1;
                XMAC(pc, pa,  W_1_16, +W_0_16); pc++; pa += 1;
                XMAC(pc, pa,  W_0_16, +W_1_16); pc++; pa += 1;
            case 3:
                XMAC(pc, pa,  W_0_08, -W_0_08); pc++; pa += 1;
                XMAC(pc, pa, -W_0_08, -W_0_08); pc++; pa += 1;
                pa -= 2;
                XMAC(pc, pa, -W_0_08, +W_0_08); pc++; pa += 1;
                XMAC(pc, pa,  W_0_08, +W_0_08); pc++; pa += 1;
            case 2:
                XMAC(pc, pa, -W_0_04, -W_0_02); pc++; pa += 1;
                pa -= 1;
                XMAC(pc, pa,  W_0_04, +W_0_02); pc++; pa += 1;
            case 1:
                XMAC(pc, pa, -W_0_02, +W_0_04); pc++; pa += 1;
                pa -= 1;
            case 0:
                XMAC(pc, pa,  W_0_02, +W_0_04); pc++; pa += 1;
            }
        }
    }
}
void main(int argc, char* argv[])
{
    int i;
    for (i = 0; i < (1 << LN_FFT); i++) {
        XY[0][i].r = atof( argv[i % (argc - 1) + 1] );  
        XY[0][i].i = 0;
    }
    fft_caterpillar();
    for(i = 0; i < (1 << LN_FFT); i++) {
        printf("%.4f %.4f\n", pOut[OutOrder[i]].r, pOut[OutOrder[i]].i);
    }
}

这是一个32点FFT的例子。复数MAC的数量是N。二维数组XY以交换方案组织，代表4个RAM内存库。每个`XMAC`同时访问3个内存库：一个用于读取，两个用于写入。该程序是为DSP平台准备的。

这种方法的主要优点是简化了地址算术。构建有效FFT算法的传统方法是在数据数组上使用复杂的地址算术构造循环。这种算术的图形表示看起来像蝴蝶翅膀。在这个例子中，蝴蝶转置被分解了，所以它是*毛毛虫方法*。长开关操作符也保留了毛毛虫的特征。

进一步的机器优化

表方法部分中描述的几乎所有标准优化都可以应用于毛毛虫方法。毛毛虫方法的主要优点是能够将其优化为一组抽象的 MAC 指令。上面的 `XMAC` 语句实现了 c = a + ω * b 复数方程。它可以由以下 4 个 MAC 表示：

reg = a.real + w.real * b.real;
c.real = reg - w.imag * b.imag;
reg = a.imag + w.real * b.imag;
c.imag = reg + w.imag * b.real;

实际上，XMAC 可以通过 Intel AVX 内在函数表示为

#define XMAC(c, a, b, wr, wi) \
   vec1 = _mm_setr_pd(wr, wi); \
   vec2 = *( __m128d*)b; \
   vec1 = _mm_hadd_pd(_mm_mul_pd(vec1, _mm_mul_pd(vec2, neg)), \
           _mm_mul_pd(vec1, _mm_shuffle_pd(vec2, vec2, 1))); \
   *( __m128d*)c = _mm_add_pd(*( __m128d*)a, vec1);

所提供的毛毛虫方法示例同时操作两个 64 位双字。对于 AVX3，四个最近的 XMAC 可以合并为一个语句，以类似的方式操作八个 64 位双字。对于 32 位浮点或定点实现，这可能会更多。

结论

一个精心设计的工具可以根据客户的具体需求生成快速的数据驱动算法实现。这些算法实现可以为编译器优化做好准备，或者它们可以基于汇编代码。一般来说，ROM 中的大型常量表对于低成本平台是可以接受的，因此可以实现最高效率。另外，FFT 毛毛虫方法可以适用于大多数 DSP 架构。这种方法可以扩展到广泛的滤波和变换算法。

参考文献

[1] 递归FFT（快速傅里叶变换）算法笔记，2010年秋季，COSC 511，Susan Haynes http://www.emunix.emich.edu/~haynes/Papers/FFT/NotesOnRecursiveFFT.pdf

[2] FFT 笔记，C. S. Burrus，休斯顿莱斯大学电气与计算机工程系，邮编：77251-1892，http://www.jjj.de/fft/fftnote.txt

[3] 一些 fft_caterpillar 示例，https://github.com/r35382/fft_caterpillar