Finite Impulse Filter and decimation test program

Dependencies:   mbed

Revision:
0:3425cd32678e
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main.cpp	Fri Mar 26 12:53:06 2010 +0000
@@ -0,0 +1,127 @@
+
+// fir_int.cpp  Tm Owen  25/3/2010
+
+// FIR decimation filter implementation using integers on 3 channel data.
+// Test/timing program for implementing a simple FIR filter 
+// This is a low pass filter and decimator implemented on the MBED
+// as an integer calculation in native 32 bit precision.
+// The normal filter coefficcients ( eg, see program SCOPEFIR at XXXX.com - free use for 60 days) are
+// given as floats that add up to about 1.  For an integer implementation
+// the coefficients must be multiplied by 2^N to give reasonable integers.
+// Given 32 bit arithmetic, we can handle data and coefficients as follows
+// 32 = COEF  +  DATA + FILTER_LENGTH   all expressed as bits 
+// We are dealing with 13 bit data and a 24 coef long filter ( 5 bit )
+// so we can have coefficients up to 14 bits  - data and coef resolution is 
+// balanced. Anything bigger will overflow 32 bits somewhere.
+// SO;- multiply the coefs by 2^14 ( = 16384 ), do the sums and divide the
+// answer by 2^14 or shift it right 14 places.
+// I put in the multiplied coefs, could be done once in software. 
+// DECIM is the decimation ratio - the number of input samples that will be 
+// consumed for each output sample.  Normally this is a power of 2, but
+// it can be any integer?  You need to design the coefficients based on the
+// sampling rate and output rate - the filter  cutoff frequency must be less than
+// half the output sample frequency.
+
+// I get this to run at about 8 microsecs per  3 output values - there are 72 int mults and around 125 adds.
+// in the for loop - quite fast enough for most data logging functions but a bit slow for audio. 
+
+// input data is in a .csv file, data is output to a similar file.
+
+#include "mbed.h"
+Serial pc(USBTX, USBRX);                // tx, rx
+DigitalOut led(LED1);                 // useful indicator of progress
+FILE *fp, *fp1;
+Timer t;                               // output write file - bin
+LocalFileSystem local("local");
+#define SHIFT       15    // coefficient multiplier as a power of 2   
+#define COEF        24    // number of coefficients in fir
+#define DECIM        3    // decimation factor required
+#define BUFSIZE    400    // number of input samples to process 
+
+int xarr[COEF], yarr[COEF], zarr[COEF];
+
+// The FIR filter values shown below are good for 13 bit data at an oputput cutoff of 1/8 input rate
+// and a -80db cutoff at 1/4 of input rate, e.g for 800 s/s gives 200 s/s with 100Hz -3db point.
+int coef[] = {-4,
+33,
+208,
+528,
+720,
+3121,
+-802,
+-1784,
+-1113,
+2081,
+6756,
+10245,
+10245,
+6756,
+2081,
+-1113,
+-1784,
+-802,
+312,
+720,
+528,
+208,
+33,
+-4
+ };
+int x,y,z, count =0;
+int xout=0,yout=0,zout=0;
+int once =1;
+int bufx[BUFSIZE],bufy[BUFSIZE],bufz[BUFSIZE],buf_sump =0,buf_curp;
+int obufx[BUFSIZE/DECIM],obufy[BUFSIZE/DECIM],obufz[BUFSIZE/DECIM];
+int elapse_us;
+main()
+    {
+    int ind;
+    pc.printf("\n\rFIR_INT is an integer implementation of a low pass FIR decimation filter");
+    pc.printf("\n\r...................... Tim Owen .... March 2010..........................\n\r\n\r");
+    pc.printf("\n\r Filter length = %d,  Decimation factor = %d, Shift factor = 2^%d\n\r", COEF, DECIM,SHIFT);  
+    pc.printf("\n\rOpen two files - one in and one out.........");
+    fp = fopen("/local/TXT2.TXT", "r")  ;
+    if(!fp) pc.printf("\n\rNo go input file");
+    fp1 = fopen("/local/TXTFILT.TXT","w");
+     if(!fp1) pc.printf("\n\rNo go output file");
+     pc.printf("\n\r coefficients %d %d %d %d %d %d %d %d %d %d %d %d", coef[0],coef[1],coef[2],coef[3],coef[4], coef[5], coef[6],coef[7],coef[8],coef[9],coef[10],coef[11]);
+     t.start();
+      
+      for(ind = 0; ind < BUFSIZE ; ind++)
+        { // copy data to array so we can time without reads
+           if( fscanf(fp, "%d,%d,%d",&bufx[ind],&bufy[ind],&bufz[ind]) != EOF){};   // fscanf > ptrs
+        }
+    t.reset();  // start timing
+        
+    while((count < BUFSIZE/DECIM ) && (feof(fp) == 0))
+        { // go round this loop once for each OUTPUT  value
+            
+             obufx[count] =0; // zero the OUTPUT value
+             obufy[count] =0; 
+             obufz[count] =0;   
+             
+             buf_sump = buf_curp;  // set sumation index  from current start position    
+             
+          for( ind = 0; ind < COEF; ind++)
+            {// do the summation
+               
+               obufx[count]  +=(bufx[buf_sump] * coef[ind]);
+               obufy[count]  +=(bufy[buf_sump] * coef[ind]);
+               obufz[count]  +=(bufz[buf_sump] * coef[ind]);
+               buf_sump++; 
+            }
+              buf_curp += DECIM ;  // set up start for next sumation 
+             count++;
+         }     
+          elapse_us = t.read_us();  
+         for( ind =0; ind < BUFSIZE/DECIM ; ind++)
+            {          
+             fprintf(fp1,"%d,%d,%d\r",(obufx[ind]>>SHIFT),(obufy[ind]>>SHIFT),(obufz[ind]>>SHIFT));
+            } 
+         fclose(fp);
+         fclose(fp1);  
+        pc.printf("\n\rSamples = %d, elapse time = %dus  us/count = %d\n\r BFN..... \n\r", count,elapse_us, elapse_us/count);
+    }
+    
+    
+    
\ No newline at end of file