Finite Impulse Filter and decimation test program

Dependencies:   mbed

Files at this revision

API Documentation at this revision

Comitter:
timowen
Date:
Fri Mar 26 12:53:06 2010 +0000
Commit message:

Changed in this revision

main.cpp Show annotated file Show diff for this revision Revisions of this file
mbed.bld Show annotated file Show diff for this revision Revisions of this file
diff -r 000000000000 -r 3425cd32678e main.cpp
--- /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
diff -r 000000000000 -r 3425cd32678e mbed.bld
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mbed.bld	Fri Mar 26 12:53:06 2010 +0000
@@ -0,0 +1,1 @@
+http://mbed.org/users/mbed_official/code/mbed/builds/49a220cc26e0