如何实现 FFT 算法

• 下载源代码 - 26.9 Kb

引言

本文主要描述了一种针对复数样本数组实现一维 FFT 算法的方法。本文旨在展示一种高效且快速的 FFT 算法，用户可以根据自己的需求轻松修改。在我开发一个用于对捕获声音样本进行频率分析的软件时，我研究了 FFT 算法。

背景

首先，这段代码的 95% 不是我写的。这实际上是 1982 年《C 语言数值算法》一书中所描述的代码！！！是的，20 多年前了！！！但在我看来，它仍然非常非常好。但这与原始代码略有不同。当我研究该算法时，我注意到了一种可以利用的模式，并基于此，我通过代码中的微小更改成功改进了算法，使得算法的 O()（大 O 表示法）（衡量算法复杂度的单位）减少了 N 次计算。（成功实现改进方法后，我在网上进行了一些研究，发现我并没有发现任何新东西。我只是注意到了一些别人已经发现的东西 :-P）

我不会深入探讨“理论”，因此如果您想了解“为什么”，我强烈建议阅读第 12 章。书中还可以找到其他形式的 FFT，例如 2D 或 3D FFT。

FFT

快速傅里叶变换是一种优化的计算算法，用于对 2^N 个样本的数组实现离散傅里叶变换。它允许确定离散信号的频率，在频域中表示信号，进行卷积等。该算法的复杂度为 O(N*log2(N))。实际上，该算法的复杂度略高，因为数据需要通过一种称为位反转的操作进行准备。在《C 语言数值算法》中，位反转部分的复杂度为 O(2N)。通过我对这里提供的代码进行的微小更改，它将其复杂度降低到 O(N)。这表示性能提高了约 8%。

频域中信号的示例。

FFT 分两部分计算。第一部分通过应用位反转方法将原始数据数组转换为位反转顺序数组。这使得第二部分的数学计算“容易得多”。第二部分通过 N*log2(N) 次操作（Danielson-Lanzcos 算法的应用）处理 FFT。

让我们从一个复数数据数组开始。例如，在这种情况下，这个数组可以是浮点数数组，其中 data[偶数索引] 是实部，data[奇数索引] 是虚部。（这可以适应实数数据数组，只需将复数部分填充为 0，或使用书中实现的实数数组 FFT。）数组的大小必须是 N^2 阶（2、4、8、16、32、64 等）。如果样本大小不匹配，只需将其放入下一个 2^N 大小的数组中，并用 0 填充剩余空间。

一个小而不重要的考虑：原始代码使用的数据数组认为信息的开头在索引 1 -> data[1]，而 data[0] 被忽略。我的代码将其修改为从 0 -> data[0] 开始。

首先，我们定义 FFT 函数

//data -> float array that represent the array of complex samples
//number_of_complex_samples -> number of samples (N^2 order number) 
//isign -> 1 to calculate FFT and -1 to calculate Reverse FFT
float FFT (float data[], unsigned long number_of_complex_samples, int isign)
{
    //variables for trigonometric recurrences
    unsigned long n,mmax,m,j,istep,i;
    double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;

位反转方法

首先，必须对原始数组进行转换，以便执行 Danielson-Lanzcos 算法。例如，complex[索引] 将与 complex[位反转顺序索引] 交换位置。如果索引（二进制）是 0b00001，则位反转顺序索引将是 0b10000。在图 1 中，您可以看到转换后数据数组的变化。

根据《C 语言数值算法》实现此方法如下

    //the complex array is real+complex so the array 
    //as a size n = 2* number of complex samples
    // real part is the data[index] and the complex part is the data[index+1]
    n=number_of_complex_samples * 2; 

    //binary inversion (note that 
    //the indexes start from 1 witch means that the
    //real part of the complex is on the odd-indexes
    //and the complex part is on the even-indexes
    j=1;
    for (i=1;i<n;i+=2) { 
        if (j > i) {
            //swap the real part
            SWAP(data[j],data[i]); 
            //swap the complex part
            SWAP(data[j+1],data[i+1]);
        }
        m=n/2;
        while (m >= 2 && j > m) {
            j -= m;
            m = m/2;
        }
        j += m;
    }

SWAP 在函数外部，它大致是这样的

#defineSWAP(a,b)tempr=(a);(a)=(b);(b)=tempr
//tempr is a variable from our FFT function

图 1（8 个长度的数据数组）

如果你仔细观察图1，你会发现一个模式。让我们看看：将数组用一面镜子分成两半。如果你看着镜子的反射面，你会看到与另一边完全相同的东西。这种镜像效应允许你在数组的前半部分执行位反转方法，并几乎直接在后半部分使用它。但你必须小心一件事。只有当更改只发生在前一半时，你才能应用这种效应。这意味着如果更改发生在前半部分的某个索引和后半部分的某个索引之间，这是无效的，否则你将进行交换然后再次撤销（在16个长度的数据数组上执行此操作，你就会明白我在说什么）。所以代码变成这样

    //the complex array is real+complex so the array 
    //as a size n = 2* number of complex samples
    // real part is the data[index] and 
    //the complex part is the data[index+1]
    n=number_of_complex_samples * 2; 

    //binary inversion (note that the indexes 
    //start from 0 witch means that the
    //real part of the complex is on the even-indexes 
    //and the complex part is on the odd-indexes
    j=0;
    for (i=0;i<n/2;i+=2) {
        if (j > i) {
            //swap the real part
            SWAP(data[j],data[i]);
            //swap the complex part
            SWAP(data[j+1],data[i+1]);
            // checks if the changes occurs in the first half
            // and use the mirrored effect on the second half
            if((j/2)<(n/4)){
                //swap the real part
                SWAP(data[(n-(i+2))],data[(n-(j+2))]);
                //swap the complex part
                SWAP(data[(n-(i+2))+1],data[(n-(j+2))+1]);
            }
        }
        m=n/2;
        while (m >= 2 && j >= m) {
            j -= m;
            m = m/2;
        }
        j += m;
    }

Danielson-Lanzcos

代码的后半部分与书中描述的完全一致。这会将 N*log2(N) 次三角递归应用于数据。我将在这里展示的代码是我的版本（数据从索引 0 开始）

    //Danielson-Lanzcos routine 
    mmax=2;
    //external loop
    while (n > mmax)
    {
        istep = mmax<<  1;
        theta=sinal*(2*pi/mmax);
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        //internal loops
        for (m=1;m<mmax;m+=2) {
            for (i= m;i<=n;i+=istep) {
                j=i+mmax;
                tempr=wr*data[j-1]-wi*data[j];
                tempi=wr*data[j]+wi*data[j-1];
                data[j-1]=data[i-1]-tempr;
                data[j]=data[i]-tempi;
                data[i-1] += tempr;
                data[i] += tempi;
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }

如何使用 FFT

现在让我们看看如何使用 FFT。想象一下，您想收集一个信号或函数的样本，并想知道它的基频。这个样本可能来自任何来源：您插入代码中的函数、一段捕获的声音等。

假设信号是实数数组信号（就像声音捕捉缓冲区），我如何使用 FFT??? 首先，您需要选择要使用的 FFT 变体。有一个专门用于实数数组的变体，但在本例中，我将使用此变体。它不是最有效的，但易于使用。

接下来要考虑的是要发送到 FFT 的数据量和采样率。采样率必须是 2^N 的数字，但您不需要发送 2^N 个样本的数组进行处理（阅读 NR 以了解不同的实现）。您可以只发送 50 个样本，例如，并用 0 填充剩余数组。但请记住，您发送的计算数据越多，FFT 就越精确。

实数数组已转换为复数数组，其中复数部分等于 0 后，即可计算 FFT。

现在来看结果

FFT 计算完成后，您可以使用 FFT 得到的复数数组来得出结论。

如果您想知道信号的基频，请找到数组的绝对最大值，频率将由数组的索引给出。复数的绝对值是实部平方加虚部平方的平方根。

如果绝对最大值出现在索引 [102][103]（实部、虚部）给出的复数中，那么您的基频将是 (102/2)=61Hz。您必须将其除以二，因为数组的长度是两倍（记住：实部、虚部），因此结果不是索引位置 (102)，而是它的一半 (61)。

如果您打算绘制傅里叶信号，方法也一样。频率 30 由复数 [30][31] 的绝对值给出，依此类推……请记住。由于奈奎斯特冗余（最小采样率必须是信号最高频率的两倍），计算出的 FFT 数组的后半部分必须忽略。它只是前半部分的镜像。如果您想测量高达 6000 的频率，您将需要下一个 2^N 的数字，即 2*6000。

请看我在 CChildView 类的 OnPaint 函数中将 FFT 应用于正弦信号的示例。FFT 在 CFourier 类中实现。这个示例是一个愚蠢的示例，结构也很愚蠢，但我认为它很容易理解。更改参数，玩一下，尝试不同的东西，然后查看结果。这样，您就可以得出自己的结论。

注意事项

如果您打算使用此 FFT 实现，请阅读 NR 许可。

关于这个问题有很多文档。理解起来不容易，但我认为这是一个非常有趣的话题。;)