65.9K
CodeProject 正在变化。 阅读更多。
Home

快速 FFT

Mosi_62
starIconstarIconstarIconstarIconstarIcon

5.00/5 (18投票s)

2013年7月11日

CPOL

4分钟阅读

viewsIcon

73962

downloadIcon

6671

一个比标准实现运行速度稍快的 FFT 算法。

基本 Radix-2-FFT 算法（递归）

源自傅里叶变换的主要公式

 

 

对于 N = 2^j 个样本，Radix-2-FFT 算法基于以下公式：

 

 

其中 M = N/2。

通过替换 u[m] = x[2m] 和 v[m] = x[2m + 1]，我们得到

 

 

在这个公式中，两个子 DFT

 

和

是可见的，这意味着基本的傅里叶变换可以分解成两个子变换，而这两个子变换又可以各自再分解成两个子子变换，以此类推，直到每个变换只剩下两个样本。

如果我们从 8 个样本开始进行第一次分割，看起来会是这样

如果我们继续分割，流程最终看起来会是这样

乍一看，它看起来相当复杂。输入数据是混杂的，并且有很多箭头指向各个方向。但是，如果我们只看包含一个步骤的第一个图形，我们会得到一个简单的系统。在上面的公式中，我们看到 [2m] 元素是偶数，而 [2m+1] 元素是奇数。这就是为什么奇数元素在这里与偶数元素分离。如果我们看一下第一个公式，并假设一个子变换的输出覆盖了该变换的输入变量，那么我们只需要执行子变换，然后计算复数和（假设 x[k] 是子变换的输出，而 X[k] 是这个变换的输出）

并且使用

 

我们可以替换

这就是第一个图形所展示的。

通过递归方法，可以将它放入一个相当小的算法中，该算法比非递归算法更容易理解。这里需要认识到的一个重要细节可能是，我们不需要所谓的位反转。我们只需要在每次递归中将奇数样本与偶数样本分离，正如在图形中看到的那样。通常在对 Radix-2 算法的描述中没有提到这一点。

当然，这里也可以进行位反转。它必须在算法的开始处完成……有很多不同的方法来实现这种转换。

递归过程看起来会是这样

public void CalcSubFFT(TKomplex[] a, int n, int lo)
{
    int i, m;
    TKomplex w;
    TKomplex v;
    TKomplex h;
    if (n > 1)
    {
        Shuffle(a, n, lo);  // separate odd from even elements
        m = n / 2;
        CalcSubFFT(a, m, lo);  //recursive call for u(k)
        CalcSubFFT(a, m, lo + m); //recursive call for v(k)
        for (i = lo; i < lo + m; i++)
        {
            w.real = Math.Cos(2.0 * Math.PI * (double)(i - lo) / (double)(n));
            w.imag = Math.Sin(2.0 * Math.PI * (double)(i - lo) / (double)(n));
            h = kprod(a[i + m], w);
            v = a[i];
            a[i] = ksum(a[i], h);
            a[i + m] = kdiff(v, h);
        }
    }
}

将奇数和偶数样本分开

public void Shuffle(TKomplex[] a, int n, int lo)
{
    if (n > 2)
    {
        int i, m = n / 2;
        TKomplex[] b = new TKomplex[m];
        for (i = 0; i < m; i++)
        b[i] = a[i * 2 + lo + 1];
        for (i = 0; i < m; i++)
        a[i + lo] = a[i * 2 + lo];
        for (i = 0; i < m; i++)
        a[i + lo + m] = b[i];
    }
}

最后，我们必须将所有这些放在一个函数中，并将每个元素除以 2N，并将元素 0 除以 4N（参见傅里叶变换的主要公式）。

public void CalcFFT()
{
    int i;
    CalcSubFFT(y, N, 0);
    for (i = 0; i < N; i++)
    {
        y[i].imag = y[i].imag / (double)N * 2.0;
        y[i].real = y[i].real / (double)N * 2.0;
    }
    y[0].imag = y[0].imag / 2.0;
    y[0].real = y[0].real / 2.0;
}

减少计算时间

正如我在我关于快速 DFT 的文章中已经发现的那样，我们也必须在这些算法中计算许多正弦和余弦值。如果我们想在一个应用程序中重复执行 FFT，那么将所有的正弦和余弦值放入一个查找表中，并从那里获取它们是有意义的。这大大加快了计算速度。

我们这样创建查找表

for (i = 0; i < (N / 2); i++)
    we[i].real = Math.Cos(2.0 * Math.PI * (double)(i) / (double)(N));
    we[i].imag = Math.Sin(2.0 * Math.PI * (double)(i) / (double)(N));
}

并替换以下行

w.real = Math.Cos(2.0 * Math.PI * (double)(i - lo) / (double)(n));
w.imag = Math.Sin(2.0 * Math.PI * (double)(i - lo) / (double)(n));

by

w.real = we[(i - lo)*N/n].real;
w.imag = we[(i - lo)*N/n].imag;

对算法的这种修改将我的计算机上 4096 个样本的计算时间从 5 毫秒减少到 3 毫秒。

在 DFT 算法中，我们必须执行取模运算才能获得正确的索引。这里没有必要，因为无论如何我们只有 0 到 2Pi 的正弦和余弦值。

基本 Radix-2-FFT 算法（非递归）

易于理解的递归方法的缺点是，并非所有环境都允许递归函数调用，并且时间消耗略高。总而言之，使用这些洗牌和递归需要完成更多的操作。因此，我也实现了一个非递归算法。它看起来有点复杂，但它是值得的。

在这里我们需要 Bitreversal 函数

public void BitInvert(TKomplex[] a, int n)
{
    int i, j, mv = n/2;
    int k, rev = 0;
    TKomplex b;
    for (i = 1; i < n; i++) // run tru all the indexes from 1 to n
    {
        k = i;
        mv = n / 2;
        rev = 0;
        while (k > 0) // invert the actual index
        {
            if ((k % 2) > 0)
            rev = rev + mv;
            k = k / 2;
            mv = mv / 2; 
        }
        {  // switch the actual sample and the bitinverted one
            if (i < rev)
            {
                b = a[rev];
                a[rev] = a[i];
                a[i] = b;
            }
        }
    }
}

对于每个索引“I”，通过“k”运行所有位并将它们反转。取模运算告诉我最低有效位是 1 还是 0，并且是否必须将其放入反转的索引中。从最高有效位开始。将 k 除以 2 将第二最低有效位向右移动一位。通过这种方式，我们从最低有效位到最高有效位获取所有位，并在“rev”中进行反转。现在我们只需要切换索引“i”和“rev”。但只有在“I”小于“rev”时才替换。否则，如果“I”变得大于“rev”，我们会将所有索引更改回。

在此位反转之后，我们可以通过以下方式计算子 FFT

public void CalcSubFFT(TKomplex[] a, int n)
{
    int i, k, m;
    TKomplex w;
    TKomplex v;
    TKomplex h;
    k = 1;
    while (k <= n/2)  // we start at stage 0 whit 2 sample per transformation
    {
        m = 0;
        while (m <= (n-2*k))
        {
            for (i = m; i < m + k; i++)
            {
                w.real = Math.Cos( Math.PI * (double)(i-m)/(double)(k));
                w.imag = Math.Sin( Math.PI * (double)(i-m)/(double)(k));
                h = kprod(a[i + k], w);
                v = a[i];
                a[i] = ksum(a[i], h);
                a[i + k] = kdiff(v, h);
            }
            m = m + 2 * k;
        }
        k = k * 2;
    }
}

这两个函数必须放在主调用中，并且所有元素必须除以 2N，并且元素 0 必须除以 4N。

public void CalcFFT()
{
    int i;
    BitInvert(y, N);
    CalcSubFFT(y, N);
    for (i = 0; i < N; i++)
    {
        y[i].imag = y[i].imag / (double)N * 2.0;
        y[i].real = y[i].real / (double)N * 2.0;
    }
    y[0].imag = y[0].imag / 2.0;
    y[0].real = y[0].real / 2.0;
}

减少计算时间

如果我们最终将正弦和余弦值放入一个查找表中

we = new TKomplex[N / 2];
for (i = 0; i < (N / 2); i++)  // Init look up table for sine and cosine values
{
    we[i].real = Math.Cos(2* Math.PI * (double)(i) / (double)(N));
    we[i].imag = Math.Sin(2* Math.PI * (double)(i) / (double)(N));
}

并替换

w.real = Math.Cos( Math.PI * (double)(i-m)/(double)(k));
w.imag = Math.Sin( Math.PI * (double)(i-m)/(double)(k));

by

w.real = we[((i-m)*N / k/ 2)].real;
w.imag = we[((i-m)*N / k / 2)].imag;

这也将转换的计算时间减少了大约 30%。只要我们只运行一次转换，就没有可见的与递归算法相比的时间消耗差异。但是如果我们运行 10 次，差异就会变得可见

10 次变换，每个变换 4096 个样本时的时间消耗比较

 

计算时间

递归标准

45 毫秒

带查找表的递归

30 毫秒

非递归标准

41 毫秒

带查找表的非递归

24 毫秒

更多数学知识请访问 www.mosismath.com :-)

© . All rights reserved.