實時做FFT

Question

我想實時地對音頻信號進行FFT，即當人在麥克風中講話時。我將獲取數據（我使用portaudio來做這個，如果wavein會更容易，我會很樂意使用它 - 如果你能告訴我是怎麼做的）。接下來，我使用FFTW庫 - 我知道如何執行一維，二維（實數&複雜）FFT，但我不太確定如何做到這一點，因爲我必須做3D FFT來獲得頻率，幅度（這將確定顏色梯度）和時間。或者它只是一個2D FFT，我得到的幅度和頻率？

Answer 1

如果您需要一幅圖中的幅度，頻率和時間，則該變換被稱爲時頻分解。最流行的稱爲短時傅立葉變換。它的工作原理如下：
1.取一小部分信號（比如說1秒）
2.用小窗口（比如說5毫秒）打開窗口
3.計算加窗信號的一維傅立葉變換。
4.移動窗口少量（2.5毫秒）
5.重複上述步驟直到信號結束。
6.所有這些數據都被輸入到一個矩陣中，然後用它來創建信號的三維表示，顯示其沿頻率，幅度和時間的分解。

窗口的長度將決定您能夠在頻率和時間域中獲得的分辨率。查詢here瞭解更多關於STFT的細節，並尋找「Robi Polikar」關於小波變換的教程，以供外行介紹上述內容。

編輯1：
你把窗口函數（有數不清的功能在那裏 - here is a list最直觀的是一個長方形的窗口，但最常用的是漢明/漢寧窗功能，您可以按照步驟。如果你手中拿着一支紙筆並將它畫出來，

假設你獲得的信號是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.你不明白上述步驟讓我覺得你應該重新訪問時頻分解也是。

Answer 2

您可以通過選擇一個短時間跨度並僅分析（FFT'ing）該時間跨度來創建實時FFT。你可能只需選擇100-500毫秒的非重疊時間段即可脫身;分析更純粹的方法是使用滑動窗口（再次例如100-500毫秒），但這通常是不必要的，並且您可以在沒有太多處理能力的情況下使用不重疊的時間片來顯示漂亮的圖形。

Answer 3

我使用滑動DFT，這是許多倍快速的情況下，你需要做一個傅里葉變換每次樣品到達輸入緩衝區的情況下。

這是基於這樣一個事實，即一旦您對最後N個樣本進行傅里葉變換，並且新樣本到達，您可以「撤消」最老樣本的效果，並應用最新樣本的效果，在單通通過fourier數據！這意味着FFT的滑動DFT性能爲O（n）與O（Log2（n）乘以n）。而且，對於緩衝區大小來保持性能沒有限制。

下面的完整測試程序將滑動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; 
} 

Answer 4

實時FFT指從你剛纔描述的內容完全不同。這意味着對於給定的N和X [N]，您的算法給出Fx [i]，同時遞增值i。意思是，直到當前值計算完成後，計算值纔會計算。這與你剛剛描述的完全不同。

硬件通常使用大約1k-16k點的FFT。固定N，不是實時計算。按照以前的答案所述移動窗口FFT。