6

私はおそらくこれを間違って尋ねて、自分自身を非常に愚かに見えるようにするつもりですが、ここに行きます：

.wavファイルでオーディオの操作と処理を行おうとしています。これで、すべてのデータ（ヘッダーを含む）を読み取ることができますが、データの頻度が必要です。これを行うには、FFTを使用する必要があります。

インターネットを上下に検索して見つけたところ、「Cの数値レシピ」の本から例を取り出しましたが、配列の代わりにベクトルを使用するように修正しました。さて、ここに問題があります：

私は（使用例として）一連の数値とサンプリングレートを与えられました：

X = {50, 206, -100, -65, -50, -6, 100, -135}

サンプリングレート：8000サンプル数：8

したがって、これに答える必要があります。

  0Hz     A=0       D=1.57079633
  1000Hz     A=50      D=1.57079633
  2000HZ     A=100     D=0
  3000HZ     A=100     D=0
  4000HZ     A=0       D=3.14159265

私が書き直したコードはコンパイルされますが、これらの数値を方程式（関数）に入力しようとすると、セグメンテーション違反が発生します。コードに問題がありますか、またはサンプリングレートが高すぎますか？（はるかに小さいサンプリングレートを使用する場合、アルゴリズムはセグメント化しません）。コードは次のとおりです。

#include <iostream>
#include <math.h>
#include <vector>
using namespace std;

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr;
#define pi 3.14159

void ComplexFFT(vector<float> &realData, vector<float> &actualData, unsigned long sample_num, unsigned int sample_rate, int sign)
{
    unsigned long n, mmax, m, j, istep, i;
    double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;

    // CHECK TO SEE IF VECTOR IS EMPTY;

    actualData.resize(2*sample_rate, 0);

    for(n=0; (n < sample_rate); n++)
    {
        if(n < sample_num)
        {
            actualData[2*n] = realData[n];
        }else{
            actualData[2*n] = 0;
            actualData[2*n+1] = 0;
        }
    }

    // Binary Inversion
    n = sample_rate << 1;
    j = 0;

    for(i=0; (i< n /2); i+=2)
    {
        if(j > i)
        {
            SWAP(actualData[j], actualData[i]);
            SWAP(actualData[j+1], actualData[i+1]);
            if((j/2)<(n/4))
            {
                SWAP(actualData[(n-(i+2))], actualData[(n-(j+2))]);
                SWAP(actualData[(n-(i+2))+1], actualData[(n-(j+2))+1]);
            }
        }
        m = n >> 1;
         while (m >= 2 && j >= m) {
          j -= m;
          m >>= 1;
         }
         j += m;
     }
     mmax=2;

     while(n > mmax) {

        istep = mmax << 1;
        theta = sign * (2*pi/mmax);
        wtemp = sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi = sin(theta);
        wr = 1.0;
        wi = 0.0;

        for(m=1; (m < mmax); m+=2) {
            for(i=m; (i <= n); i += istep)
            {
                j = i*mmax;
                tempr = wr*actualData[j-1]-wi*actualData[j];
                tempi = wr*actualData[j]+wi*actualData[j-1];

                actualData[j-1] = actualData[i-1] - tempr;
                actualData[j] = actualData[i]-tempi;
                actualData[i-1] += tempr;
                actualData[i] += tempi;
            }
            wr = (wtemp=wr)*wpr-wi*wpi+wr;
            wi = wi*wpr+wtemp*wpi+wi;
        }
        mmax = istep;
    }

    // determine if the fundamental frequency
    int fundemental_frequency = 0;
    for(i=2; (i <= sample_rate); i+=2)
    {
        if((pow(actualData[i], 2)+pow(actualData[i+1], 2)) > pow(actualData[fundemental_frequency], 2)+pow(actualData[fundemental_frequency+1], 2)) {
            fundemental_frequency = i;
        }

    }
}
int main(int argc, char *argv[]) {

    vector<float> numbers;
    vector<float> realNumbers;

    numbers.push_back(50);
    numbers.push_back(206);
    numbers.push_back(-100);
    numbers.push_back(-65);
    numbers.push_back(-50);
    numbers.push_back(-6);
    numbers.push_back(100);
    numbers.push_back(-135);

    ComplexFFT(numbers, realNumbers, 8, 8000, 0);

    for(int i=0; (i < realNumbers.size()); i++)
    {
        cout << realNumbers[i] << "\n";
    }
}

もう1つは（これはばかげているように聞こえますが）、ComplexFFT関数を介して渡される「intsign」に何が期待されるのかはよくわかりません。ここで問題が発生する可能性があります。

誰かがこの問題に対する提案や解決策を持っていますか？

ありがとうございました ：）

4

3 に答える 3

4

問題は、アルゴリズムの翻訳方法の誤りにあると思います。

  • ではなくに初期化jするつもりでしたか?10

  • for(i = 0; (i < n/2); i += 2)おそらくあるはずfor (i = 1; i < n; i += 2)です。

  • あなたSWAPの s はおそらく

    SWAP(actualData[j - 1], actualData[i - 1]);
SWAP(actualData[j], actualData[i]);

  • 次SWAPの s は何のためですか? それらは必要ないと思います。

    if((j/2)<(n/4))
{
    SWAP(actualData[(n-(i+2))], actualData[(n-(j+2))]);
    SWAP(actualData[(n-(i+2))+1], actualData[(n-(j+2))+1]);
}

  • ビット反転を行う場合は、おそらくj >= min にするwhile (m >= 2 && j >= m)必要があります。j > m

  • Danielson-Lanczos セクションを実装するコードではj = i*mmax;、追加する必要はありませんでしたj = i + mmax;か?

それとは別に、コードを単純化するためにできることはたくさんあります。

SWAPマクロを使用できる場合は、マクロの使用をお勧めしませんstd::swap... を提案するつもりでしたstd::swap_rangesが、データはすべて実数であるため、実数部分を交換するだけでよいことに気付きました (時系列の虚数部分はすべて0):

std::swap(actualData[j - 1], actualData[i - 1]);

を使用して全体を単純化することもできstd::complexます。

于 2012-08-07T18:30:43.193 に答える
2

Cの数値レシピのFFTはCooley-Tukeyアルゴリズムを使用しているため、最後の質問に答えると、int sign渡されることで同じルーチンを使用して順方向（sign=-1）と逆方向（）の両方のsign=1FFTを計算できます。signこれは、を定義するときに使用している方法と一致しているようですtheta = sign * (2*pi/mmax)。

于 2012-08-07T18:29:28.590 に答える
2

ベクターのサイズ変更が原因だと思います。

1つの可能性:サイズを変更すると、スタックに一時オブジェクトが作成されてから、ヒープに戻ると思います。

于 2012-08-07T17:57:54.510 に答える