实时进行 FFT答案

【问题标题】：Doing FFT in realtime实时进行 FFT
【发布时间】：2011-10-03 13:00:40
【问题描述】：

我想实时对音频信号进行 FFT，这意味着当这个人在麦克风里说话时。我将获取数据（我使用 portaudio 进行此操作，如果使用 wavein 会更容易，我很乐意使用它 - 如果您能告诉我如何使用）。接下来我使用 FFTW 库 - 我知道如何执行 1D、2D（实数和复数）FFT，但我不太确定如何执行此操作，因为我必须执行 3D FFT 才能获得频率、幅度（这将确定颜色渐变）和时间。或者它只是一个 2D FFT，我得到幅度和频率？

【问题讨论】：

Audio spectrum analysis using FFT algorithm in Java的可能重复
@rave：一直问同样的问题是不礼貌的stackoverflow.com/questions/6635400/…
不是拼贴作业，大学毕业2个月了。是的，它在 matlab 中要简单得多，但我想用 C 来做。@msw 对不起..我很绝望......这两天有问题。我习惯在 Redhat 中编程，但是当我切换到 Windows 时，我不得不经历学习使用 MSVC++ 及其所有链接的麻烦。
@Rave：如果你想用 C 语言来做，到目前为止你做了什么？给我们看看代码？
@Rave：没有冒犯的意思：只是 2 个 FFT 问题一起弹出！ :)

标签： c real-time signal-processing fft

【解决方案1】：

如果您需要一张图中的幅度、频率和时间，那么这种变换称为时频分解。最流行的一种称为短时傅里叶变换。它的工作原理如下：
1. 取一小部分信号（比如 1 秒）
2. 用一个小窗口（比如 5 毫秒）打开它
3. 计算加窗信号的一维傅里叶变换。
4. 将窗口移动少量（2.5 毫秒）
5. 重复上述步骤直到信号结束。
6. 所有这些数据都输入到矩阵中，然后用于创建信号的 3D 表示，显示信号沿频率、幅度和时间的分解。

窗口的长度将决定您能够在频域和时域中获得的分辨率。查看here 了解有关 STFT 的更多详细信息，并搜索“Robi Polikar”关于小波变换的教程，以获得外行对上述内容的介绍。

编辑 1：
您使用一个窗口函数（那里有无数函数 - here is a list。最直观的是矩形窗口，但最常用的是 Hamming/Hanning 窗口函数。如果您有纸笔，可以按照以下步骤操作用手画出来。

假设您获得的信号长度为 1 秒，名称为x[n]。窗口函数的长度为 5 毫秒，命名为 w[n]。将窗口放在信号的开头（因此窗口的结尾与信号的 5ms 点重合）并将x[n] 和 w[n] 相乘，如下所示：
y[n] = x[n] * w[n] - 逐点相乘信号。
对y[n] 进行 FFT。

然后您将窗口移动少量（例如 2.5 毫秒）。所以现在窗口从信号x[n] 的 2.5 毫秒延伸到 7.5 毫秒。重复乘法和 FFT 生成步骤。换句话说，您有 2.5 毫秒的重叠。您将看到更改窗口的长度和重叠会在时间轴和频率轴上为您提供不同的分辨率。

完成此操作后，您需要将所有数据输入一个矩阵，然后将其显示出来。重叠是为了最大限度地减少边界处可能出现的错误，并在如此短的时间范围内获得更一致的测量结果。

PS：如果您了解 STFT 和信号的其他时频分解，那么您应该对第 2 步和第 4 步没有任何问题。您没有理解上述步骤让我觉得您应该重新审视时间-频率分解也。

【讨论】：

嘿，我了解 STFT 背后的理论，但我有点困惑你将如何在 C 中做到这一点，你的意思是用一个小窗口窗口它是什么意思，以及第 4 步，移动窗口？第2步和第4步没看懂，能解释一下吗，谢谢
@Rave：我已经解释了一点。看看吧。
嗨，Sriram，谢谢，我找到了这个库，名为 siglib:numerix-dsp.com/free。它具有汉明和矩形等窗口功能。如果有人对 siglib 有任何经验，我很想听听
恐怕我对 siglib 没有任何经验。但是非常欢迎您尝试一些代码并提出问题。此外，如果您发现此答案有帮助，您可能希望将其标记为正确答案。
@Sriram - 是否可以不只根据每个音频回调的缓冲区元素计算窗口？

【解决方案2】：

实时 FFT 的意思与你刚才描述的完全不同。这意味着对于给定的 N 和 X[N]，您的算法在增加值 i 的同时给出 Fx[i]。意思是，在当前值计算完成之前，不会计算后续值。这和你刚才描述的完全不同。

硬件通常使用大约 1k-16k 点的 FFT。固定 N，不是实时计算。如先前答案所述移动窗口 FFT。

【讨论】：

【解决方案3】：

我使用滑动 DFT，在每次样本到达输入缓冲区时都需要进行傅立叶变换的情况下，它比 FFT 快许多倍。

这是基于这样一个事实，即一旦您对最后 N 个样本执行了傅立叶变换，并且新样本到达，您可以“撤消”最旧样本的效果，并应用最新样本的效果，在单次传递傅立叶数据！这意味着滑动 DFT 的性能为 O(n)，而 FFT 的性能为 O(Log2(n) 乘以 n)。此外，为了保持性能，缓冲区大小没有 2 次方的限制。

下面的完整测试程序比较了滑动DFT和fftw。在我的生产代码中，我将以下代码优化为不可读，使其速度提高了三倍。

    #include <complex>
#include <iostream>
#include <time.h>
#include <math_defines.h>
#include <float.h>

#define DO_FFTW // libfftw
#define DO_SDFT

#if defined(DO_FFTW)
#pragma comment( lib, "d:\\projects\\common\\fftw\\libfftw3-3.lib" )

namespace fftw {
#include <fftw/fftw3.h>
}

fftw::fftw_plan plan_fwd;
fftw::fftw_plan plan_inv;

#endif

typedef std::complex<double> complex;

// Buffer size, make it a power of two if you want to improve fftw
const int N = 750;


// input signal
complex in[N];

// frequencies of input signal after ft
// Size increased by one because the optimized sdft code writes data to freqs[N]
complex freqs[N+1];

// output signal after inverse ft of freqs
complex out1[N];
complex out2[N];

// forward coeffs -2 PI e^iw -- normalized (divided by N)
complex coeffs[N];
// inverse coeffs 2 PI e^iw
complex icoeffs[N];

// global index for input and output signals
int idx;


// these are just there to optimize (get rid of index lookups in sdft)
complex oldest_data, newest_data;

//initilaize e-to-the-i-thetas for theta = 0..2PI in increments of 1/N
void init_coeffs()
{
    for (int i = 0; i < N; ++i) {
        double a = -2.0 * PI * i  / double(N);
        coeffs[i] = complex(cos(a)/* / N */, sin(a) /* / N */);
    }
    for (int i = 0; i < N; ++i) {
        double a = 2.0 * PI * i  / double(N);
        icoeffs[i] = complex(cos(a),sin(a));
    }
}


// initialize all data buffers
void init()
{
    // clear data
    for (int i = 0; i < N; ++i)
        in[i] = 0;
    // seed rand()
    srand(857);
    init_coeffs();
    oldest_data = newest_data = 0.0;
    idx = 0;
}

// simulating adding data to circular buffer
void add_data()
{
    oldest_data = in[idx];
    newest_data = in[idx] = complex(rand() / double(N));

}


// sliding dft
void sdft()
{
    complex delta = newest_data - oldest_data;
    int ci = 0;
    for (int i = 0; i < N; ++i) {
        freqs[i] += delta * coeffs[ci];
        if ((ci += idx) >= N)
            ci -= N;
    }
}

// sliding inverse dft
void isdft()
{
    complex delta = newest_data - oldest_data;
    int ci = 0;
    for (int i = 0; i < N; ++i) {
        freqs[i] += delta * icoeffs[ci];
        if ((ci += idx) >= N)
            ci -= N;
    }
}

// "textbook" slow dft, nested loops, O(N*N)
void ft()
{
    for (int i = 0; i < N; ++i) {
        freqs[i] = 0.0;
        for (int j = 0; j < N; ++j) {
            double a = -2.0 * PI * i * j / double(N);
            freqs[i] += in[j] * complex(cos(a),sin(a));
        }
    }
}

double mag(complex& c)
{
    return sqrt(c.real() * c.real() + c.imag() * c.imag());
}

void powr_spectrum(double *powr)
{
    for (int i = 0; i < N/2; ++i) {
        powr[i] = mag(freqs[i]);
    }

}

int main(int argc, char *argv[])
{
    const int NSAMPS = N*10;
    clock_t start, finish;

#if defined(DO_SDFT)
// ------------------------------ SDFT ---------------------------------------------
    init();

    start = clock();
    for (int i = 0; i < NSAMPS; ++i) {

        add_data();

        sdft();
        // Mess about with freqs[] here
        //isdft();

        if (++idx == N) idx = 0; // bump global index

        if ((i % 1000) == 0)
            std::cerr << i << " iters..." << '\r';
    }
    finish = clock();

    std::cout << "SDFT: " << NSAMPS / ((finish-start) / (double)CLOCKS_PER_SEC) << " fts per second." << std::endl;

    double powr1[N/2];
    powr_spectrum(powr1);
#endif

#if defined(DO_FFTW)

// ------------------------------ FFTW ---------------------------------------------
    plan_fwd = fftw::fftw_plan_dft_1d(N, (fftw::fftw_complex *)in, (fftw::fftw_complex *)freqs, FFTW_FORWARD, FFTW_MEASURE);
    plan_inv = fftw::fftw_plan_dft_1d(N, (fftw::fftw_complex *)freqs, (fftw::fftw_complex *)out2, FFTW_BACKWARD, FFTW_MEASURE);

    init();

    start = clock();
    for (int i = 0; i < NSAMPS; ++i) {

        add_data();

        fftw::fftw_execute(plan_fwd);
        // mess about with freqs here
        //fftw::fftw_execute(plan_inv);

        if (++idx == N) idx = 0; // bump global index

        if ((i % 1000) == 0)
            std::cerr << i << " iters..." << '\r';
    }
    // normalize fftw's output 
    for (int j = 0; j < N; ++j)
        out2[j] /= N;

    finish = clock();

    std::cout << "FFTW: " << NSAMPS / ((finish-start) / (double)CLOCKS_PER_SEC) << " fts per second." << std::endl;
    fftw::fftw_destroy_plan(plan_fwd);
    fftw::fftw_destroy_plan(plan_inv);

    double powr2[N/2];
    powr_spectrum(powr2);
#endif
#if defined(DO_SDFT) && defined(DO_FFTW)
// ------------------------------      ---------------------------------------------
    const double MAX_PERMISSIBLE_DIFF = 1e-11; // DBL_EPSILON;
    double diff;
    // check my ft gives same power spectrum as FFTW
    for (int i = 0; i < N/2; ++i)
        if ( (diff = abs(powr1[i] - powr2[i])) > MAX_PERMISSIBLE_DIFF)
            printf("Values differ by more than %g at index %d.  Diff = %g\n", MAX_PERMISSIBLE_DIFF, i, diff);

#endif
    return 0;
}

【讨论】：

嗨乔希，这个优化看起来非常有前途！我的场景涉及在每次 FFT 之前通过汉宁窗口对输入数据进行加权（即我没有像您的代码那样使用矩形窗口），您是否碰巧知道您的算法是否有一个变体可以与非- 矩形窗口数据？非常感谢您！ =)
没有这样做。但是Hamming(idx) = 0.54 - 0.46 * cos(2PI * idx / (N-1))。因此，与 FT 本身一样，它可以在 idx 扫描缓冲区时逐步运行。要更改的行是complex delta = newest_data - oldest_data; - 不是取最旧和最新数据值之间的差异，而是需要在取差异之前对这两个值应用汉明函数。当我有时间时，我会尝试修改我的代码来做到这一点。 编辑： 你说的是汉宁——你是说汉宁还是汉明？汉恩是Hann(n) = 0.5 - 0.5 * cos(2PI * n / (N-1))
Sebastian，我已经有一段时间没有看过这段代码了，我从第一原理“手工”推导出来的。我当时不知道我重新发明了 SDFT，直到我实现了它。其实我有点失望！我对奇怪的索引感到有些困惑：我认为这是我使用循环缓冲区、碰撞全局指针和避免代码中的任何模运算的副产品。所以我猜数学是教科书，索引不是。
@Sebastian：我也想知道同样的事情。（即，为什么不将freqs[i] 相乘，如果你看一下数学就应该如此）正如乔希已经说过的，这是因为循环缓冲区，但我会澄清一下：代码不是' t 计算从窗口直接看到的信号的 DFT，但也旋转了i % N 步长。 IE。对于[1, 2, 3, 4, 5, 6, 7, 8, 9] 和N=5，当i=2 计算DFT 时为[6, 7, 3, 4, 5] 而不是[3, 4, 5, 6, 7]。当像这样旋转时，除一个之外的所有数字都保持不变，freqs[i] 本身不会相乘。
@KenVernaillen 我会挖掘它并在找到它时发布一个链接

【解决方案4】：

您可以通过选择较短的时间跨度并分析（FFT'ing）该时间跨度来创建实时 FFT。您可能只需选择 100-500 毫秒的非重叠时间跨度即可；更纯粹的分析方法是使用滑动窗口（同样为 100-500 毫秒），但这通常是不必要的，您可以在没有太多处理能力的情况下显示具有非重叠时间跨度的漂亮图形。

【讨论】：