From 9b738dc4928fbf0b985ab811b8b9157f60066d43 Mon Sep 17 00:00:00 2001
From: Mark Borgerding <mark@borgerding.net>
Date: Fri, 30 Jan 2004 00:17:01 +0000
Subject: [PATCH] slight fast fir changes

---
 README               |  5 ++-
 _kiss_fft_guts.h     |  2 +-
 kiss_fft.c           |  1 -
 tools/kiss_fastfir.c | 81 ++++++++++++++++++++++++++++++++++++++++++--
 4 files changed, 82 insertions(+), 7 deletions(-)

diff --git a/README b/README
index 2660033..fc881c9 100644
--- a/README
+++ b/README
@@ -72,9 +72,8 @@ LICENSE:
     BSD, see COPYING for details. Basically, "free to use, give credit where due, no guarantees"
 
 TODO:
-    *) Add 2d real optimized FFT
-    *) Make a better self-test program(s).  Should report snr & timing for short,float, 
-        or double.
+    *) Add real optimization for odd length FFTs
+    *) Add real optimization to the n-dimensional FFT
     *) Add simple windowing function, e.g. Hamming : w(i)=.54-.46*cos(2pi*i/(n-1))
     *) Make the fixed point scaling and bit shifts more easily configurable.
     *) Document/revisit the input/output fft scaling
diff --git a/_kiss_fft_guts.h b/_kiss_fft_guts.h
index b17c159..9071c06 100644
--- a/_kiss_fft_guts.h
+++ b/_kiss_fft_guts.h
@@ -70,7 +70,7 @@ typedef struct {
 #define C_SUBFROM( res , a)\
     do {    (res).r -= (a).r;  (res).i -= (a).i;  }while(0)
 
-static
+static inline
 kiss_fft_cpx kf_cexp(double phase) /* returns e ** (j*phase)   */
 {
     kiss_fft_cpx x;
diff --git a/kiss_fft.c b/kiss_fft.c
index 669e605..8a2d330 100644
--- a/kiss_fft.c
+++ b/kiss_fft.c
@@ -297,7 +297,6 @@ void kf_factor(int n,int * facbuf)
 }
 
 /*
- *      void * kiss_fft_alloc(int nfft,int inverse_fft)
  *
  * User-callable function to allocate all necessary scratch space for the fft.
  *
diff --git a/tools/kiss_fastfir.c b/tools/kiss_fastfir.c
index 4c97dbb..62d9425 100644
--- a/tools/kiss_fastfir.c
+++ b/tools/kiss_fastfir.c
@@ -231,6 +231,73 @@ size_t kiss_fastfir(
 #include <unistd.h>
 #include <sys/types.h>
 #include <sys/mman.h>
+#include <assert.h>
+
+void direct_file_filter(
+        FILE * fin,
+        FILE * fout,
+        const kffsamp_t * imp_resp,
+        size_t n_imp_resp,
+        size_t nfft )
+{
+    size_t nlag = n_imp_resp - 1;
+
+    const kffsamp_t *tmph;
+    kffsamp_t *buf, *lagbuf;
+    kffsamp_t outval;
+    size_t nread;
+    size_t nbuf;
+    size_t oldestlag = 0;
+    size_t i, ii;
+
+    nbuf = 4096;
+    buf = (kffsamp_t *) malloc ( sizeof (kffsamp_t) * nbuf);
+    lagbuf = (kffsamp_t *) malloc (sizeof (kffsamp_t) * nlag);
+    if (!lagbuf || !buf) {
+        perror("lagbuf allocation");
+        exit(1);
+    }
+
+    if ( fread (lagbuf, sizeof (kffsamp_t), nlag, fin) !=  nlag ) {
+        perror ("insufficient data to overcome transient");
+        exit (1);
+    }
+
+    do {
+        nread = fread (buf, sizeof (kffsamp_t), nbuf, fin);
+        if (nread <= 0)
+            break;
+
+        for (i = 0; i < nread; ++i) {
+            outval = 0;
+            tmph = imp_resp;
+
+
+            for (ii = oldestlag; ii < nlag; ++ii)
+                outval += lagbuf[ii] * *tmph++;
+
+
+            for (ii = 0; ii < oldestlag; ++ii)
+                outval += lagbuf[ii] * *tmph++;
+
+            outval += buf[i] * *tmph++;
+
+
+            lagbuf[oldestlag] = buf[i];
+            buf[i] = outval;
+
+            if (++oldestlag == nlag)
+                oldestlag = 0;
+        }
+
+        if (fwrite (buf, sizeof (buf[0]), nread, fout) != nread) {
+            perror ("short write");
+            exit (1);
+        }
+    } while (nread);
+    free (buf);
+    free (lagbuf);
+}
 
 void do_file_filter(
         FILE * fin,
@@ -283,12 +350,13 @@ void do_file_filter(
 int main(int argc,char**argv)
 {
     kffsamp_t * h;
+    int use_direct=0;
     size_t nh,nfft=0;
     FILE *fin=stdin;
     FILE *fout=stdout;
     FILE *filtfile=NULL;
     while (1) {
-        int c=getopt(argc,argv,"n:h:i:o:v");
+        int c=getopt(argc,argv,"n:h:i:o:vd");
         if (c==-1) break;
         switch (c) {
             case 'v':
@@ -318,10 +386,16 @@ int main(int argc,char**argv)
                     exit(1);
                 }
                 break;
+            case 'd':
+                use_direct=1;
+                break;
             case '?':
                      fprintf(stderr,"usage options:\n"
+                            "\t-n nfft: fft size to use\n"
+                            "\t-d : use direct FIR filtering, not fast convolution\n"
                             "\t-i filename: input file\n"
                             "\t-o filename: output(filtered) file\n"
+                            "\t-n nfft: fft size to use\n"
                             "\t-h filename: impulse response\n");
                      exit (1);
             default:fprintf(stderr,"bad %c\n",c);break;
@@ -339,7 +413,10 @@ int main(int argc,char**argv)
     fread(h,sizeof(kffsamp_t),nh,filtfile);
     fclose(filtfile);
  
-    do_file_filter( fin, fout, h,nh,nfft);
+    if (use_direct)
+        direct_file_filter( fin, fout, h,nh,nfft);
+    else
+        do_file_filter( fin, fout, h,nh,nfft);
 
     if (fout!=stdout) fclose(fout);
     if (fin!=stdin) fclose(fin);