ДПФ с частотой Диапазон

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

EDIT: Я мог бы не объяснил проблему для ДПФ продуманно, тем не менее, я рад предоставить ответ ниже:

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