# Shawn's DSP Tutorials

tutorials on practical aspects of digital signal processing

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

In this lesson I will show how to code a finite impulse response (FIR) digital filter in the C programming language using floating point operations. Understanding the basic FIR filter is a stepping stone to understanding more complicated filters such as polyphase FIR filters and adaptive FIR filters.

A FIR filter is one in which the output depends only on the inputs (and not on previous outputs). The process is a discrete-time convolution between the input signal and the impulse response of the filter to apply. I will call this impulse response “the filter coefficients.” The following equivalent equations describe the convolution process:

equation 1:

equation 2:

where h is the impulse response, x is the input, y[n] is the output at time n and N is the length of the impulse response (the number of filter coefficients). For each output calculated, there are N multiply-accumulate operations (MACs). There is a choice of going through the filter coefficients forward and the input backwards (equation 1), or through the filter coefficients backwards and the input forwards (equation 2). In the example that follows, I will use the form in equation 1.

In my code example, samples will be stored in a working array, with lower indices being further back in time (at the “low” end). The input samples will be arranged in an array long enough to apply all of the filter coefficients. The length of the array required to process one input sample is N (the number of coefficients in the filter). The length required for M input samples is N – 1 + M.

A function will be defined that filters several input samples at a time. For example, a realtime system for processing telephony quality speech might process eighty samples at a time. This is done to improve the efficiency of the processing, since the overhead of calling the function is eighty times smaller than if it were called for every single sample. The drawback is that additional delay is introduced in the output.

The basic steps for applying a FIR filter are the following:

1. Arrange the new samples at the high end of the input sample buffer.
2. Loop through an outer loop that produces each output sample.
3. Loop through an inner loop that multiplies each filter coefficient by an input sample and adds to a running sum.
4. Shift the previous input samples back in time by the number of new samples that were just processed.

The code example defines an initialization function and a filtering function, and includes a test program. Download the PDF from the following link to see the code, or view it in-line below:

firFloat

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 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 | `#include <stdio.h>` `#include <stdint.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` `double` `insamp[ BUFFER_LEN ];` `// FIR init` `void` `firFloatInit(` `void` `)` `{` `` `memset``( insamp, 0,` `sizeof``( insamp ) );` `}` `// the FIR filter function` `void` `firFloat(` `double` `*coeffs,` `double` `*input,` `double` `*output,` `` `int` `length,` `int` `filterLength )` `{` `` `double` `acc;` `// accumulator for MACs` `` `double` `*coeffp;` `// pointer to coefficients` `` `double` `*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``(``double``) );` `` `// apply the filter to each input sample` `` `for` `( n = 0; n < length; n++ ) {` `` `// calculate output n` `` `coeffp = coeffs;` `` `inputp = &insamp[filterLength - 1 + n];` `` `acc = 0;` `` `for` `( k = 0; k < filterLength; k++ ) {` `` `acc += (*coeffp++) * (*inputp--);` `` `}` `` `output[n] = acc;` `` `}` `` `// shift input samples back in time for next time` `` `memmove``( &insamp[0], &insamp[length],` `` `(filterLength - 1) *` `sizeof``(``double``) );` `}` `//////////////////////////////////////////////////////////////` `// Test program` `//////////////////////////////////////////////////////////////` `// bandpass filter centred around 1000 Hz` `// sampling rate = 8000 Hz` `#define FILTER_LEN 63` `double` `coeffs[ FILTER_LEN ] =` `{` `` `-0.0448093, 0.0322875, 0.0181163, 0.0087615, 0.0056797,` `` `0.0086685, 0.0148049, 0.0187190, 0.0151019, 0.0027594,` `` `-0.0132676, -0.0232561, -0.0187804, 0.0006382, 0.0250536,` `` `0.0387214, 0.0299817, 0.0002609, -0.0345546, -0.0525282,` `` `-0.0395620, 0.0000246, 0.0440998, 0.0651867, 0.0479110,` `` `0.0000135, -0.0508558, -0.0736313, -0.0529380, -0.0000709,` `` `0.0540186, 0.0766746, 0.0540186, -0.0000709, -0.0529380,` `` `-0.0736313, -0.0508558, 0.0000135, 0.0479110, 0.0651867,` `` `0.0440998, 0.0000246, -0.0395620, -0.0525282, -0.0345546,` `` `0.0002609, 0.0299817, 0.0387214, 0.0250536, 0.0006382,` `` `-0.0187804, -0.0232561, -0.0132676, 0.0027594, 0.0151019,` `` `0.0187190, 0.0148049, 0.0086685, 0.0056797, 0.0087615,` `` `0.0181163, 0.0322875, -0.0448093` `};` `void` `intToFloat( int16_t *input,` `double` `*output,` `int` `length )` `{` `` `int` `i;` `` `for` `( i = 0; i < length; i++ ) {` `` `output[i] = (``double``)input[i];` `` `}` `}` `void` `floatToInt(` `double` `*input, int16_t *output,` `int` `length )` `{` `` `int` `i;` `` `for` `( i = 0; i < length; i++ ) {` `` `if` `( input[i] > 32767.0 ) {` `` `input[i] = 32767.0;` `` `}` `else` `if` `( input[i] < -32768.0 ) {` `` `input[i] = -32768.0;` `` `}` `` `// convert` `` `output[i] = (int16_t)input[i];` `` `}` `}` `// number of samples to read per loop` `#define SAMPLES 80` `int` `main(` `void` `)` `{` `` `int` `size;` `` `int16_t input[SAMPLES];` `` `int16_t output[SAMPLES];` `` `double` `floatInput[SAMPLES];` `` `double` `floatOutput[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"``);` `` `return``;` `` `}` `` `// open the output waveform file` `` `out_fid =` `fopen``(` `"outputFloat.pcm"``,` `"wb"` `);` `` `if` `( out_fid == 0 ) {` `` `printf``(``"couldn't open outputFloat.pcm"``);` `` `return``;` `` `}` `` `// initialize the filter` `` `firFloatInit();` `` `// process all of the samples` `` `do` `{` `` `// read samples from file` `` `size =` `fread``( input,` `sizeof``(int16_t), SAMPLES, in_fid );` `` `// convert to doubles` `` `intToFloat( input, floatInput, size );` `` `// perform the filtering` `` `firFloat( coeffs, floatInput, floatOutput, size,` `` `FILTER_LEN );` `` `// convert to ints` `` `floatToInt( floatOutput, output, size );` `` `// write samples to file` `` `fwrite``( output,` `sizeof``(int16_t), size, out_fid );` `` `}` `while` `( size != 0 );` `` `fclose``( in_fid );` `` `fclose``( out_fid );` `` `return` `0;` `}`

First is the definition of the working array to hold the input samples. I have defined the array length to handle up to 80 input samples per function call, and a filter length up to 63 coefficients. In other words, the filter function that uses this working array will support between 1 to 80 input samples, and a filter length between 1 and 63.

Next is the initialization function, which simply zeroes the entire input sample array. This insures that all of the inputs that come before the first actual input are set to zero.

Then comes the FIR filtering function itself. The function takes pointers to the filter coefficients, the new input samples, and the output sample array. The “length” argument specifies the number of input samples to process (must be between 0 and 80 inclusive). The “filterLength” argument specifies the number of coefficients in the filter (must be between 1 and 63 inclusive). Note that the same filter coefficients and filterLength should be passed in each time the firFloat function is called, until the entire input is processed.

The first step in the function is the storage of the new samples. Note closely where the samples are placed in the array.

Next comes the outer loop that processes each input sample to produce one output sample. The outer loop sets up pointers to the coefficients and the input samples, and zeroes an accumulator value for the multiply-accumulate operation in the inner loop. It also stores each calculated output. The inner loop simply performs the multiply-accumulate operation.

Finally comes the step to move the input samples back in time. This is probably the most difficult step in the FIR filter and the part that always takes me the longest to work out. Pay close attention to the position of the samples to move, and the number of samples moved. The amount of time shift is equal to the number of input samples just processed. The number of samples to move is one less than the number of coefficients in the filter. Also note the use of the memmove function rather than the memcpy function, since there is a potential for moving overlapping data (whether or not the data overlaps depends on the relative lengths of the input samples and number of coefficients).

That concludes the code example for the FIR filter. In Part 2, I will show how to do the same thing using fixed point operations.

### Like this:

Like Loading...

### 18 Comments on “Implementation of FIR Filtering in C (Part 1)”

1. walterwhite73 Says:
   June 19, 2012 at 10:27 am

   Hi,
   thanks for the tutorial. just a quick question about the moving the input back by 63 samples. why not move it back 80 samples and discard the first 80-63 of the block. are you using overlap add or overlap save method ?

   Reply
2. walterwhite73 Says:
   June 20, 2012 at 2:33 am

   ok i see now. you move from the 80th, 63 samples back to begining of buffer, makes sense now

   Reply
3. Klevis Says:
   October 10, 2012 at 9:13 am

   Hi,
   Thanks it helped me a lot

   Reply
4. valentin Says:
   December 1, 2012 at 1:19 pm

   There is some bullshit in your example:

   for ( i = 0; i 32767.0 ) {
   095 input[i] = 32767.0;
   096 } else if ( input[i] < -32768.0 ) {
   097 input[i] = -32768.0;
   098 }

   I do not understand the loop and if. Might be it is not necessary an all in this example

   Reply
   - Shawn Says:
     December 2, 2012 at 12:06 am

     Some of the code got deleted somehow. It should make more sense now.

     Reply
5. Deepak Bhatia Says:
   August 14, 2013 at 12:54 am

   How do we generate the input file of sine wav with integer values ?

   Reply
   - Shawn Says:
     September 13, 2013 at 7:25 pm

     Sorry for the late reply. There are some comments about generating the input on this page:

     https://sestevenson.wordpress.com/2010/01/04/implementation-of-fir-filtering-in-c-part-3/

     You can use programs such as Adobe Audition or Audacity to generate sine waves with the correct sampling rate, format, etc.

     Reply
6. Alina Says:
   September 13, 2013 at 7:14 am

   How did you find the coefficients?

   Reply
   - Shawn Says:
     September 13, 2013 at 7:30 pm

     You should look at the comments on the following page (June 18, 2010):

     https://sestevenson.wordpress.com/2010/01/04/implementation-of-fir-filtering-in-c-part-3/

     I used a program called Scilab to generate the coefficients. The Scilab code snippet in the comment shows how to generate fixed point coefficients, so if you want floating point, then leave out the multiplication by 32768 and rounding.

     Reply
     - Alina Says:
       September 14, 2013 at 4:57 am

       That was really helpful Thanks a lot!

7. nig Says:
   December 19, 2013 at 5:14 am

   In case of .wav files .. how can I get wav file header information before starting the filtering?

   Reply
   - Shawn Says:
     December 21, 2013 at 9:57 pm

     Hi. Take a look at the comments on this page starting at March 13, 2012:

     https://sestevenson.wordpress.com/2010/01/04/implementation-of-fir-filtering-in-c-part-3/

     The summary is, write code to read the header information (details of which are here:http://www.sonicspot.com/guide/wavefiles.html) or try using libsndfile, which handles reading and writing of audio files.

     I hope that helps.

     Reply
     - nig Says:
       December 22, 2013 at 5:56 am

       Thanks a lot Shawn,
       A quick look at the links provided seem to be very helpful.
       Will be tackling this problem soon, hoping for the best
       Keep up the good work and thanks again

8. WAV file handling and filtering Says:
   December 29, 2013 at 12:09 pm

   […] wav file using C programming. Till now I managed to tackle the filter part following this link: Implementation of FIR Filtering in C (Part 1) | Shawn's DSP Tutorials (I have my own filter co-efficients which have been produced using FDAtool on MATLAB and some other […]

   Reply
9. Andy Says:
   April 15, 2014 at 6:38 pm

   i have input value in txt file, in this program i got output file in txt but its not readable ……
   help please

   Reply
   - Shawn Says:
     April 15, 2014 at 9:30 pm

     The program does not read nor write text files, but uses binary files instead. The format is raw PCM, 16 bits per sample, 8000 Hz sampling rate. Take a look at some of the previous comments on this page. I if ever get around to it, I’ll add some code to read and write .wav files, which would be a bit easier to use than raw PCM.

     Reply
10. hamid Says:
    December 20, 2015 at 5:09 am

    hi dear
    my samplig rate is 1Mhz
    lenght of sampling is 256 sample
    i want to detect 145khz with FIR.
    can i use ur code?
    tanks in advance

    Reply
    - Shawn Says:
      December 21, 2015 at 8:53 pm

      Hamid,

      The code on this page cannot be used to detect a 145 kHz tone. But take a look at my pages about the TKEO starting at:

      https://sestevenson.wordpress.com/tone-detection-using-the-teager-kaiser-energy-operator-part-1/

      I have a couple of pages that describe how to do tone detection using the “energy operator”, which should be preceded by a bandpass filtering of the input signal (the code on this page). You would need to calculate your own filter coefficients with the 145 kHz centre frequency and the 1 MHz sampling frequency.

      -Shawn

      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: