Finite Impulse Filter and decimation test program

Dependencies:   mbed

main.cpp

Committer:
timowen
Date:
2010-03-26
Revision:
0:3425cd32678e

File content as of revision 0:3425cd32678e:


// 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);
    }