From ba105600b458d0f35667958c2caf24e3b8cf955b Mon Sep 17 00:00:00 2001 From: Mark Borgerding Date: Mon, 9 May 2005 01:43:49 +0000 Subject: [PATCH] getting ready for v 1_2_2 --- CHANGELOG | 9 +++++---- _kiss_fft_guts.h | 41 ++++++++++++++++++++++++++++++++++++----- tools/kiss_fftr.c | 14 +++++++++----- 3 files changed, 50 insertions(+), 14 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index 16a4f83..95c86b6 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,9 +1,10 @@ 1.2.2 (May 6, 2005) The Matthew release - Replaced fixed point division with multiply&shift. Thanks to Jean-Marc Valin for the prodding. - - Added C_FIXDIV calls to real fft routines to prevent possible overflow when using fixed point. - Credit to Robert Oschler of robodance for pointing me towards the bug. + Replaced fixed point division with multiply&shift. Thanks to Jean-Marc Valin for + discussions regarding. Considerable speedup. + Corrected overflow protection in real fft routines when using fixed point. + Finder's Credit goes to Robert Oschler of robodance for pointing me at the bug. + This also led to the CHECK_OVERFLOW_OP macro. 1.2.1 (April 4, 2004) compiles cleanly with just about every -W warning flag under the sun diff --git a/_kiss_fft_guts.h b/_kiss_fft_guts.h index d61bbce..f295254 100644 --- a/_kiss_fft_guts.h +++ b/_kiss_fft_guts.h @@ -43,6 +43,13 @@ struct kiss_fft_state{ * */ #ifdef FIXED_POINT +#if defined(CHECK_OVERFLOW) +# define CHECK_OVERFLOW_OP(a,op,b) \ + if ( (long)(a) op (long)(b) > 32767 || (long)(a) op (long)(b) < -32768 ) { \ + fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(long)(a) op (long)(b) ); } +#endif + + # define smul(a,b) ( (long)(a)*(b) ) # define sround( x ) (short)( ( (x) + (1<<14) ) >>15 ) @@ -53,7 +60,7 @@ struct kiss_fft_state{ (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0) # define DIVSCALAR(x,k) \ - (x) = sround( smul( x, 32768/k ) ) + (x) = sround( smul( x, 32767/k ) ) # define C_FIXDIV(c,div) \ do { DIVSCALAR( (c).r , div); \ @@ -75,14 +82,38 @@ struct kiss_fft_state{ (c).i *= (s); }while(0) #endif +#ifndef CHECK_OVERFLOW_OP +# define CHECK_OVERFLOW_OP(a,op,b) /* noop */ +#endif + #define C_ADD( res, a,b)\ - do { (res).r=(a).r+(b).r; (res).i=(a).i+(b).i; }while(0) + do { \ + CHECK_OVERFLOW_OP((a).r,+,(b).r)\ + CHECK_OVERFLOW_OP((a).i,+,(b).i)\ + (res).r=(a).r+(b).r; (res).i=(a).i+(b).i; \ + }while(0) #define C_SUB( res, a,b)\ - do { (res).r=(a).r-(b).r; (res).i=(a).i-(b).i; }while(0) + do { \ + CHECK_OVERFLOW_OP((a).r,-,(b).r)\ + CHECK_OVERFLOW_OP((a).i,-,(b).i)\ + (res).r=(a).r-(b).r; (res).i=(a).i-(b).i; \ + }while(0) #define C_ADDTO( res , a)\ - do { (res).r += (a).r; (res).i += (a).i; }while(0) + do { \ + CHECK_OVERFLOW_OP((res).r,+,(a).r)\ + CHECK_OVERFLOW_OP((res).i,+,(a).i)\ + (res).r += (a).r; (res).i += (a).i;\ + }while(0) + #define C_SUBFROM( res , a)\ - do { (res).r -= (a).r; (res).i -= (a).i; }while(0) + do {\ + CHECK_OVERFLOW_OP((res).r,-,(a).r)\ + CHECK_OVERFLOW_OP((res).i,-,(a).i)\ + (res).r -= (a).r; (res).i -= (a).i; \ + }while(0) + + + static void kf_cexp(kiss_fft_cpx * x,double phase) /* returns e ** (j*phase) */ diff --git a/tools/kiss_fftr.c b/tools/kiss_fftr.c index 681d242..859ac90 100644 --- a/tools/kiss_fftr.c +++ b/tools/kiss_fftr.c @@ -65,6 +65,7 @@ void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *fr { /* input buffer timedata is stored row-wise */ int k,N; + kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc; if ( st->substate->inverse) { fprintf(stderr,"kiss fft usage error: improper alloc\n"); @@ -76,12 +77,15 @@ void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *fr /*perform the parallel fft of two real signals packed in real,imag*/ kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf ); - freqdata[0].r = st->tmpbuf[0].r + st->tmpbuf[0].i; + tdc.r = st->tmpbuf[0].r; + tdc.i = st->tmpbuf[0].i; + C_FIXDIV(tdc,2); + + CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i); + freqdata[0].r = tdc.r + tdc.i; freqdata[0].i = 0; - C_FIXDIV(freqdata[0],2); for (k=1;k <= N/2 ; ++k ) { - kiss_fft_cpx fpnk,fpk,f1k,f2k,tw; fpk = st->tmpbuf[k]; fpnk.r = st->tmpbuf[N-k].r; @@ -100,9 +104,9 @@ void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *fr freqdata[N-k].r = (f1k.r - tw.r)/2; freqdata[N-k].i = - (f1k.i - tw.i)/2; } - freqdata[N].r = st->tmpbuf[0].r - st->tmpbuf[0].i; + CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i); + freqdata[N].r = tdc.r - tdc.i; freqdata[N].i = 0; - C_FIXDIV(freqdata[N],2); } void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)