C语言 实时进行FFT

Doing FFT in realtime

提问by Rave

I want to do the FFT of an audio signal in real time, meaning while the person is speaking in the microphone. I will fetch the data (I do this with portaudio, if it would be easier with wavein I would be happy to use that - if you can tell me how). Next I am using the FFTW library - I know how to perform 1D, 2D (real&complex) FFT, but I am not so sure how to do this, since I would have to do a 3D FFT to get frequency, amplitude (this would determine the color gradient) and time. Or is it just a 2D FFT, and I get amplitude and frequency?

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

回答by Josh Greifer

I use a Sliding DFT, which is manytimes faster than an FFT in the case where you need to do a fourier transform each time a sample arrives in the input buffer.

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

It's based on the fact that once you have performed a fourier transform for the last N samples, and a new sample arrives, you can "undo" the effect of the oldest sample, and apply the effect of the latest sample, in a single passthrough the fourier data! This means that the sliding DFT performance is O(n)compared with O(Log2(n) times n)for the FFT. Also, there's no restriction to powers of two for the buffer size to maintain performance.

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

The complete test program below compares the sliding DFT with fftw. In my production code I've optimized the below code to unreadibility, to make it three times faster.

下面的完整测试程序将滑动 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;
}

回答by Sriram

If you need amplitude, frequency and time in one graph, then the transform is known as a Time-Frequency decomposition. The most popular one is called the Short Time Fourier Transform. It works as follows:
1. Take a small portion of the signal (say 1 second)
2. Window it with a small window (say 5 ms)
3. Compute the 1D fourier transform of the windowed signal.
4. Move the window by a small amount (2.5 ms)
5. Repeat above steps until end of signal.
6. All of this data is entered into a matrix that is then used to create the kind of 3D representation of the signal that shows its decomposition along frequency, amplitude and time.

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

The length of the window will decide the resolution you are able to obtain in frequency and time domains. Check herefor more details on STFT and search for "Robi Polikar"'s tutorials on wavelet transforms for a layman's introduction to the above.

窗口的长度将决定您能够在频域和时域中获得的分辨率。请点击这里关于STFT更多详细信息，搜索“ROBI Polikar”的上一个外行的介绍上述小波变换的教程。

Edit 1:
You take a windowing function (there are innumerable functions out there - here is a list. Most intuitive is a rectangular window but the most commonly used are the Hamming/Hanning window functions. You can follow the steps below if you have a paper-pencil in hand and draw it along.

编辑 1：
你使用一个窗口函数（那里有无数的函数 -这里是一个列表。最直观的是矩形窗口，但最常用的是汉明/汉宁窗口函数。如果你有一个，你可以按照下面的步骤手里拿着纸笔，沿着它画。

Assume that the signal that you have obtained is 1 sec long and is named x[n]. The windowing function is 5 msec long and is named w[n]. Place the window at the start of the signal (so the end of the window coincides with the 5ms point of the signal) and multiply the x[n]and w[n]like so:
y[n] = x[n] * w[n]- point by point multiplication of the signals.
Take an FFT of y[n].

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

Then you shift the window by a small amount (say 2.5 msec). So now the window stretches from 2.5ms to 7.5 ms of the signal x[n]. Repeat the multiplication and FFT generation steps. In other words, you have an overlap of 2.5 msec. You will see that changing the length of the window and the overlap gives you different resolutions on the time and Frequency axis.

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

Once you do this, you need to feed all the data into a matrix and then have it displayed. The overlap is for minimising the errors that might arise at boundaries and also to get more consistent measurements over such short time frames.

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

P.S: If you had understood STFT and other time-frequency decompositions of a signal, then you should have had no problems with steps 2 and 4. That you have not understood the above mentioned steps makes me feel like you should revisit time-frequency decompositions also.

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

回答by Bernd Elkemann

You can create a realtime FFT by choosing a short time-span and analysing (FFT'ing) just that time-span. You can probably get away with just selecting non-overlapping timespans of say 100-500 milliseconds; the analytically purer way to do this would be using a sliding-window (again of e.g. 100-500 ms), but that is often unnecessary and you can show nice graphics with the non-overlapping timespans without much processing power.

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

回答by sanaris

Real-time FFTmeans completely different from what you just described. It means that for given N and X[N] your algorithm gives Fx[i]while incrementing value i. Meaning, proceeding value does not compute until current value computation completed. This is completely different from what you just described.

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

Hardware usually uses FFT with around 1k-16k points. Fixed N, not real-time computation. Moving window FFT as described with previous answers.

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