# Shawn's DSP Tutorials

tutorials on practical aspects of digital signal processing

## Implementation of FIR Filtering in C (Part 2)

In Part 1 I showed how to code a FIR filter in C using floating point. In this lesson I will show how to do the same thing using fixed point operations. The code example below will demonstrate the application of fixed point multiplication, rounding and saturation. The code has definitions for the FIR filtering function, followed by an example test program.

The following link is a PDF version of the code example:

firFixed

And here is the code example:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 | `#include <stdio.h>` `#include <stdint.h>` `#include <stdlib.h>` `#include <string.h>` `//////////////////////////////////////////////////////////////` `// Filter Code Definitions` `//////////////////////////////////////////////////////////////` `// maximum number of inputs that can be handled` `// in one function call` `#define MAX_INPUT_LEN 80` `// maximum length of filter than can be handled` `#define MAX_FLT_LEN 63` `// buffer to hold all of the input samples` `#define BUFFER_LEN (MAX_FLT_LEN - 1 + MAX_INPUT_LEN)` `// array to hold input samples` `int16_t insamp[ BUFFER_LEN ];` `// FIR init` `void` `firFixedInit(` `void` `)` `{` `` `memset``( insamp, 0,` `sizeof``( insamp ) );` `}` `// the FIR filter function` `void` `firFixed( int16_t *coeffs, int16_t *input, int16_t *output,` `` `int` `length,` `int` `filterLength )` `{` `` `int32_t acc;` `// accumulator for MACs` `` `int16_t *coeffp;` `// pointer to coefficients` `` `int16_t *inputp;` `// pointer to input samples` `` `int` `n;` `` `int` `k;` `` `// put the new samples at the high end of the buffer` `` `memcpy``( &insamp[filterLength - 1], input,` `` `length *` `sizeof``(int16_t) );` `` `// apply the filter to each input sample` `` `for` `( n = 0; n < length; n++ ) {` `` `// calculate output n` `` `coeffp = coeffs;` `` `inputp = &insamp[filterLength - 1 + n];` `` `// load rounding constant` `` `acc = 1 << 14;` `` `// perform the multiply-accumulate` `` `for` `( k = 0; k < filterLength; k++ ) {` `` `acc += (int32_t)(*coeffp++) * (int32_t)(*inputp--);` `` `}` `` `// saturate the result` `` `if` `( acc > 0x3fffffff ) {` `` `acc = 0x3fffffff;` `` `}` `else` `if` `( acc < -0x40000000 ) {` `` `acc = -0x40000000;` `` `}` `` `// convert from Q30 to Q15` `` `output[n] = (int16_t)(acc >> 15);` `` `}` `` `// shift input samples back in time for next time` `` `memmove``( &insamp[0], &insamp[length],` `` `(filterLength - 1) *` `sizeof``(int16_t) );` `}` `//////////////////////////////////////////////////////////////` `// Test program` `//////////////////////////////////////////////////////////////` `// bandpass filter centred around 1000 Hz` `// sampling rate = 8000 Hz` `// gain at 1000 Hz is about 1.13` `#define FILTER_LEN 63` `int16_t coeffs[ FILTER_LEN ] =` `{` `` `-1468, 1058, 594, 287, 186, 284, 485, 613,` `` `495, 90, -435, -762, -615, 21, 821, 1269,` `` `982, 9, -1132, -1721, -1296, 1, 1445, 2136,` `` `1570, 0, -1666, -2413, -1735, -2, 1770, 2512,` `` `1770, -2, -1735, -2413, -1666, 0, 1570, 2136,` `` `1445, 1, -1296, -1721, -1132, 9, 982, 1269,` `` `821, 21, -615, -762, -435, 90, 495, 613,` `` `485, 284, 186, 287, 594, 1058, -1468` `};` `// number of samples to read per loop` `#define SAMPLES 80` `int` `main(` `void` `)` `{` `` `int` `size;` `` `int16_t input[SAMPLES];` `` `int16_t output[SAMPLES];` `` `FILE` `*in_fid;` `` `FILE` `*out_fid;` `` `// open the input waveform file` `` `in_fid =` `fopen``(` `"input.pcm"``,` `"rb"` `);` `` `if` `( in_fid == 0 ) {` `` `printf``(``"couldn't open input.pcm"``);` `` `exit``(EXIT_FAILURE);` `` `}` `` `// open the output waveform file` `` `out_fid =` `fopen``(` `"outputFixed.pcm"``,` `"wb"` `);` `` `if` `( out_fid == 0 ) {` `` `printf``(``"couldn't open outputFixed.pcm"``);` `` `exit``(EXIT_FAILURE);` `` `}` `` `// initialize the filter` `` `firFixedInit();` `` `// process all of the samples` `` `do` `{` `` `// read samples from file` `` `size =` `fread``( input,` `sizeof``(int16_t), SAMPLES, in_fid );` `` `// perform the filtering` `` `firFixed( coeffs, input, output, size, FILTER_LEN );` `` `// write samples to file` `` `fwrite``( output,` `sizeof``(int16_t), size, out_fid );` `` `}` `while` `( size != 0 );` `` `fclose``( in_fid );` `` `fclose``( out_fid );` `` `return` `0;` `}`

The first thing to notice is that the definitions for the input sample storage and handling are nearly the same as for the code in Part 1. The only difference is that the storage type is a 16 bit integer instead of double precision floating point.

The next difference is the inclusion of rounding in the calculation of each output. Rounding is used when converting the calculated result from a Q30 format number to Q15. Notice that I have loaded the rounding constant into the accumulator value (acc) at the beginning of the loop rather than adding it at the end. This is a small optimization commonly seen in code for FIR filters. If you are coding in assembly language, and your chip has a rounding instruction, it may be better to do the rounding at the end (depending on what the instruction actually does).

The multiplication itself is now an integer multiplication of two 16 bit values, each of which is a Q15 number. The accumulator variable is 32 bits, and holds a Q1.30 format number. There is one bit for the sign, one integer bit, and thirty fractional bits. Notice that I have cast each multiplier to a 32 bit value. Failure to do so should result in a 16 bit product and produce wrong results.

Next comes the overflow handling for converting from Q1.30 to Q30. The code checks for values beyond the limits of the largest/smallest Q30 number (no integer bits), and saturates if necessary.

Finally, a right shift by 15 is used to convert the MAC result from Q30 to Q15 and the result is stored to the output array.

The test program is simpler than the one in Part 1 because there is no need to convert the input and output samples to or from floating point. The most important thing to note is the change in the filter coefficient array. To generate these coefficients, I took the floating point coefficients from Part 1, multiplied by 32768, and then rounded to the nearest integer. The coefficients are in Q15 format, and note that none of the original floating point coefficients are close to one. Multiplying by 32768 would cause a problem for any coefficients larger than 32767/32768 or less than -1.

As in Part 1, the test input file should be 16 bit samples at a sampling rate of 8000 Hz. Try using a 1000 Hz tone with noise added to it. The bandpass filter should remove a large portion of the noise.

The filter has greater than unity gain at 1000 Hz and the gain is about 1.13. Admittedly it is not a great filter design, but it tests the saturation code if a full scale 1000 Hz sine wave is used as an input. Try a full scale 1000 Hz tone input with and without the saturation check and see what the difference is. There should be a small amount of distortion with the saturation code present (due to slight flattening of tone) and a large amount of distortion without it (due to overflow).

A fixed point FIR filter is fairly easy to implement. Personally I find the sample storage and movement to be the most difficult part. In Part 3 of this lesson I will demonstrate how to run multiple FIR filters on the same input.

### Like this:

Like Loading...

### 6 Comments on “Implementation of FIR Filtering in C (Part 2)”

1. alonso Says:
   June 9, 2013 at 6:20 am

   Hi!
   Thank you for the very helpfull website.

   I have a quation. i get the input data as unsigned int16 from (16bit ADC). this input should be filtered with a dsp. the coeff of the filter are in fixed-point ( Q16.15). should a transformation take place between the input int16 and fixed point .

   Thanks again

   Reply
   - Shawn Says:
     June 9, 2013 at 8:13 pm

     Hi Alonso,

     I’m glad you find the website helpful. The ADC input may need to be transformed. If the ADC has a mode where it produces two’s complement numbers, then no transformation is needed. If the ADC is using a different number representation (such as signed magnitude) than a transformation would be required. Take a look at the documentation for the ADC to see what format the sampled values will have.

     Hope that helps.

     Reply
2. alonso Says:
   June 10, 2013 at 10:04 am

   Hi, Thanks for the answer.
   The ADC delivers unsigned int16. the DSP works with fixed point Q16.15.

   what about overflow? i read somewhe that the coeff should be scaled? do you have any idea about avoiding overflow ?

   sincerly

   Reply
   - Shawn Says:
     June 12, 2013 at 8:45 pm

     To avoid overflow, you can scale the coefficients so that the largest output is always less than one for any possible input. So when designing the filter, ensure the magnitude response is slightly less than one over the whole frequency range. Do this before converting the coefficients from floating point to fixed point.

     Also, I’ve written the code so that the result will saturate at either extreme rather than overflowing. But if you scale your coefficients, you can save some cycles by removing the overflow checks.

     Reply
3. Ben Morse (@Benathon) Says:
   September 1, 2015 at 12:49 pm

   You need to add #include string.h to both examples or else it wont compile with gcc. also you should “return 1;” on lines 102,109

   Reply
   - Shawn Says:
     September 1, 2015 at 10:25 pm

     Thanks Ben. I did compile and test the code with GCC, but that was many years ago. I added string.h and stdlib.h, and changed the error returns to “exit” statements. Let me know if you see any more problems. I had some issues with the formatting concerning , and & characters. I haven’t updated the pdf example, and may just delete the link. I should check/fix my examples on other pages as well…

     Reply

### Leave a Reply Cancel reply

Enter your comment here...

Fill in your details below or click an icon to log in:

Email (required) (Address never made public)
Name (required)
Website

You are commenting using your WordPress.com account. ( Log Out / Change )

You are commenting using your Twitter account. ( Log Out / Change )

You are commenting using your Facebook account. ( Log Out / Change )

You are commenting using your Google+ account. ( Log Out / Change )

Cancel

Connecting to %s

Notify me of new comments via email.

Follow

### Follow “Shawn's DSP Tutorials”

Get every new post delivered to your Inbox.

Join 25 other followers

Build a website with WordPress.com
%d bloggers like this: