1

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 はサンプル/入力サイズに対して入力ベクトルを 2 回繰り返します。ただし、異なる周波数ビンを構築する必要があります。たとえば、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;
    }
}
4

1 に答える 1

0

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;
}
于 2012-01-16T18:40:09.597 に答える