实时进行FFT
我想实时地对音频信号进行FFT，这意味着当人在麦克风中讲话时。 我将获取数据（我使用portaudio执行此操作，如果使用wavein会更容易，我会乐意使用它 – 如果您可以告诉我如何）。 接下来我使用的是FFTW库 – 我知道如何执行1D，2D（实数和复数）FFT，但我不太清楚如何做到这一点，因为我必须进行3D FFT以获得频率，幅度（这将决定颜色渐变）和时间。 或者它只是一个2D FFT，我得到幅度和频率？
如果在一个图中需要幅度，频率和时间，则变换称为时频分解。 最受欢迎的称为短时傅里叶变换。 它的工作原理如下：
1.取一小部分信号（比如1秒）
2.用一个小窗口（比如5毫秒）窗口
3.计算窗口信号的1D傅立叶变换。
4.将窗口移动少量（2.5毫秒）
5.重复上述步骤，直到信号结束。
6.所有这些数据都输入矩阵，然后用于创建信号的3D表示，显示其沿频率，幅度和时间的分解。
窗口的长度将决定您在频域和时域中可以获得的分辨率。 在这里查看有关STFT的更多详细信息，并搜索“Robi Polikar”关于小波变换的教程，以便外行人介绍上述内容。
编辑1：
你有一个窗口function（有无数的function – 这里是一个列表 。最直观的是一个矩形窗口，但最常用的是汉明/汉宁窗口function。如果你有一支纸笔，你可以按照下面的步骤在手，并绘制它。
假设您获得的信号长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.5ms延伸到7.5ms。 重复乘法和FFT生成步骤。 换句话说，您有2.5毫秒的重叠。 您将看到更改窗口长度和重叠会在时间和频率轴上提供不同的分辨率。
完成此操作后，您需要将所有数据输入矩阵然后显示。 重叠用于最小化在边界处可能出现的误差，并且还在这样的短时间帧内获得更一致的测量。
PS：如果你已经理解了信号的STFT和其他时频分解，那么你应该对第2步和第4步没有任何问题。你还没有理解上面提到的步骤让我觉得你应该重新考虑时频分解也。
我使用Sliding DFT，在每次样本到达输入缓冲区时需要进行傅里叶变换的情况下，它比FFT快许多倍。
这是基于以下事实：一旦您对最后N个样本执行了傅里叶变换，并且新样本到达，您可以“撤消”最旧样本的效果，并在一次通过中应用最新样本的效果通过傅里叶数据！ 这意味着对于FFT，滑动DFT性能是O（n）与O（Log2（n）乘以n）相比。 此外，对于保持性能的缓冲区大小，对2的幂没有限制。
下面的完整测试程序将滑动DFT与fftw进行比较。 在我的生产代码中，我将以下代码优化为不可读性，使其速度提高了三倍。
#include #include #include #include #include #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::fftw_plan plan_fwd; fftw::fftw_plan plan_inv; #endif typedef std::complex 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）仅时间跨度来创建实时FFT。 你可以选择不重叠的时间跨度，例如100-500毫秒; 分析上更纯粹的方法是使用滑动窗口（例如100-500毫秒），但这通常是不必要的，你可以用非重叠的时间跨度显示漂亮的图形而没有太多的处理能力。
实时FFT意味着与您刚才描述的完全不同。 这意味着对于给定的N和X [N]，您的算法在递增值i时给出Fx [i] 。 意思是，在当前值计算完成之前，不会计算进行值。 这与您刚才描述的完全不同。
硬件通常使用FFT，大约1k-16k点。 修正了N，而不是实时计算。 移动窗口FFT，如前面的答案所述。