DFT与频率范围

我们需要更改/重新实现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;
    }
}