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; 
}