0

Accelerate vDSP フレームワークを使用して、既存の FFT ベースのローパス フィルターを iOS に移植しようとしています。

サンプルの最初の約 1/4 については、FFT が期待どおりに機能しているようです。しかし、その後、結果は間違っているように見え、さらに奇妙なことがミラーリングされます (信号の後半が前半の大部分をミラーリングします)。

以下のテスト アプリケーションの結果を確認できます。最初に元のサンプリングされたデータがプロットされ、次に予想されるフィルター処理された結果の例 (15Hz を超える信号をフィルター処理)、最後に現在の FFT コードの結果 (目的の結果と FFT の結果の例は異なるスケールであることに注意してください)元のデータ):

FFT結果

ローパス フィルターの実際のコードは次のとおりです。

double *lowpassFilterVector(double *accell, uint32_t sampleCount, double lowPassFreq, double sampleRate )
{
    double stride = 1;

    int ln = log2f(sampleCount);
    int n = 1 << ln;

    // So that we get an FFT of the whole data set, we pad out the array to the next highest power of 2.
    int fullPadN = n * 2;
    double *padAccell = malloc(sizeof(double) * fullPadN);
    memset(padAccell, 0, sizeof(double) * fullPadN);
    memcpy(padAccell, accell, sizeof(double) * sampleCount);

    ln = log2f(fullPadN);
    n = 1 << ln;

    int nOver2 = n/2;

    DSPDoubleSplitComplex A;
    A.realp = (double *)malloc(sizeof(double) * nOver2);
    A.imagp = (double *)malloc(sizeof(double) * nOver2);

    // This can be reused, just including it here for simplicity.
    FFTSetupD setupReal = vDSP_create_fftsetupD(ln, FFT_RADIX2);

    vDSP_ctozD((DSPDoubleComplex*)padAccell,2,&A,1,nOver2);

    // Use the FFT to get frequency counts
    vDSP_fft_zripD(setupReal, &A, stride, ln, FFT_FORWARD);


    const double factor = 0.5f;
    vDSP_vsmulD(A.realp, 1, &factor, A.realp, 1, nOver2);
    vDSP_vsmulD(A.imagp, 1, &factor, A.imagp, 1, nOver2);
    A.realp[nOver2] = A.imagp[0];
    A.imagp[0] = 0.0f;
    A.imagp[nOver2] = 0.0f;

    // Set frequencies above target to 0.

    // This tells us which bin the frequencies over the minimum desired correspond to
    NSInteger binLocation = (lowPassFreq * n) / sampleRate;

    // We add 2 because bin 0 holds special FFT meta data, so bins really start at "1" - and we want to filter out anything OVER the target frequency
    for ( NSInteger i = binLocation+2; i < nOver2; i++ )
    {
        A.realp[i] = 0;
    }

    // Clear out all imaginary parts
    bzero(A.imagp, (nOver2) * sizeof(double));
    //A.imagp[0] = A.realp[nOver2];


    // Now shift back all of the values
    vDSP_fft_zripD(setupReal, &A, stride, ln, FFT_INVERSE);

    double *filteredAccell = (double *)malloc(sizeof(double) * fullPadN);

    // Converts complex vector back into 2D array
    vDSP_ztocD(&A, stride, (DSPDoubleComplex*)filteredAccell, 2, nOver2);

    //  Have to scale results to account for Apple's FFT library algorithm, see:
    // http://developer.apple.com/library/ios/#documentation/Performance/Conceptual/vDSP_Programming_Guide/UsingFourierTransforms/UsingFourierTransforms.html#//apple_ref/doc/uid/TP40005147-CH202-15952
    double scale = (float)1.0f / fullPadN;//(2.0f * (float)n);
    vDSP_vsmulD(filteredAccell, 1, &scale, filteredAccell, 1, fullPadN);

    // Tracks results of conversion
    printf("\nInput & output:\n");
    for (int k = 0; k < sampleCount; k++)
    {
        printf("%3d\t%6.2f\t%6.2f\t%6.2f\n", k, accell[k], padAccell[k], filteredAccell[k]);
    }


    // Acceleration data will be replaced in-place.
    return filteredAccell;
}

元のコードでは、ライブラリは 2 のべき乗でないサイズの入力データを処理していました。私の Accelerate コードでは、最も近い 2 のべき乗まで入力をパディングしています。以下のサンプル テストの場合、元のサンプル データは 1000 サンプルなので、1024 にパディングされます。これが結果に影響を与えるとは思いませんが、可能な違いのために含めます。

ソリューションを試してみたい場合は、グラフを生成するサンプル プロジェクトをここ (FFTTest フォルダー内) からダウンロードできます。

FFT サンプル プロジェクト コード

洞察に感謝します。私は以前にFFTを使用したことがないため、重要な何かが欠けているように感じます。

4

1 に答える 1

3

厳密に実数 (複素数ではない) の結果が必要な場合は、IFFT の前のデータが共役対称でなければなりません。結果をミラー対称にしたくない場合は、IFFT の前に虚数成分をゼロにしないでください。IFFT が通過帯域に大量のリップルを持つフィルターを作成する前にビンをゼロにするだけです。

Accelerate フレームワークは、2 のべき乗だけでなく、より多くの FFT 長もサポートしています。

于 2013-05-28T05:12:38.720 に答える