tim owen
Finite Impulse Filter and decimation test program
main.cpp@0:3425cd32678e, 2010-03-26 (annotated)
- Committer:
- timowen
- Date:
- Fri Mar 26 12:53:06 2010 +0000
- Revision:
- 0:3425cd32678e
Who changed what in which revision?
| User | Revision | Line number | New contents of line |
|---|---|---|---|
| timowen | 0:3425cd32678e | 1 | |
| timowen | 0:3425cd32678e | 2 | // fir_int.cpp Tm Owen 25/3/2010 |
| timowen | 0:3425cd32678e | 3 | |
| timowen | 0:3425cd32678e | 4 | // FIR decimation filter implementation using integers on 3 channel data. |
| timowen | 0:3425cd32678e | 5 | // Test/timing program for implementing a simple FIR filter |
| timowen | 0:3425cd32678e | 6 | // This is a low pass filter and decimator implemented on the MBED |
| timowen | 0:3425cd32678e | 7 | // as an integer calculation in native 32 bit precision. |
| timowen | 0:3425cd32678e | 8 | // The normal filter coefficcients ( eg, see program SCOPEFIR at XXXX.com - free use for 60 days) are |
| timowen | 0:3425cd32678e | 9 | // given as floats that add up to about 1. For an integer implementation |
| timowen | 0:3425cd32678e | 10 | // the coefficients must be multiplied by 2^N to give reasonable integers. |
| timowen | 0:3425cd32678e | 11 | // Given 32 bit arithmetic, we can handle data and coefficients as follows |
| timowen | 0:3425cd32678e | 12 | // 32 = COEF + DATA + FILTER_LENGTH all expressed as bits |
| timowen | 0:3425cd32678e | 13 | // We are dealing with 13 bit data and a 24 coef long filter ( 5 bit ) |
| timowen | 0:3425cd32678e | 14 | // so we can have coefficients up to 14 bits - data and coef resolution is |
| timowen | 0:3425cd32678e | 15 | // balanced. Anything bigger will overflow 32 bits somewhere. |
| timowen | 0:3425cd32678e | 16 | // SO;- multiply the coefs by 2^14 ( = 16384 ), do the sums and divide the |
| timowen | 0:3425cd32678e | 17 | // answer by 2^14 or shift it right 14 places. |
| timowen | 0:3425cd32678e | 18 | // I put in the multiplied coefs, could be done once in software. |
| timowen | 0:3425cd32678e | 19 | // DECIM is the decimation ratio - the number of input samples that will be |
| timowen | 0:3425cd32678e | 20 | // consumed for each output sample. Normally this is a power of 2, but |
| timowen | 0:3425cd32678e | 21 | // it can be any integer? You need to design the coefficients based on the |
| timowen | 0:3425cd32678e | 22 | // sampling rate and output rate - the filter cutoff frequency must be less than |
| timowen | 0:3425cd32678e | 23 | // half the output sample frequency. |
| timowen | 0:3425cd32678e | 24 | |
| timowen | 0:3425cd32678e | 25 | // I get this to run at about 8 microsecs per 3 output values - there are 72 int mults and around 125 adds. |
| timowen | 0:3425cd32678e | 26 | // in the for loop - quite fast enough for most data logging functions but a bit slow for audio. |
| timowen | 0:3425cd32678e | 27 | |
| timowen | 0:3425cd32678e | 28 | // input data is in a .csv file, data is output to a similar file. |
| timowen | 0:3425cd32678e | 29 | |
| timowen | 0:3425cd32678e | 30 | #include "mbed.h" |
| timowen | 0:3425cd32678e | 31 | Serial pc(USBTX, USBRX); // tx, rx |
| timowen | 0:3425cd32678e | 32 | DigitalOut led(LED1); // useful indicator of progress |
| timowen | 0:3425cd32678e | 33 | FILE *fp, *fp1; |
| timowen | 0:3425cd32678e | 34 | Timer t; // output write file - bin |
| timowen | 0:3425cd32678e | 35 | LocalFileSystem local("local"); |
| timowen | 0:3425cd32678e | 36 | #define SHIFT 15 // coefficient multiplier as a power of 2 |
| timowen | 0:3425cd32678e | 37 | #define COEF 24 // number of coefficients in fir |
| timowen | 0:3425cd32678e | 38 | #define DECIM 3 // decimation factor required |
| timowen | 0:3425cd32678e | 39 | #define BUFSIZE 400 // number of input samples to process |
| timowen | 0:3425cd32678e | 40 | |
| timowen | 0:3425cd32678e | 41 | int xarr[COEF], yarr[COEF], zarr[COEF]; |
| timowen | 0:3425cd32678e | 42 | |
| timowen | 0:3425cd32678e | 43 | // The FIR filter values shown below are good for 13 bit data at an oputput cutoff of 1/8 input rate |
| timowen | 0:3425cd32678e | 44 | // and a -80db cutoff at 1/4 of input rate, e.g for 800 s/s gives 200 s/s with 100Hz -3db point. |
| timowen | 0:3425cd32678e | 45 | int coef[] = {-4, |
| timowen | 0:3425cd32678e | 46 | 33, |
| timowen | 0:3425cd32678e | 47 | 208, |
| timowen | 0:3425cd32678e | 48 | 528, |
| timowen | 0:3425cd32678e | 49 | 720, |
| timowen | 0:3425cd32678e | 50 | 3121, |
| timowen | 0:3425cd32678e | 51 | -802, |
| timowen | 0:3425cd32678e | 52 | -1784, |
| timowen | 0:3425cd32678e | 53 | -1113, |
| timowen | 0:3425cd32678e | 54 | 2081, |
| timowen | 0:3425cd32678e | 55 | 6756, |
| timowen | 0:3425cd32678e | 56 | 10245, |
| timowen | 0:3425cd32678e | 57 | 10245, |
| timowen | 0:3425cd32678e | 58 | 6756, |
| timowen | 0:3425cd32678e | 59 | 2081, |
| timowen | 0:3425cd32678e | 60 | -1113, |
| timowen | 0:3425cd32678e | 61 | -1784, |
| timowen | 0:3425cd32678e | 62 | -802, |
| timowen | 0:3425cd32678e | 63 | 312, |
| timowen | 0:3425cd32678e | 64 | 720, |
| timowen | 0:3425cd32678e | 65 | 528, |
| timowen | 0:3425cd32678e | 66 | 208, |
| timowen | 0:3425cd32678e | 67 | 33, |
| timowen | 0:3425cd32678e | 68 | -4 |
| timowen | 0:3425cd32678e | 69 | }; |
| timowen | 0:3425cd32678e | 70 | int x,y,z, count =0; |
| timowen | 0:3425cd32678e | 71 | int xout=0,yout=0,zout=0; |
| timowen | 0:3425cd32678e | 72 | int once =1; |
| timowen | 0:3425cd32678e | 73 | int bufx[BUFSIZE],bufy[BUFSIZE],bufz[BUFSIZE],buf_sump =0,buf_curp; |
| timowen | 0:3425cd32678e | 74 | int obufx[BUFSIZE/DECIM],obufy[BUFSIZE/DECIM],obufz[BUFSIZE/DECIM]; |
| timowen | 0:3425cd32678e | 75 | int elapse_us; |
| timowen | 0:3425cd32678e | 76 | main() |
| timowen | 0:3425cd32678e | 77 | { |
| timowen | 0:3425cd32678e | 78 | int ind; |
| timowen | 0:3425cd32678e | 79 | pc.printf("\n\rFIR_INT is an integer implementation of a low pass FIR decimation filter"); |
| timowen | 0:3425cd32678e | 80 | pc.printf("\n\r...................... Tim Owen .... March 2010..........................\n\r\n\r"); |
| timowen | 0:3425cd32678e | 81 | pc.printf("\n\r Filter length = %d, Decimation factor = %d, Shift factor = 2^%d\n\r", COEF, DECIM,SHIFT); |
| timowen | 0:3425cd32678e | 82 | pc.printf("\n\rOpen two files - one in and one out........."); |
| timowen | 0:3425cd32678e | 83 | fp = fopen("/local/TXT2.TXT", "r") ; |
| timowen | 0:3425cd32678e | 84 | if(!fp) pc.printf("\n\rNo go input file"); |
| timowen | 0:3425cd32678e | 85 | fp1 = fopen("/local/TXTFILT.TXT","w"); |
| timowen | 0:3425cd32678e | 86 | if(!fp1) pc.printf("\n\rNo go output file"); |
| timowen | 0:3425cd32678e | 87 | 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]); |
| timowen | 0:3425cd32678e | 88 | t.start(); |
| timowen | 0:3425cd32678e | 89 | |
| timowen | 0:3425cd32678e | 90 | for(ind = 0; ind < BUFSIZE ; ind++) |
| timowen | 0:3425cd32678e | 91 | { // copy data to array so we can time without reads |
| timowen | 0:3425cd32678e | 92 | if( fscanf(fp, "%d,%d,%d",&bufx[ind],&bufy[ind],&bufz[ind]) != EOF){}; // fscanf > ptrs |
| timowen | 0:3425cd32678e | 93 | } |
| timowen | 0:3425cd32678e | 94 | t.reset(); // start timing |
| timowen | 0:3425cd32678e | 95 | |
| timowen | 0:3425cd32678e | 96 | while((count < BUFSIZE/DECIM ) && (feof(fp) == 0)) |
| timowen | 0:3425cd32678e | 97 | { // go round this loop once for each OUTPUT value |
| timowen | 0:3425cd32678e | 98 | |
| timowen | 0:3425cd32678e | 99 | obufx[count] =0; // zero the OUTPUT value |
| timowen | 0:3425cd32678e | 100 | obufy[count] =0; |
| timowen | 0:3425cd32678e | 101 | obufz[count] =0; |
| timowen | 0:3425cd32678e | 102 | |
| timowen | 0:3425cd32678e | 103 | buf_sump = buf_curp; // set sumation index from current start position |
| timowen | 0:3425cd32678e | 104 | |
| timowen | 0:3425cd32678e | 105 | for( ind = 0; ind < COEF; ind++) |
| timowen | 0:3425cd32678e | 106 | {// do the summation |
| timowen | 0:3425cd32678e | 107 | |
| timowen | 0:3425cd32678e | 108 | obufx[count] +=(bufx[buf_sump] * coef[ind]); |
| timowen | 0:3425cd32678e | 109 | obufy[count] +=(bufy[buf_sump] * coef[ind]); |
| timowen | 0:3425cd32678e | 110 | obufz[count] +=(bufz[buf_sump] * coef[ind]); |
| timowen | 0:3425cd32678e | 111 | buf_sump++; |
| timowen | 0:3425cd32678e | 112 | } |
| timowen | 0:3425cd32678e | 113 | buf_curp += DECIM ; // set up start for next sumation |
| timowen | 0:3425cd32678e | 114 | count++; |
| timowen | 0:3425cd32678e | 115 | } |
| timowen | 0:3425cd32678e | 116 | elapse_us = t.read_us(); |
| timowen | 0:3425cd32678e | 117 | for( ind =0; ind < BUFSIZE/DECIM ; ind++) |
| timowen | 0:3425cd32678e | 118 | { |
| timowen | 0:3425cd32678e | 119 | fprintf(fp1,"%d,%d,%d\r",(obufx[ind]>>SHIFT),(obufy[ind]>>SHIFT),(obufz[ind]>>SHIFT)); |
| timowen | 0:3425cd32678e | 120 | } |
| timowen | 0:3425cd32678e | 121 | fclose(fp); |
| timowen | 0:3425cd32678e | 122 | fclose(fp1); |
| timowen | 0:3425cd32678e | 123 | pc.printf("\n\rSamples = %d, elapse time = %dus us/count = %d\n\r BFN..... \n\r", count,elapse_us, elapse_us/count); |
| timowen | 0:3425cd32678e | 124 | } |
| timowen | 0:3425cd32678e | 125 | |
| timowen | 0:3425cd32678e | 126 | |
| timowen | 0:3425cd32678e | 127 |