mirror/kissfft
1
0
Fork 0
mirror of https://github.com/mborgerding/kissfft.git synced 2025-06-04 01:28:23 -04:00
Code Issues Packages Projects Releases Wiki Activity

better, reduced mutliply in real fft

Browse Source
This commit is contained in:
Mark Borgerding 2003-12-10 05:48:42 +00:00
parent b880c599ef
commit 8a992d0e98
3 changed files with 31 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
fft.py
View File

@ -77,10 +77,11 @@ def real_fft( f,inv ):
F1k = Fp[k] + Fp[N-k].conjugate() F1k = Fp[k] + Fp[N-k].conjugate()
F2k = Fp[k] - Fp[N-k].conjugate() F2k = Fp[k] - Fp[N-k].conjugate()
F[k] = ( F1k + tw * F2k ) * .5 F2k *= tw
F[k] = ( F1k + F2k ) * .5
F[N-k] = ( F1k - F2k ).conjugate() * .5
#F[N-k] = ( F1kp + e ** ( -j*pi*(.5+float(N-k)/N ) ) * F2kp ) * .5 #F[N-k] = ( F1kp + e ** ( -j*pi*(.5+float(N-k)/N ) ) * F2kp ) * .5
F[N-k] = ( F1k.conjugate() - tw.conjugate() * F2k.conjugate() ) * .5 #F[N-k] = ( F1k.conjugate() - tw.conjugate() * F2k.conjugate() ) * .5
F[N] = complex( Fp[0].real - Fp[0].imag , 0 ) F[N] = complex( Fp[0].real - Fp[0].imag , 0 )
return F return F

15
kiss_fftr.c
View File

@ -79,20 +79,25 @@ void kiss_fftr(const void * cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *fr
freqdata[0].r = st->tmpbuf[0].r + st->tmpbuf[0].i; freqdata[0].r = st->tmpbuf[0].r + st->tmpbuf[0].i;
freqdata[0].i = 0; freqdata[0].i = 0;
for (k=1;k<N;++k) { for (k=1;k <= N/2 ; ++k ) {
kiss_fft_cpx fpnk,fpk,f1k; kiss_fft_cpx fpnk,fpk,f1k,f2k,tw;
fpk = st->tmpbuf[k]; fpk = st->tmpbuf[k];
fpnk.r = st->tmpbuf[N-k].r; fpnk.r = st->tmpbuf[N-k].r;
fpnk.i = -st->tmpbuf[N-k].i; fpnk.i = -st->tmpbuf[N-k].i;
C_ADD( f1k, fpk , fpnk ); C_ADD( f1k, fpk , fpnk );
C_SUBFROM( fpk , fpnk ); C_SUB( f2k, fpk , fpnk );
C_MUL( freqdata[k], fpk , st->super_twiddles[k] ); C_MUL( tw , f2k , st->super_twiddles[k]);
C_ADDTO(freqdata[k],f1k);
C_ADD( freqdata[k] , f1k ,tw);
freqdata[k].r /= 2; freqdata[k].r /= 2;
freqdata[k].i /= 2; freqdata[k].i /= 2;
freqdata[N-k].r = f1k.r - tw.r;
freqdata[N-k].i = - (f1k.i - tw.i);
freqdata[N-k].r /= 2;
freqdata[N-k].i /= 2;
} }
freqdata[N].r = st->tmpbuf[0].r - st->tmpbuf[0].i; freqdata[N].r = st->tmpbuf[0].r - st->tmpbuf[0].i;
freqdata[N].i = 0; freqdata[N].i = 0;

34
test/test_real.c
View File

@ -34,9 +34,7 @@ double snr_compare( kiss_fft_cpx * vec1,kiss_fft_cpx * vec2, int n)
printf( "\npoor snr, try a scaling factor %f\n" , scale ); printf( "\npoor snr, try a scaling factor %f\n" , scale );
return snr; return snr;
} }
#define RANDOM #define RANDOM
#ifndef RANDOM #ifndef RANDOM
#define NFFT 8 #define NFFT 8
#else #else
@ -80,6 +78,22 @@ int main()
kiss_fftr(kiss_fftr_state,sin,sout); kiss_fftr(kiss_fftr_state,sin,sout);
printf( "nfft=%d, inverse=%d, snr=%g\n", printf( "nfft=%d, inverse=%d, snr=%g\n",
NFFT,0, snr_compare(cout,sout,(NFFT/2)+1) ); NFFT,0, snr_compare(cout,sout,(NFFT/2)+1) );
#ifdef RANDOM
ts = cputime();
for (i=0;i<NUMFFTS;++i) {
kiss_fft(kiss_fft_state,cin,cout);
}
tfft = cputime() - ts;
ts = cputime();
for (i=0;i<NUMFFTS;++i) {
kiss_fftr( kiss_fftr_state, sin, cout );
/* kiss_fftri(kiss_fftr_state,cout,sin); */
}
trfft = cputime() - ts;
printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft);
#endif
free(kiss_fft_state); free(kiss_fft_state);
free(kiss_fftr_state); free(kiss_fftr_state);
@ -99,22 +113,6 @@ int main()
printf( "nfft=%d, inverse=%d, snr=%g\n", printf( "nfft=%d, inverse=%d, snr=%g\n",
NFFT,1, snr_compare(cin,cin,NFFT/2) ); NFFT,1, snr_compare(cin,cin,NFFT/2) );
#ifdef RANDOM
ts = cputime();
for (i=0;i<NUMFFTS;++i) {
kiss_fft(kiss_fft_state,cin,cout);
}
tfft = cputime() - ts;
ts = cputime();
for (i=0;i<NUMFFTS;++i) {
/* kiss_fftr(kiss_fftr_state,sin,cout); */
kiss_fftri(kiss_fftr_state,cout,sin);
}
trfft = cputime() - ts;
printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft);
#endif
free(kiss_fft_state); free(kiss_fft_state);
free(kiss_fftr_state); free(kiss_fftr_state);