频率范围的 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;
}