tim owen
/
fir_int
Finite Impulse Filter and decimation test program
Revision 0:3425cd32678e, committed 2010-03-26
- 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