DFT與頻率範圍

Question

我們需要改變/重新實現在GSL標準DFT實現，這是

int 
FUNCTION(gsl_dft_complex,transform) (const BASE data[], 
            const size_t stride, const size_t n, 
            BASE result[], 
            const gsl_fft_direction sign) 
{ 

    size_t i, j, exponent; 

    const double d_theta = 2.0 * ((int) sign) * M_PI/(double) n; 

    /* FIXME: check that input length == output length and give error */ 

    for (i = 0; i < n; i++) 
    { 
     ATOMIC sum_real = 0; 
     ATOMIC sum_imag = 0; 

     exponent = 0; 

     for (j = 0; j < n; j++) 
     { 
      double theta = d_theta * (double) exponent; 
      /* sum = exp(i theta) * data[j] */ 

      ATOMIC w_real = (ATOMIC) cos (theta); 
      ATOMIC w_imag = (ATOMIC) sin (theta); 

      ATOMIC data_real = REAL(data,stride,j); 
      ATOMIC data_imag = IMAG(data,stride,j); 

      sum_real += w_real * data_real - w_imag * data_imag; 
      sum_imag += w_real * data_imag + w_imag * data_real; 

      exponent = (exponent + i) % n; 
     } 
     REAL(result,stride,i) = sum_real; 
     IMAG(result,stride,i) = sum_imag; 
    } 

    return 0; 
} 

在該實現中，GSL迭代輸入矢量兩次採樣/輸入大小。但是，我們需要構建不同的頻率分檔。例如，我們有4096個樣本，但是我們需要計算128個不同頻率的DFT。你能幫我定義或實施必要的DFT行爲嗎？提前致謝。

編輯：我們不搜索第一個m頻率。

其實，以下方法正確找到DFT結果與給定的頻率斌號碼？ N =樣本大小 B =頻率區間大小

k = 0,...,127 X[k] = SUM(0,N){x[i]*exp(-j*2*pi*k*i/B)} 

編輯：我可能沒有精心解釋了DFT的問題，不過，我很樂意爲您提供以下答案：

void compute_dft(const std::vector<double>& signal, 
       const std::vector<double>& frequency_band, 
       std::vector<double>& result, 
       const double sampling_rate) 
{ 
    if(0 == result.size() || result.size() != (frequency_band.size() << 1)){ 
     result.resize(frequency_band.size() << 1, 0.0); 
    } 

    //note complex signal assumption 
    const double d_theta = -2.0 * PI * sampling_rate; 

    for(size_t k = 0; k < frequency_band.size(); ++k){ 
     const double f_k = frequency_band[k]; 
     double real_sum = 0.0; 
     double imag_sum = 0.0; 

     for(size_t n = 0; n < (signal.size() >> 1); ++n){ 
      double theta = d_theta * f_k * (n + 1); 

      double w_real = cos(theta); 
      double w_imag = sin(theta); 

      double d_real = signal[2*n]; 
      double d_imag = signal[2*n + 1]; 

      real_sum += w_real * d_real - w_imag * d_imag; 
      imag_sum += w_real * d_imag + w_imag * d_real; 
     } 

     result[2*k] = real_sum; 
     result[2*k + 1] = imag_sum; 
    } 
} 

Answer 1

假設你只是想要第一個m輸出頻率：

int 
FUNCTION(gsl_dft_complex,transform) (const BASE data[], 
            const size_t stride, 
            const size_t n, // input size 
            const size_t m, // output size (m <= n) 
            BASE result[], 
            const gsl_fft_direction sign) 
{ 

    size_t i, j, exponent; 

    const double d_theta = 2.0 * ((int) sign) * M_PI/(double) n; 

    /* FIXME: check that m <= n and give error */ 

    for (i = 0; i < m; i++) // for each of m output bins 
    { 
     ATOMIC sum_real = 0; 
     ATOMIC sum_imag = 0; 

     exponent = 0; 

     for (j = 0; j < n; j++) // for each of n input points 
     { 
      double theta = d_theta * (double) exponent; 
      /* sum = exp(i theta) * data[j] */ 

      ATOMIC w_real = (ATOMIC) cos (theta); 
      ATOMIC w_imag = (ATOMIC) sin (theta); 

      ATOMIC data_real = REAL(data,stride,j); 
      ATOMIC data_imag = IMAG(data,stride,j); 

      sum_real += w_real * data_real - w_imag * data_imag; 
      sum_imag += w_real * data_imag + w_imag * data_real; 

      exponent = (exponent + i) % n; 
     } 
     REAL(result,stride,i) = sum_real; 
     IMAG(result,stride,i) = sum_imag; 
    } 

    return 0; 
}