[xiph-commits] r12645 - trunk/speex/libspeex

jm at svn.xiph.org jm at svn.xiph.org
Mon Mar 5 04:32:19 PST 2007


Author: jm
Date: 2007-03-05 04:32:16 -0800 (Mon, 05 Mar 2007)
New Revision: 12645

Modified:
   trunk/speex/libspeex/fftwrap.c
   trunk/speex/libspeex/kiss_fftr.c
   trunk/speex/libspeex/kiss_fftr.h
Log:
Removed one step of unnecessary copying of data and instead changed the 
packing of the kiss_fft real fft.


Modified: trunk/speex/libspeex/fftwrap.c
===================================================================
--- trunk/speex/libspeex/fftwrap.c	2007-03-05 09:14:02 UTC (rev 12644)
+++ trunk/speex/libspeex/fftwrap.c	2007-03-05 12:32:16 UTC (rev 12645)
@@ -169,14 +169,7 @@
    int shift;
    struct kiss_config *t = (struct kiss_config *)table;
    shift = maximize_range(in, in, 32000, t->N);
-   kiss_fftr(t->forward, in, t->freq_data);
-   out[0] = t->freq_data[0].r;
-   for (i=1;i<t->N>>1;i++)
-   {
-      out[(i<<1)-1] = t->freq_data[i].r;
-      out[(i<<1)] = t->freq_data[i].i;
-   }
-   out[(i<<1)-1] = t->freq_data[i].r;
+   kiss_fftr2(t->forward, in, out);
    renorm_range(in, in, shift, t->N);
    renorm_range(out, out, shift, t->N);
 }
@@ -189,14 +182,9 @@
    float scale;
    struct kiss_config *t = (struct kiss_config *)table;
    scale = 1./t->N;
-   kiss_fftr(t->forward, in, t->freq_data);
-   out[0] = scale*t->freq_data[0].r;
-   for (i=1;i<t->N>>1;i++)
-   {
-      out[(i<<1)-1] = scale*t->freq_data[i].r;
-      out[(i<<1)] = scale*t->freq_data[i].i;
-   }
-   out[(i<<1)-1] = scale*t->freq_data[i].r;
+   kiss_fftr2(t->forward, in, out);
+   for (i=0;i<t->N;i++)
+      out[i] *= scale;
 }
 #endif
 
@@ -204,17 +192,7 @@
 {
    int i;
    struct kiss_config *t = (struct kiss_config *)table;
-   t->freq_data[0].r = in[0];
-   t->freq_data[0].i = 0;
-   for (i=1;i<t->N>>1;i++)
-   {
-      t->freq_data[i].r = in[(i<<1)-1];
-      t->freq_data[i].i = in[(i<<1)];
-   }
-   t->freq_data[i].r = in[(i<<1)-1];
-   t->freq_data[i].i = 0;
-
-   kiss_fftri(t->backward, t->freq_data, out);
+   kiss_fftri2(t->backward, in, out);
 }
 
 

Modified: trunk/speex/libspeex/kiss_fftr.c
===================================================================
--- trunk/speex/libspeex/kiss_fftr.c	2007-03-05 09:14:02 UTC (rev 12644)
+++ trunk/speex/libspeex/kiss_fftr.c	2007-03-05 12:32:16 UTC (rev 12645)
@@ -132,7 +132,7 @@
     }
 }
 
-void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
+void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata, kiss_fft_scalar *timedata)
 {
     /* input buffer timedata is stored row-wise */
     int k, ncfft;
@@ -168,3 +168,91 @@
     }
     kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata);
 }
+
+void kiss_fftr2(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_scalar *freqdata)
+{
+   /* input buffer timedata is stored row-wise */
+   int k,ncfft;
+   kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc;
+
+   if ( st->substate->inverse) {
+      speex_error("kiss fft usage error: improper alloc\n");
+   }
+
+   ncfft = st->substate->nfft;
+
+   /*perform the parallel fft of two real signals packed in real,imag*/
+   kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf );
+    /* The real part of the DC element of the frequency spectrum in st->tmpbuf
+   * contains the sum of the even-numbered elements of the input time sequence
+   * The imag part is the sum of the odd-numbered elements
+   *
+   * The sum of tdc.r and tdc.i is the sum of the input time sequence. 
+   *      yielding DC of input time sequence
+   * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... 
+   *      yielding Nyquist bin of input time sequence
+    */
+ 
+   tdc.r = st->tmpbuf[0].r;
+   tdc.i = st->tmpbuf[0].i;
+   C_FIXDIV(tdc,2);
+   CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i);
+   CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i);
+   freqdata[0] = tdc.r + tdc.i;
+   freqdata[2*ncfft-1] = tdc.r - tdc.i;
+
+   for ( k=1;k <= ncfft/2 ; ++k ) {
+      fpk    = st->tmpbuf[k]; 
+      fpnk.r =   st->tmpbuf[ncfft-k].r;
+      fpnk.i = - st->tmpbuf[ncfft-k].i;
+      C_FIXDIV(fpk,2);
+      C_FIXDIV(fpnk,2);
+
+      C_ADD( f1k, fpk , fpnk );
+      C_SUB( f2k, fpk , fpnk );
+      C_MUL( tw , f2k , st->super_twiddles[k]);
+
+      freqdata[2*k-1] = HALF_OF(f1k.r + tw.r);
+      freqdata[2*k] = HALF_OF(f1k.i + tw.i);
+      freqdata[2*(ncfft-k)-1] = HALF_OF(f1k.r - tw.r);
+      freqdata[2*(ncfft-k)] = HALF_OF(tw.i - f1k.i);
+   }
+}
+
+void kiss_fftri2(kiss_fftr_cfg st,const kiss_fft_scalar *freqdata,kiss_fft_scalar *timedata)
+{
+   /* input buffer timedata is stored row-wise */
+   int k, ncfft;
+
+   if (st->substate->inverse == 0) {
+      speex_error ("kiss fft usage error: improper alloc\n");
+   }
+
+   ncfft = st->substate->nfft;
+
+   st->tmpbuf[0].r = freqdata[0] + freqdata[2*ncfft-1];
+   st->tmpbuf[0].i = freqdata[0] - freqdata[2*ncfft-1];
+   /*C_FIXDIV(st->tmpbuf[0],2);*/
+
+   for (k = 1; k <= ncfft / 2; ++k) {
+      kiss_fft_cpx fk, fnkc, fek, fok, tmp;
+      fk.r = freqdata[2*k-1];
+      fk.i = freqdata[2*k];
+      fnkc.r = freqdata[2*(ncfft - k)-1];
+      fnkc.i = -freqdata[2*(ncfft - k)];
+        /*C_FIXDIV( fk , 2 );
+      C_FIXDIV( fnkc , 2 );*/
+
+      C_ADD (fek, fk, fnkc);
+      C_SUB (tmp, fk, fnkc);
+      C_MUL (fok, tmp, st->super_twiddles[k]);
+      C_ADD (st->tmpbuf[k],     fek, fok);
+      C_SUB (st->tmpbuf[ncfft - k], fek, fok);
+#ifdef USE_SIMD        
+      st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0);
+#else
+      st->tmpbuf[ncfft - k].i *= -1;
+#endif
+   }
+   kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata);
+}

Modified: trunk/speex/libspeex/kiss_fftr.h
===================================================================
--- trunk/speex/libspeex/kiss_fftr.h	2007-03-05 09:14:02 UTC (rev 12644)
+++ trunk/speex/libspeex/kiss_fftr.h	2007-03-05 12:32:16 UTC (rev 12645)
@@ -32,6 +32,8 @@
  output freqdata has nfft/2+1 complex points
 */
 
+void kiss_fftr2(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_scalar *freqdata);
+
 void kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata);
 /*
  input freqdata has  nfft/2+1 complex points



More information about the commits mailing list