mise en œuvre de radix-4 FFT

J'ai implémenté une FFT radix-4 à 4 points et j'ai constaté que je devais manipuler les termes de sortie pour l'adapter à un dft.

Mon code est une implémentation assez directe de la formulation matricielle, donc je ne sais pas quel est le problème

//                                | 
// 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)  
//                                | 

Quelqu'un peut-il repérer où je me suis trompé?

Merci,

-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()

Je viens de porter un fft radix-4 DIF du code S. Burrus Fortran vers Java. En fait, il manque plusieurs optimisations, tout d'abord le facteur twiddle piloté par table (les facteurs sin et cos doivent être pré-calculés). Cela devrait accélérer le fft un peu plus (peut-être 50%). Je dois pirater un peu pour cela, mais si quelqu'un a la bonne réponse, je serai très heureux et reconnaissant. Je posterai le code optimisé dès que possible, j'espère peut-être avec des tests de vitesse vs l'algorithme radix-2.

De plus, les multiplications par 1 et sqrt (-1) ne sont pas supprimées. Les supprimer accélérera un peu plus. Mais globalement IMHO radix-4 ne semble pas être plus de 25% plus rapide qu'un radix-2, donc je ne sais pas si le rapport vitesse / complexité vaut vraiment la peine. Gardez à l'esprit que des bibliothèques très optimisées comme FFTW sont largement disponibles et utilisées, donc cet effort pourrait être juste un «détournement» personnel!

Voici le code java. Le porter en C, C ++ ou C # devrait être très facile.

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;
    }
}

Voici le code original Radix-4 DIF FORTRAN de Sidney Burrus:

Radix-4, DIF, un papillon FFT

Tout d'abord, votre supposé «papillon radix-4» est un DFT à 4 points, pas un FFT. Il a 16 opérations complexes (ie: N au carré). Une FFT typique à 4 points aurait seulement Nlog (base 2) N (= 8 pour N = 4). Deuxièmement, vous avez des facteurs d'échelle supposés w [] .r et w [] .i qui n'appartiennent pas. Peut-être les avez-vous obtenus à partir d'un papillon Radix-4 montré dans un graphique plus grand. Un tel papillon aurait des twiddles inter-étages, mais ils ne font pas partie du papillon. Une FFT à 4 points n'a qu'un papillon interne de -j lorsqu'elle est conçue pour une FFT à exposant négatif.

Plutôt que d'essayer de corriger votre code, il est tout aussi facile d'écrire le mien, comme indiqué ci-dessous (compilateur DevC ++; sorties ajoutées à la fin du code):

#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

Tout d'abord, les données d'entrée (4 réelles, 4 imaginaires) sont imprimées. Ensuite, une DFT à 4 points est prise. Les résultats (yr [] et yi [] plus ampli / phase) sont imprimés. Étant donné que les données d'origine r [] et i [] n'ont pas été écrasées lors de l'exécution de la DFT, ces entrées sont réutilisées en tant qu'entrées dans la FFT à 4 points. Notez que ce dernier a moins d'opérations +/- que le DFT.

Le code de la FFT n'est pas particulièrement élégant ni efficace - il existe de nombreuses façons de faire des papillons. Le code ci-dessus correspond aux quatre papillons radix-2 présentés dans le livre de Rabiner et Gold «Théorie et application du traitement numérique du signal» (p. 580, Fig. 10.9), avec des twiddles modifiés pour refléter un exposant négatif (ceux utilisés pour le chiffre dans le livre étaient positifs). Notez qu'il n'y a qu'un seul twiddle de -j dans le code, et cela ne nécessite pas de multiplication (c'est un changement de swap / signe).

Après la FFT, les résultats sont imprimés. Ils sont les mêmes que le DFT

Et enfin, les résultats mis à l'échelle de la FFT sont utilisés comme entrées d'une FFT inverse. Ceci est accompli via la méthode «échange» ou «inverser la liste» (c'est-à-dire: si la FFT (r, i) est une FFT vers l'avant, alors la FFT (i, r) est un inverse - à condition, bien sûr, que la FFT soit capable de gérer des entrées / sorties complexes - en d'autres termes - pas de routines «réelles uniquement», qui supposent généralement que les entrées imaginaires sont nulles). Cette méthode a été décrite il y a près de 25 ans dans:

P. Duhamel, B. Piron, JM Etcheto, «On Computing the Inverse DFT», IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 36, février 1988, p. 285-286.

Le résultat de l'inverse est ensuite imprimé. C'est la même chose que les données d'entrée d'origine.