radix-4 FFT-Implementierung

Ich habe eine 4-Punkt-Radix-4-FFT implementiert und festgestellt, dass ich einige Manipulationen an den Ausgangstermen vornehmen muss, damit sie mit einer dft übereinstimmen.

Mein Code ist eine ziemlich direkte Implementierung der Matrixformulierung, daher ist mir nicht klar, wo das Problem liegt

//                                | 
// radix-4 butterfly matrix form  |  complex multiplication
//                                | 
//        +-          -+ +-  -+   |    a+ib
// X[0] = | 1  1  1  1 | |x[0]|   |  * c+id
// X[1] = | 1 -i -1  i | |x[1]|   |    -------
// X[2] = | 1 -1  1 -1 | |x[2]|   |    ac + ibc
// X[3] = | 1  i -1 -i | |x[3]|   |         iad - bd
//        +-          -+ +-  -+   |    ------------------
//                                |    (ac-bd) + i(bc+ad)  
//                                | 

Kann jemand erkennen, wo ich falsch gelaufen bin?

Vielen Dank,

-David

typedef double fp; // base floating-point type


// naiive N-point DFT implementation as reference to check fft implementation against
//
void dft(int inv, struct cfp *x, struct cfp *y, int N) {

  long int i, j;
  struct cfp w;
  fp ang;

  for(i=0; i<N; i++) { // do N-point FFT/IFFT
    y[i].r = y[i].i = 0;
    if (inv) ang =  2*PI*(fp)i/(fp)N;
    else     ang = -2*PI*(fp)i/(fp)N;
    for (j=0; j<N; j++) {
      w.r = cos(j*ang);
      w.i = sin(j*ang);
      y[i].r += (x[j].r * w.r - x[j].i * w.i);
      y[i].i += (x[j].r * w.i + x[j].i * w.r);
    }
  }

  // scale output in the case of an IFFT
  if (inv) {  
    for (i=0; i<N; i++) {
      y[i].r = y[i].r/(fp)N;
      y[i].i = y[i].i/(fp)N;
    }
  }

} // dft()


void r4fft4(int inv, int reorder, struct cfp *x, struct cfp *y) {
  struct cfp x1[4], w[4];
  fp         ang, temp;
  int        i;

  //                                | 
  // radix-4 butterfly matrix form  |  complex multiplication
  //                                | 
  //        +-          -+ +-  -+   |    a+ib
  // y[0] = | 1  1  1  1 | |x[0]|   |  * c+id
  // y[1] = | 1 -i -1  i | |x[1]|   |    -------
  // y[2] = | 1 -1  1 -1 | |x[2]|   |    ac + ibc
  // y[3] = | 1  i -1 -i | |x[3]|   |         iad - bd
  //        +-          -+ +-  -+   |    ------------------
  //                                |    (ac-bd) + i(bc+ad)  
  //                                | 

  if (inv) ang =  2*PI/(fp)4; // invert sign for IFFT
  else     ang = -2*PI/(fp)4;
  //
  w[1].r = cos(ang*1); w[1].i = sin(ang*1); // twiddle1 = exp(-2*pi/4 * 1);
  w[2].r = cos(ang*2); w[2].i = sin(ang*2); // twiddle2 = exp(-2*pi/4 * 2);
  w[3].r = cos(ang*3); w[3].i = sin(ang*3); // twiddle3 = exp(-2*pi/4 * 3);

  //         *1       *1       *1       *1
  y[0].r  = x[0].r + x[1].r + x[2].r + x[3].r;
  y[0].i  = x[0].i + x[1].i + x[2].i + x[3].i;
  //         *1       *-i      *-1      *i
  x1[1].r = x[0].r + x[1].i - x[2].r - x[3].i;               
  x1[1].i = x[0].i - x[1].r - x[2].i + x[3].r;               
  //         *1       *-1      *1       *-1
  x1[2].r = x[0].r - x[1].r + x[2].r - x[3].r;
  x1[2].i = x[0].i - x[1].i + x[2].i - x[3].i;
  //         *1       *i       *-1      *-i
  x1[3].r = x[0].r - x[1].i - x[2].r + x[3].i;
  x1[3].i = x[0].i + x[1].r - x[2].i - x[3].r;
  //
  y[1].r = x1[1].r*w[1].r - x1[1].i*w[1].i; // scale radix-4 output
  y[1].i = x1[1].i*w[1].r + x1[1].r*w[1].i;
  //
  y[2].r = x1[2].r*w[2].r - x1[2].i*w[2].i; // scale radix-4 output
  y[2].i = x1[2].i*w[2].r + x1[2].r*w[2].i;
  //
  y[3].r = x1[3].r*w[3].r - x1[3].i*w[3].i; // scale radix-4 output
  y[3].i = x1[3].i*w[3].r + x1[3].r*w[3].i;

  // reorder output stage ... mystery as to why I need this
  if (reorder) {
    temp = y[1].r; 
    y[1].r = -1*y[1].i; 
    y[1].i = temp;
    //
    y[2].r = -1*y[2].r; 
    //
    temp = y[3].r; 
    y[3].r = y[3].i; 
    y[3].i = -1*temp;
  }

  // scale output for inverse FFT
  if (inv) {
    for (i=0; i<4; i++) { // scale output by 1/N for IFFT
      y[i].r = y[i].r/(fp)4;
      y[i].i = y[i].i/(fp)4;
    }
  }

} // r4fft4()

Ich habe gerade einen Radix-4-DIF-FFT von S. Burrus Fortran-Code nach Java portiert. Tatsächlich fehlen mehrere Optimierungen, vor allem der tabellengesteuerte Twiddle-Faktor (Sin- und Cos-Faktoren sollten vorberechnet werden). Dies sollte das fft etwas beschleunigen (vielleicht 50%). Ich muss ein bisschen dafür hacken, aber wenn jemand die richtige Antwort hat, bin ich sehr glücklich und dankbar. Ich werde den optimierten Code so schnell wie möglich veröffentlichen. Ich hoffe, dass ich einige Geschwindigkeitstests gegen den Radix-2-Algorithmus durchführen kann.

Außerdem werden die Multiplikationen mit 1 und sqrt (-1) nicht entfernt. Wenn Sie sie entfernen, wird die Geschwindigkeit etwas erhöht. Aber insgesamt scheint IMHO Radix-4 nicht mehr als 25% schneller zu sein als ein Radix-2, daher weiß ich nicht, ob sich das Verhältnis Geschwindigkeit / Komplexität wirklich lohnt. Denken Sie daran, dass sehr optimierte Bibliotheken wie FFTW weitgehend verfügbar sind und verwendet werden. Diese Bemühungen können also nur eine persönliche "Umleitung" sein!

Hier ist der Java-Code. Das Portieren nach C, C ++ oder C # sollte sehr einfach sein.

public static void FFTR4(double[] X, double[] Y, int N, int M) {
    // N = 4 ^ M
    int N1,N2;
    int I1, I2, I3;
    double CO1,CO2,CO3,SI1,SI2,SI3;
    double A,B,C,E;
    double R1,R2,R3,R4;
    double S1,S2,S3,S4;
    // N = 1 << (M+M);
    N2 = N;
    I2 = 0; I3 = 0;
    for (int K=0; K<M; ++K) {
        N1 = N2;
        N2 = N2 / 4;
        E = PI2 / (double)N1;
        A = 0.0;
        for (int J=0; J < N2; ++J) {
            A = J*E;
            B = A + A;
            C = A + B;
            //Should be pre-calculated for optimization
            CO1 = Math.cos(A);
            CO2 = Math.cos(B);
            CO3 = Math.cos(C);
            SI1 = Math.sin(A);
            SI2 = Math.sin(B);
            SI3 = Math.sin(C);
            for (int I = J; I<N; I+=N1) {
                I1 = I + N2;
                I2 = I1 + N2;
                I3 = I2 + N2;
                R1 = X[I] + X[I2];
                R3 = X[I] - X[I2];
                S1 = Y[I] + Y[I2];
                S3 = Y[I] - Y[I2];
                R2 = X[I1] + X[I3];
                R4 = X[I1] - X[I3];
                S2 = Y[I1] + Y[I3];
                S4 = Y[I1] - Y[I3];
                X[I] = R1 + R2;
                R2 = R1 - R2;
                R1 = R3 - S4;
                R3 = R3 + S4;
                Y[I] = S1 + S2;
                S2 = S1 - S2;
                S1 = S3 + R4;
                S3 = S3 - R4;
                X[I1] = CO1*R3 + SI1*S3;
                Y[I1] = CO1*S3 - SI1*R3;
                X[I2] = CO2*R2 + SI2*S2;
                Y[I2] = CO2*S2 - SI2*R2;
                X[I3] = CO3*R1 + SI3*S1;
                Y[I3] = CO3*S1 - SI3*R1;
            }
        }
    }

    // Radix-4 bit-reverse
    double T;
    int J = 0;
    N2 = N>>2;
    for (int I=0; I < N-1; I++) {
        if (I < J) {
            T = X[I];
            X[I] = X[J];
            X[J] = T;
            T = Y[I];
            Y[I] = Y[J];
            Y[J] = T;
        }
        N1 = N2;
        while ( J >= 3*N1 ) {
            J -= 3*N1;
            N1 >>= 2;
        }
        J += N1;
    }
}

Hier ist der Original Radix-4 DIF FORTRAN Code von Sidney Burrus:

Radix-4, DIF, One Butterfly FFT

Erstens ist Ihr angeblicher "Radix-4-Schmetterling" eine 4-Punkt-DFT, keine FFT. Es hat 16 komplexe Operationen (dh: N im Quadrat). Eine typische 4-Punkt-FFT hätte nur Nlog (Basis 2) N (= 8 für N = 4). Zweitens haben Sie einige vermeintliche Skalierungsfaktoren w [] .r und w [] .i, die nicht dazu gehören. Vielleicht haben Sie sie von einem Radix-4-Schmetterling erhalten, der in einer größeren Grafik dargestellt ist. An einen solchen Schmetterling wären einige Zwischenstufen-Twiddles angehängt, aber sie sind eigentlich nicht Teil des Schmetterlings. Eine 4-Punkt-FFT hat nur einen internen Schmetterling von -j, wenn sie für eine negative Exponenten-FFT ausgelegt ist.

Anstatt zu versuchen, Ihren Code zu reparieren, ist es genauso einfach, meinen eigenen zu schreiben, wie unten gezeigt (DevC ++ - Compiler; Ausgaben am Ende des Codes angehängt):

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;
void fft4(double* r, double* i);    // prototype declaration
int main (int nNumberofArgs, char* pszArgs[ ] ) { // arguments needed for Dev C++ I/O

double r[4] = {1.5, -2.3, 4.65, -3.51}, i[4] = {-1.0, 2.6, 3.75, -2.32} ;
long n, k, j;      double  yr[4] = {0.}, yi[4] = {0.};
double ang, C, S, twopi = 6.2831853071795865;

cout<<"\n original real/imag data";
cout<<"\n n         r[n]            i[n]\n";
for (n = 0; n < 4; n++)  {
    printf("%2d\t%9.4f\t%9.4f\n",n,r[n],i[n]);
} //end for loop over n

// 4 point DFT
for (k = 0; k < 4; k++) {
    ang = twopi*k/4;
    for (j = 0; j < 4; j++) {
        C = cos(j*ang);       S = sin(j*ang);
        yr[k] = yr[k] + r[j]*C + i[j]*S;   // ( C - jS )*( r + ji )
        yi[k] = yi[k] + i[j]*C - r[j]*S;   // = ( rC + iS ) + j( iC - rS )
    }
}

cout<<"\n 4 point DFT results";
cout<<"\n n         yr[n]           yi[n]           amplitude       phase(radians)\n";
double amp, phase;
for (n = 0; n < 4; n++)  {
    yr[n] = yr[n]/4 ;      yi[n] = yi[n]/4 ;  // scale outputs
    amp = sqrt( yr[n]*yr[n] + yi[n]*yi[n] ) ;
    phase = atan2( yi[n], yr[n] ) ; 
    printf("%2d\t%9.4f\t%9.4f\t%9.4f\t%9.4f\n",n,yr[n],yi[n],amp,phase);
} //end for loop over n

fft4(r, i) ;

cout<<"\n 4 point FFT results";
cout<<"\n n         r[n]            i[n]            amplitude       phase(radians)\n";

for (n = 0; n < 4; n++)  {
    r[n] = r[n]/4 ;      i[n] = i[n]/4 ;  // scale outputs
    amp = sqrt( r[n]*r[n] + i[n]*i[n] ) ;
    phase = atan2( i[n], r[n] ) ; 
    printf("%2d\t%9.4f\t%9.4f\t%9.4f\t%9.4f\n",n,r[n],i[n],amp,phase);
} //end for loop over n

fft4(i, r); // this is an inverse FFT (complex in/out routine)

cout<<"\n 4 point inverse FFT results";
cout<<"\n n         r[n]            i[n]\n";
for (n = 0; n < 4; n++)  {
    printf("%2d\t%9.4f\t%9.4f\n",n,r[n],i[n]);
} //end for loop over n

system ("PAUSE");
return 0;
} // end main
//************************ fft4 **********
void fft4(double* r, double* i) {
double t;

t = r[0]; r[0] = t + r[2]; r[2] = t - r[2];
t = i[0]; i[0] = t + i[2]; i[2] = t - i[2];
t = r[1]; r[1] = t + r[3]; r[3] = t - r[3];
t = i[1]; i[1] = t + i[3]; i[3] = t - i[3];

t = r[3]; r[3] = i[3]; i[3] = -t; // (r + ji)*(-j)

t = r[0]; r[0] = t + r[1]; r[1] = t - r[1];
t = i[0]; i[0] = t + i[1]; i[1] = t - i[1];
t = r[2]; r[2] = t + r[3]; r[3] = t - r[3];
t = i[2]; i[2] = t + i[3]; i[3] = t - i[3];

t = r[1]; r[1] = r[2]; r[2] = t;  // swap 1
t = i[1]; i[1] = i[2]; i[2] = t;  //  and 2
} // end fft4




 original real/imag data
 n         r[n]            i[n]
 0         1.5000         -1.0000
 1        -2.3000          2.6000
 2         4.6500          3.7500
 3        -3.5100         -2.3200

 4 point DFT results
 n         yr[n]           yi[n]           amplitude       phase(radians)
 0         0.0850          0.7575          0.7623          1.4591
 1         0.4425         -1.4900          1.5543         -1.2821
 2         2.9900          0.6175          3.0531          0.2037
 3        -2.0175         -0.8850          2.2031         -2.7282

 4 point FFT results
 n         r[n]            i[n]            amplitude       phase(radians)
 0         0.0850          0.7575          0.7623          1.4591
 1         0.4425         -1.4900          1.5543         -1.2821
 2         2.9900          0.6175          3.0531          0.2037
 3        -2.0175         -0.8850          2.2031         -2.7282

 4 point inverse FFT results
 n         r[n]            i[n]
 0         1.5000         -1.0000
 1        -2.3000          2.6000
 2         4.6500          3.7500
 3        -3.5100         -2.3200

Zunächst werden die Eingabedaten (4 reelle, 4 imaginäre) ausgedruckt. Dann wird eine 4-Punkt-DFT genommen. Die Ergebnisse (yr [] und yi [] plus Ampere / Phase) werden ausgedruckt. Da die Originaldaten r [] und i [] bei der DFT nicht überschrieben wurden, werden diese Eingaben als Eingaben für die 4-Punkt-FFT wiederverwendet. Beachten Sie, dass letztere weniger +/- Operationen als die DFT hat.

Der Code für die FFT ist weder besonders elegant noch effizient - es gibt viele Möglichkeiten, Schmetterlinge zu machen. Der obige Code entspricht den vier Radix-2-Schmetterlingen, die in Rabiner und Golds Buch „Theorie und Anwendung der digitalen Signalverarbeitung“ (S. 580, Abb. 10.9) gezeigt sind Zahlen im Buch waren positiv). Beachten Sie, dass der Code nur ein Twiddle von -j enthält und dies keine Multiplikation erfordert (es handelt sich um einen Swap- / Vorzeichenwechsel).

Nach der FFT werden die Ergebnisse gedruckt. Sie sind die gleichen wie die DFT

Und schließlich werden die skalierten Ergebnisse der FFT als Eingaben für eine inverse FFT verwendet. Dies wird über die Methode 'Austausch' oder 'Liste umkehren' erreicht (dh wenn FFT (r, i) eine Vorwärts-FFT ist, dann ist FFT (i, r) eine Umkehrung - vorausgesetzt natürlich, dass die FFT ist fähig, komplexe Ein- / Ausgänge zu handhaben - mit anderen Worten - keine "nur realen" Routinen, die normalerweise davon ausgehen, dass imaginäre Eingänge Null sind). Diese Methode wurde vor fast 25 Jahren beschrieben in:

P. Duhamel, B. Piron, JM Etcheto, "Über die Berechnung der inversen DFT", IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 36, Februar 1988, S. 285-286.

Das Ergebnis der Umkehrung wird dann ausgedruckt. Dies entspricht den ursprünglichen Eingabedaten.