radix-4 FFT-Implementierung

  • Ich habe eine 4-Punkt-Radix-4-FFT implementiert und festgestellt, dass ich die Ausgabebegriffe ein wenig manipulieren muss, damit sie mit einem dft übereinstimmt

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

     //                                | 
    // 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 irgendjemand erkennen, wo ich falsch liege?

    Danke,

    -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()
     
    27 July 2012
    user1211582
2 answers
  • Erstens ist Ihr vermeintlicher "Radix-4-Schmetterling" eine 4-Punkt-DFT, keine FFT. Es hat 16 komplexe Operationen (dh: N-Quadrat). Eine typische 4-Punkt-FFT hätte nur Nlog (Basis 2) N (= 8 für N = 4). Zweitens haben Sie einige vermeintliche w [] .r und w [] .i 'Skalierungsfaktoren', die nicht dazugehören. Vielleicht haben Sie sie von einem Schmetterlings-4-Schmetterling erhalten, der in einer größeren Grafik dargestellt ist. Ein solcher Schmetterling hätte einige Zwischendrehungen an der Bühne, die jedoch nicht Teil des Schmetterlings sind. Eine 4-Punkt-FFT hat nur einen internen Schmetterling von -j, wenn sie für eine negative Exponent-FFT entwickelt wurde.

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

     #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
     

    Zuerst werden die Eingabedaten ( 4 real, 4 imaginär) wird ausgedruckt. Dann wird eine 4-Punkt-DFT genommen. Die Ergebnisse (yr [] und yi [] plus Amp / Phase) werden ausgedruckt. Da die Originaldaten r [] und i [] bei der DFT nicht überschrieben wurden, werden diese Eingänge als Eingänge für die 4-Punkt-FFT verwendet. 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 Basispunkten von Radix-2, die in Rabiner und Golds Buch "Theorie und Anwendung der digitalen Signalverarbeitung" (S. 580, Abb. 10.9) dargestellt sind Zahl im Buch war positiv). Beachten Sie, dass im Code nur ein Twiddle von -j vorhanden ist, und dies erfordert keine Multiplikation (es handelt sich hierbei um einen Swap / Vorzeichenwechsel).

    Nach der FFT werden die Ergebnisse angezeigt gedruckt. Sie sind mit der DFT identisch.

    Und schließlich werden die skalierten Ergebnisse der FFT als Eingaben in eine inverse FFT verwendet. Dies wird über die 'exchange'- oder' reverse the list'-Methode erreicht (dh: Wenn FFT (r, i) eine Vorwärts-FFT ist, dann ist FFT (i, r) eine Inverse-Provision

    04 November 2012
    Željko Živković
  • Ich habe gerade eine radix-4-DIF-FFT von S. Burrus Fortran-Code nach Java portiert. Tatsächlich fehlen einige Optimierungen, vor allem der tischgesteuerte Twiddle-Faktor (Sin- und Cos-Faktoren sollten dies sein) vorberechnet). Dies sollte die 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 posten. Ich hoffe, mit einigen Geschwindigkeitstests gegen den Radix-2-Algorithmus.

    Mehr, die Multiplikationen von 1 und sqrt (-1) sind nicht vorhanden nicht entfernt Wenn Sie sie entfernen, werden Sie etwas schneller. Aber insgesamt scheint IMHO radix-4 nicht mehr als 25% schneller als ein radix-2 zu sein, daher weiß ich nicht, ob das Verhältnis zwischen Geschwindigkeit und Komplexität wirklich wert ist. Denken Sie daran, dass sehr optimierte Bibliotheken wie FFTW weitgehend verfügbar sind und verwendet werden, so dass dieser Aufwand nur eine persönliche "Divertissment" sein kann.

    Hier ist der Java-Code. Die Portierung 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 ursprüngliche Radix-4-DIF FORTRAN-Code von Sidney Burrus:

    Radix-4, DIF , Eine Schmetterlings-FFT

    20 November 2012
    Yozek