DFT

Question

我々はサンプル/入力サイズ2回入力ベクトルオーバー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における標準DFTの実装を、再実装/変更する必要があります。しかし、異なる周波数ビンを構築する必要があります。たとえば、4096サンプルありますが、128の異なる周波数に対してDFTを計算する必要があります。必要なDFTの振る舞いを定義または実装する手助けをしてください。前もって感謝します。

EDIT：最初に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

：私は精巧なDFTのための問題を説明していないかもしれないが、それにもかかわらず、私は以下の答えを提供するために満足しています

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