1

私が間違っていることがわからない。Accelerateフレームワークから得られる結果は私には正しくないようです。どんな助けでも大歓迎です!

これがAForgeとvDPSを比較するいくつかのグラフです AForgeとvDPSを比較したグラフを次に示します。 これは私が実行するvDSPコードです

fftSetup = vDSP_create_fftsetup( 16, 2);

 // Convert the data into a DSPSplitComplex 
int samples = spectrumDataSize;
int samplesOver2 = samples/2;

DSPSplitComplex * complexData = new DSPSplitComplex;
float *realpart = (float *)calloc(samplesOver2, sizeof(float));
float *imagpart = (float *)calloc(samplesOver2, sizeof(float));
complexData->realp = realpart;
complexData->imagp = imagpart;    

vDSP_ctoz((DSPComplex *)realData, 2, complexData, 1,samplesOver2);

// Calculate the FFT
// ( I'm assuming here you've already called vDSP_create_fftsetup() )
vDSP_fft_zrip(fftSetup, complexData, 1, log2f(samples), FFT_FORWARD);

// Scale the data
//float scale = (float) FFT_SCALE; //scale is 32
vDSP_vsmul(complexData->realp, 1, &scale, complexData->realp, 1,samplesOver2);
vDSP_vsmul(complexData->imagp, 1, &scale, complexData->imagp, 1, samplesOver2);


vDSP_zvabs(complexData, 1, spectrumData, 1, samples);

free(complexData->realp);
free(complexData->imagp);
delete complexData;

// All done!
return spectrumData;

これは私がAForgeで行うことです

        foreach (float f in floatData)
            {
                if (i >= this.fft.Length)
                    break;
                this.fft[i++] = new Complex(f * fftSize, 0);
            }
            AForge.Math.FourierTransform.FFT(this.fft, FourierTransform.Direction.Forward);
4

6 に答える 6

2

次のサブルーチンの後

vDSP_ctoz((DSPComplex *)realData, 2, complexData, 1,samplesOver2);

実行され、要素complexDataがあります。samplesOver2でもその後すぐに君は電話する

vDSP_zvabs(complexData, 1, spectrumData, 1, samples);

つまり、2 倍の要素complexDataを持つことが期待されます。samplesこれはできません。

また、どのようにrealData配置されていますか?vDSP_ctoz最初の引数が次の形式で配置されることを期待しているため、私は尋ねます

real0, imag0, real1, imag1, ... real(n-1), imag(n-1).

あなたのデータが実際に本物imag0, imag1, ... imag(n-1)なら、すべて0になるはずです。そうでない場合は、それvDSP_ctozを期待していない可能性があります。(実際のデータを巧妙な方法でパックしていない限り、それは 2 [sic] 半分賢いでしょう!)

最後に、vDSP_create_fftsetup( 16, 2);おそらく次のように変更する必要があります

vDSP_create_fftsetup(16, 0);

================================================== =================

ポストスクリプトに追加された私のサンプルコード:

  FFTSetup fftSetup = vDSP_create_fftsetup(
                                           16,         // vDSP_Length __vDSP_log2n,
                                           kFFTRadix2  // FFTRadix __vDSP_radix
                                           // CAUTION: kFFTRadix2 is an enum that is equal to 0
                                           //          kFFTRadix5 is an enum that is equal to 2
                                           // DO NOT USE 2 IF YOU MEAN kFFTRadix2
                                           );
  NSAssert(fftSetup != NULL, @"vDSP_create_fftsetup() failed to allocate storage");

  int numSamples = 65536;  // numSamples must be an integer power of 2; in this case 65536 = 2 ^ 16
  float realData[numSamples];

  // Prepare the real data with (ahem) fake data, in this case
  // the sum of 3 sinusoidal waves representing a C major chord.
  // The fake data is rigged to have a sampling frequency of 44100 Hz (as for a CD).
  // As always, the Nyquist frequency is just half the sampling frequency, i.e., 22050 Hz.
  for (int i = 0; i < numSamples; i++)
  {
    realData[i] = sin(2 * M_PI * 261.76300048828125 * i / 44100.0)  // C4 = 261.626 Hz
                + sin(2 * M_PI * 329.72717285156250 * i / 44100.0)  // E4 = 329.628 Hz
                + sin(2 * M_PI * 392.30804443359375 * i / 44100.0); // G4 = 391.995 Hz
  }

  float splitReal[numSamples / 2];
  float splitImag[numSamples / 2];

  DSPSplitComplex splitComplex;
  splitComplex.realp = splitReal;
  splitComplex.imagp = splitImag;

  vDSP_ctoz(
            (const DSPComplex *)realData,  // const DSPComplex __vDSP_C[],
            2,                             // vDSP_Stride __vDSP_strideC,  MUST BE A MULTIPLE OF 2
            &splitComplex,                 // DSPSplitComplex *__vDSP_Z,
            1,                             // vDSP_Stride __vDSP_strideZ,
            (numSamples / 2)               // vDSP_Length __vDSP_size
            );

  vDSP_fft_zrip(
                fftSetup,                               // FFTSetup __vDSP_setup,
                &splitComplex,                          // DSPSplitComplex *__vDSP_ioData,
                1,                                      // vDSP_Stride __vDSP_stride,
                (vDSP_Length)lround(log2(numSamples)),  // vDSP_Length __vDSP_log2n,
                // IMPORTANT: THE PRECEDING ARGUMENT MUST BE LOG_BASE_2 OF THE NUMBER OF floats IN splitComplex
                // FOR OUR EXAMPLE, THIS WOULD BE (numSamples / 2) + (numSamples / 2) = numSamples
                kFFTDirection_Forward                   // FFTDirection __vDSP_direction
                );

  printf("DC component = %f\n", splitComplex.realp[0]);
  printf("Nyquist component = %f\n\n", splitComplex.imagp[0]);

  // Next, we compute the Power Spectral Density (PSD) from the FFT.
  // (The PSD is just the magnitude-squared of the FFT.)
  // (We don't bother with scaling as we are only interested in relative values of the PSD.)
  float powerSpectralDensity[(numSamples / 2) + 1];  // the "+ 1" is to make room for the Nyquist component

  // We move the Nyquist component out of splitComplex.imagp[0] and place it
  // at the end of the array powerSpectralDensity, squaring it as we go:
  powerSpectralDensity[numSamples / 2] = splitComplex.imagp[0] * splitComplex.imagp[0];

  // We can now zero out splitComplex.imagp[0] since the imaginary part of the DC component is, in fact, zero:
  splitComplex.imagp[0] = 0.0;

  // Finally, we compute the squares of the magnitudes of the elements of the FFT:
  vDSP_zvmags(
              &splitComplex,         // DSPSplitComplex *__vDSP_A,
              1,                     // vDSP_Stride __vDSP_I,
              powerSpectralDensity,  // float *__vDSP_C,
              1,                     // vDSP_Stride __vDSP_K,
              (numSamples / 2)       // vDSP_Length __vDSP_N
              );

  // We print out a table of the PSD as a function of frequency
  // Replace the "< 600" in the for-loop below with "<= (numSamples / 2)" if you want
  // the entire spectrum up to and including the Nyquist frequency:
  printf("Frequency_in_Hz    Power_Spectral_Density\n");
  for (int i = 0; i < 600; i++)  
  {
    printf("%f,          %f\n", (i / (float)(numSamples / 2)) * 22050.0, powerSpectralDensity[i]);
    // Recall that the array index i = 0 corresponds to zero frequency
    // and that i = (numSamples / 2) corresponds to the Nyquist frequency of 22050 Hz.
    // Frequency values intermediate between these two limits are scaled proportionally (linearly).
  }

  // The output PSD should be zero everywhere except at the three frequencies
  // corresponding to the C major triad.  It should be something like this:

/***************************************************************************
 DC component = -0.000000
 Nyquist component = -0.000000

 Frequency_in_Hz    Power_Spectral_Density
 0.000000,          0.000000
 0.672913,          0.000000
 1.345825,          0.000000
 2.018738,          0.000000
 2.691650,          0.000000
 .
 .
 .
 260.417175,          0.000000
 261.090088,          0.000000
 261.763000,          4294967296.000000
 262.435913,          0.000000
 263.108826,          0.000000
 .
 .
 .
 328.381348,          0.000000
 329.054260,          0.000000
 329.727173,          4294967296.000000
 330.400085,          0.000000
 331.072998,          0.000000
 .
 .
 .
 390.962219,          0.000000
 391.635132,          0.000000
 392.308044,          4294966784.000000
 392.980957,          0.000000
 393.653870,          0.000000
 .
 .
 .
***************************************************************************/

  vDSP_destroy_fftsetup(fftSetup);
于 2012-06-24T02:15:05.437 に答える
2

この行:

vDSP_zvabs(complexData, 1, spectrumData, 1, samples);

次のようにする必要があります。

float cr = complexData->realp[0], ci = complexData->imagp[0];
vDSP_zvabs(complexData, 1, spectrumData, 1, samplesOver2);
spectrumData[0] = cr*cr;
spectrumData[samplesOver2] = ci*ci; // See remarks below.

これは、N サンプルの実数から複素数への FFT が N/2+1 の結果を返すためです。結果のうちの 2 つは実数であり、complexData->realp[0] と complexData->imagp[0] にパックされます。残りの N/2-1 の結果は複素数であり、0 < i < N/2 の場合、complexData->realp[i] の実数成分と complexData->imagp[i] の虚数成分とともに通常格納されます。

vDSP_zvabs は複素数の大きさを計算しますが、最初の出力 (spectrumData[0] 内) は、2 つの数値が [0] 要素にパッキングされているために正しくありません。SpectrumData[0] を cr*cr で上書きすると、それが修正されます。スペースが提供されている場合は、他のパックされた要素 (ナイキスト周波数) の振幅を、spectrumData[samplesOver2] に書き込むこともできます。

その他の注意事項:

SpectrumDataSize は 2 の累乗でなければなりません。

2 を底とする対数を log2f(samples) として計算するのは理想的な方法ではありません。私たち (Apple) は、log2f が 2 の整数乗に対して正確に正しい結果を返すようにしたと思いますが、浮動小数点数の正確さに依存することは、十分に確実になるように注意を払っていない限り避けるべきです。

DSPSplitComplex を「new」で動的に割り当てる必要はありません。これは 2 つのポインターのみを含む POD (plain old data) であるため、単純に「DSPSplitComplex complexData」を宣言して、構造体へのポインターではなく、構造体として使用できます。

于 2012-06-29T15:22:00.223 に答える
0

長さN/2に対して1つのFFTを計算し、長さNに対して1つのFFTを計算しているように見えます。したがって、異なる長さのFFTに対して異なる結果が得られます。

于 2012-06-15T04:03:58.753 に答える
0

vDSP_ctozあなたのデータが偶数奇数ではないことをあなたが呼んでいるので、私は仮定します。その場合は、fft の後に解凍する必要もあります。

vDSPプログラミングガイドから:

実際の FFT を呼び出すアプリケーションは、FFT 呼び出しの前と後の 2 つの変換関数を使用する必要がある場合があります。これは、入力配列が偶奇分割構成でない場合に必要です。

これを説明するサンプルコード

それが役立つことを願っています。

于 2012-06-22T02:27:05.377 に答える
0

いくつかの考え..

int numberOfInputSamples = ..;
int numberOfInputSamplesOver2 = numberOfInputSamples/2;

fftSetup = vDSP_create_fftsetup( log2(numberOfInputSamples), FFT_RADIX2 );
...
Float32 scale = (Float32) 1.0 / (2 * numberOfInputSamples);                
...
float *spectrumData = (float *)calloc( numberOfInputSamplesOver2, sizeof(float));
vDSP_zvabs( complexData, 1, spectrumData, 1, numberOfInputSamplesOver2 );

最後に numberOfInputSamplesOver2 float の大きさになりますよね?

(技術的には numberOfInputSamplesOver2+1 ですが、パッキング全体は別の問題です)

于 2012-06-13T21:19:18.413 に答える
0

私は AForge や Accelerate にはまったく詳しくありませんが、2D 画像を扱う別のプロジェクトで FFT ライブラリをアップグレードするときにいくつかの問題に遭遇しました。

FFT ライブラリからの出力データ表現は一意ではないことが判明しました。一部のアプリケーションでは、出力データを「スワップ」して、低周波数をコーナーではなく中央に配置すると、はるかに便利になります。

FFT アルゴリズムに関するこのページを確認すると、http://www.eso.org/sci/software/eclipse/eug/eug/man/fft.html、両方の形式がサポートされていることがわかります。と記載されています(一番下)。

右にグラフ化したデータは、データ配列の右半分を中央で交換 (ミラーリング) すると、左のデータのように見えるように思えます。

于 2012-06-24T23:57:18.160 に答える