0

私は電気通信プロジェクトに取り組んでいるコンピューター プログラマーです。
私たちのプロジェクトでは、一連の複素数をフーリエ変換に変更する必要があるFFTため、標準の効率的なコードが必要C89です。
次のコードを使用していますが、十分に機能します。 

    short FFT(short int dir,long m,double *x,double *y)
{
   long n,i,i1,j,k,i2,l,l1,l2;
   double c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   n = 1;
   for (i=0;i<m;i++) 
      n *= 2;

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for (i=0;i<n-1;i++) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* Compute the FFT */
   c1 = -1.0; 
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0; 
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1; 
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1) 
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<n;i++) {
         x[i] /= n;
         y[i] /= n;
      }
   }

   return(true);
}

しかし、このコードは、2^m.like CLRSbook コードのサイズの配列をサポートするだけです。
変換する必要がある配列はこの範囲になく、ゼロを追加するとコストがかかるため、任意のサイズで入力できる別のソリューションを探しています。
何をするかのようなIT++ものmatlab。しかし、純粋なCでそれらを使用したいので、それらを使用することは不可能です。さらに、IT++確認したところ、コードはブロックされていました

4

3 に答える 3

2

大衆市場向けのコンピューティング プラットフォーム (Intel と Windows または OS X、iOS など) で作業している場合は、ベンダーまたはメーカーが提供する高性能 FFT 実装があります。

それ以外の場合は、 FFTWを評価する必要があります。

2 の累乗以外のサイズの高性能 FFT を作成するのは複雑な作業です。

独自の実装を使用する場合は、2 の累乗のサイズに関して:

sqrt表示する実装は、FFT 中に計算します。ほとんどの高性能 FFT 実装では、事前に定数を計算し、それらをテーブルに格納します。

x[i] /= nスケーリングには、およびの除算演算が含まれますy[i] /= n。コンパイラは、これらを除算命令として実装する可能性があります。除算は通常、一般的なプロセッサでは遅い命令です。で割るのではなく、scale = 1. / n一度計算して で掛けたほうがよいでしょう。scalen

さらに良いのは、スケールを完全に省略することです。多くの場合、変換はスケールなしで役立ちます。または、個々の変換からスケールを省略して、集約されたスケールとして 1 回だけ適用することもできます。(たとえば、順方向 FFT と逆方向 FFT の 2 つのスケール操作を実行する代わりに、スケール操作を FFT ルーチンから除外し、順方向 FFT と逆方向 FFT の両方の後に一度だけ実行します。)

周波数領域データをビット反転した順序にすることが許容される場合は、ビット反転順列を省略できる場合があります。

ビット反転順列を維持すると、最適化できます。これを行うための手法は、プラットフォームに依存します。一部のプラットフォームには、整数のビットを反転するための命令があります (たとえば、ARM has rbit)。プラットフォームがそうでない場合は、ビット反転インデックスをテーブルに保持するか、現在のコードよりも高速に計算する方法を調査することをお勧めします。

ビット反転順列とスケーリングの両方を維持する場合は、それらを同時に行うことを検討する必要があります。順列は大量のメモリ モーションを使用し、スケーリングはプロセッサの演算ユニットを使用します。最新のプロセッサのほとんどは、これらの両方を一度に実行できるため、操作を重ねることである程度のメリットが得られます。

現在のコードは基数 2 のバタフライを使用しています。通常は基数 4 の方が優れています。これは、使用するデータの一部を変更したり、加算を減算に変更したり、その逆を変更したりするだけで、i による乗算を実行できるという事実から、いくつかの利点が得られるためです。

配列の長さがプロセッサのレベル 1 メモリ キャッシュのサイズに近づくと、これに対処する適切なコードを設計しない限り (多くの場合、配列の一部をバッファに一時的にコピーすることによって)、FFT 実装の一部がキャッシュをスラッシングし、大幅に速度が低下します。 ）。

ターゲット プロセッサに SIMD 機能がある場合は、FFT でそれらを絶対に使用する必要があります。これらは、FFT パフォーマンスを大幅に高速化します。

上記は、効率的なFFTを書くのは複雑な作業であることを伝えているはずです。多大な労力を費やしたくない場合は、FFTW またはその他の既存の実装を使用する方がよいでしょう。

于 2013-08-20T18:50:16.710 に答える
0

FFTW の代替として、非 2^N FFT 長も処理する私の mix-radix FFT をチェックしてください。C ソースはhttp://www.corix.dk/Mix-FFT/mix-fft.htmlから入手できます。

于 2015-01-10T22:36:49.523 に答える