ДПФ с частотным диапазоном

Question

Нам нужно изменить / переопределить стандартную реализацию DFT в GSL, которая составляет

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 сэмплов, но нам нужно рассчитать DFT для 128 различных частот.Не могли бы вы помочь мне определить или реализовать требуемое поведение 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;
}