main.cpp

00001 
00002 // fir_int.cpp  Tm Owen  25/3/2010
00003 
00004 // FIR decimation filter implementation using integers on 3 channel data.
00005 // Test/timing program for implementing a simple FIR filter 
00006 // This is a low pass filter and decimator implemented on the MBED
00007 // as an integer calculation in native 32 bit precision.
00008 // The normal filter coefficcients ( eg, see program SCOPEFIR at XXXX.com - free use for 60 days) are
00009 // given as floats that add up to about 1.  For an integer implementation
00010 // the coefficients must be multiplied by 2^N to give reasonable integers.
00011 // Given 32 bit arithmetic, we can handle data and coefficients as follows
00012 // 32 = COEF  +  DATA + FILTER_LENGTH   all expressed as bits 
00013 // We are dealing with 13 bit data and a 24 coef long filter ( 5 bit )
00014 // so we can have coefficients up to 14 bits  - data and coef resolution is 
00015 // balanced. Anything bigger will overflow 32 bits somewhere.
00016 // SO;- multiply the coefs by 2^14 ( = 16384 ), do the sums and divide the
00017 // answer by 2^14 or shift it right 14 places.
00018 // I put in the multiplied coefs, could be done once in software. 
00019 // DECIM is the decimation ratio - the number of input samples that will be 
00020 // consumed for each output sample.  Normally this is a power of 2, but
00021 // it can be any integer?  You need to design the coefficients based on the
00022 // sampling rate and output rate - the filter  cutoff frequency must be less than
00023 // half the output sample frequency.
00024 
00025 // I get this to run at about 8 microsecs per  3 output values - there are 72 int mults and around 125 adds.
00026 // in the for loop - quite fast enough for most data logging functions but a bit slow for audio. 
00027 
00028 // input data is in a .csv file, data is output to a similar file.
00029 
00030 #include "mbed.h"
00031 Serial pc(USBTX, USBRX);                // tx, rx
00032 DigitalOut led(LED1);                 // useful indicator of progress
00033 FILE *fp, *fp1;
00034 Timer t;                               // output write file - bin
00035 LocalFileSystem local("local");
00036 #define SHIFT       15    // coefficient multiplier as a power of 2   
00037 #define COEF        24    // number of coefficients in fir
00038 #define DECIM        3    // decimation factor required
00039 #define BUFSIZE    400    // number of input samples to process 
00040 
00041 int xarr[COEF], yarr[COEF], zarr[COEF];
00042 
00043 // The FIR filter values shown below are good for 13 bit data at an oputput cutoff of 1/8 input rate
00044 // and a -80db cutoff at 1/4 of input rate, e.g for 800 s/s gives 200 s/s with 100Hz -3db point.
00045 int coef[] = {-4,
00046 33,
00047 208,
00048 528,
00049 720,
00050 3121,
00051 -802,
00052 -1784,
00053 -1113,
00054 2081,
00055 6756,
00056 10245,
00057 10245,
00058 6756,
00059 2081,
00060 -1113,
00061 -1784,
00062 -802,
00063 312,
00064 720,
00065 528,
00066 208,
00067 33,
00068 -4
00069  };
00070 int x,y,z, count =0;
00071 int xout=0,yout=0,zout=0;
00072 int once =1;
00073 int bufx[BUFSIZE],bufy[BUFSIZE],bufz[BUFSIZE],buf_sump =0,buf_curp;
00074 int obufx[BUFSIZE/DECIM],obufy[BUFSIZE/DECIM],obufz[BUFSIZE/DECIM];
00075 int elapse_us;
00076 main()
00077     {
00078     int ind;
00079     pc.printf("\n\rFIR_INT is an integer implementation of a low pass FIR decimation filter");
00080     pc.printf("\n\r...................... Tim Owen .... March 2010..........................\n\r\n\r");
00081     pc.printf("\n\r Filter length = %d,  Decimation factor = %d, Shift factor = 2^%d\n\r", COEF, DECIM,SHIFT);  
00082     pc.printf("\n\rOpen two files - one in and one out.........");
00083     fp = fopen("/local/TXT2.TXT", "r")  ;
00084     if(!fp) pc.printf("\n\rNo go input file");
00085     fp1 = fopen("/local/TXTFILT.TXT","w");
00086      if(!fp1) pc.printf("\n\rNo go output file");
00087      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]);
00088      t.start();
00089       
00090       for(ind = 0; ind < BUFSIZE ; ind++)
00091         { // copy data to array so we can time without reads
00092            if( fscanf(fp, "%d,%d,%d",&bufx[ind],&bufy[ind],&bufz[ind]) != EOF){};   // fscanf > ptrs
00093         }
00094     t.reset();  // start timing
00095         
00096     while((count < BUFSIZE/DECIM ) && (feof(fp) == 0))
00097         { // go round this loop once for each OUTPUT  value
00098             
00099              obufx[count] =0; // zero the OUTPUT value
00100              obufy[count] =0; 
00101              obufz[count] =0;   
00102              
00103              buf_sump = buf_curp;  // set sumation index  from current start position    
00104              
00105           for( ind = 0; ind < COEF; ind++)
00106             {// do the summation
00107                
00108                obufx[count]  +=(bufx[buf_sump] * coef[ind]);
00109                obufy[count]  +=(bufy[buf_sump] * coef[ind]);
00110                obufz[count]  +=(bufz[buf_sump] * coef[ind]);
00111                buf_sump++; 
00112             }
00113               buf_curp += DECIM ;  // set up start for next sumation 
00114              count++;
00115          }     
00116           elapse_us = t.read_us();  
00117          for( ind =0; ind < BUFSIZE/DECIM ; ind++)
00118             {          
00119              fprintf(fp1,"%d,%d,%d\r",(obufx[ind]>>SHIFT),(obufy[ind]>>SHIFT),(obufz[ind]>>SHIFT));
00120             } 
00121          fclose(fp);
00122          fclose(fp1);  
00123         pc.printf("\n\rSamples = %d, elapse time = %dus  us/count = %d\n\r BFN..... \n\r", count,elapse_us, elapse_us/count);
00124     }
00125     
00126     
00127