implementasi FFT radix-4

Saya menerapkan 4-point radix-4 FFT dan menemukan bahwa saya perlu melakukan beberapa manipulasi persyaratan output untuk mendapatkannya agar cocok dengan DFT.

Kode saya adalah implementasi langsung dari formulasi matriks jadi saya tidak jelas tentang apa masalahnya

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

Adakah yang bisa melihat kesalahan saya?

Terima kasih,

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

Saya baru saja porting radix-4 DIF fft dari kode S. Burrus Fortran ke Jawa. Sebenarnya itu tidak memiliki beberapa optimasi, pertama-tama faktor twiddle didorong oleh tabel (faktor dosa dan cos harus dihitung sebelumnya). Ini akan mempercepat fft lebih banyak (mungkin 50%). Saya harus sedikit meretas untuk itu tetapi jika seseorang memiliki jawaban yang benar saya akan sangat senang dan berterima kasih. Saya akan memposting kode yang dioptimalkan secepatnya saya harap mungkin dengan beberapa tes kecepatan vs radix-2 algoritma.

Terlebih lagi, perkalian dengan 1 dan sqrt (-1) tidak dihapus. Menghapusnya akan mempercepat sedikit lagi. Tapi secara keseluruhan IMHO radix-4 tampaknya tidak lebih dari 25% lebih cepat daripada radix-2, jadi saya tidak tahu apakah rasio kecepatan / kompleksitas benar-benar bernilai. Ingatlah bahwa perpustakaan yang sangat dioptimalkan seperti FFTW sebagian besar tersedia dan digunakan, sehingga upaya ini bisa jadi hanya 'pengalihan' pribadi!

Ini adalah kode java. Porting ke C, C ++ atau C # seharusnya sangat mudah.

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

Berikut adalah kode asli Radix-4 DIF FORTRAN oleh Sidney Burrus:

Radix-4, DIF, One Butterfly FFT

Pertama, 'kupu-kupu radix-4' Anda seharusnya adalah DFT 4 poin, bukan FFT. Ini memiliki 16 operasi kompleks (yaitu: N kuadrat). FFT 4 poin biasanya hanya memiliki Nlog (basis 2) N (= 8 untuk N = 4). Kedua, Anda memiliki beberapa faktor w [] .r dan w [] .i 'skala' yang tidak termasuk. Mungkin Anda mendapatkannya dari kupu-kupu radix-4 yang ditampilkan dalam grafik yang lebih besar. Seekor kupu-kupu seperti itu akan memiliki beberapa rumbai interstage ditambahkan ke dalamnya, tetapi mereka sebenarnya bukan bagian dari kupu-kupu. Sebuah 4 poin FFT hanya memiliki kupu-kupu internal -j ketika dirancang untuk FFT eksponen negatif.

Daripada mencoba memperbaiki kode Anda, lebih mudah untuk menulis sendiri, seperti yang ditunjukkan di bawah ini (kompilator DevC ++; output ditambahkan di akhir kode):

#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

Pertama, data input (4 nyata, 4 imajiner) dicetak. Kemudian DFT 4 poin diambil. Hasilnya (yr [] dan yi [] plus amp / fase) dicetak. Karena r [] dan i [] data asli tidak ditimpa ketika melakukan DFT, input tersebut digunakan kembali sebagai input ke 4 titik FFT. Perhatikan bahwa yang terakhir memiliki lebih sedikit operasi +/- daripada DFT.

Kode untuk FFT tidak terlalu elegan dan juga tidak efisien - ada banyak cara untuk melakukan kupu-kupu. Kode di atas sesuai dengan empat kupu-kupu radix-2 yang ditunjukkan dalam buku Rabiner dan Gold "Teori dan Aplikasi Pemrosesan Sinyal Digital" (hlm. 580, Gambar 10.9), dengan twiddle yang dimodifikasi untuk mencerminkan eksponen negatif (yang digunakan untuk angka dalam buku itu positif). Perhatikan bahwa hanya ada satu twiddle dari -j dalam kode, dan ini tidak memerlukan multiply (ini adalah perubahan swap / tanda).

Setelah FFT, hasilnya dicetak. Mereka sama dengan DFT

Dan akhirnya, hasil skala dari FFT digunakan sebagai input ke FFT terbalik. Ini dicapai melalui metode 'pertukaran' atau 'membalikkan daftar' (yaitu: jika FFT (r, i) adalah FFT maju, maka FFT (i, r) adalah kebalikan - asalkan, tentu saja, bahwa FFT adalah mampu menangani input / output yang kompleks - dengan kata lain - tidak ada rutinitas 'nyata saja', yang biasanya menganggap bahwa input imajiner adalah nol). Metode ini dijelaskan hampir 25 tahun yang lalu di:

P. Duhamel, B. Piron, JM Etcheto, "Pada Komputasi DFT Invers," Transaksi IEEE pada Akustik, Pidato dan Pemrosesan Sinyal, vol. 36, Februari 1988, hlm. 285-286.

Hasil kebalikannya kemudian dicetak. Ini sama dengan data input asli.