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: y(n)=\sum_{k=0}^{N-1}h(k)x(n-k)

equation 2: y(n)=\sum_{k=n-(N-1)}^{n}h(n-k)x(k)

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

#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  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.

Advertisements
Explore posts in the same categories: FIR filters

Tags: , ,

You can comment below, or link to this permanent URL from your own site.

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

  1. valdez Says:

    Hello,

    this is ecellent material. But could you provide some test input.pcm file? I don’t know how that file should look like. How to generate it?

    • Shawn Says:

      I haven’t provided the test files because I don’t have the permissions to post them at this site. In the future (depending on me finding time) I may store them at another site and provide a link. You can generate test files in Matlab, or with an audio program such as Adobe Audition, or there is a free program called Audacity. Audition will read and write headerless files like I use in the test program (I think it calls them “raw pcm”). Audacity probably can as well. You need to make the audio 16 bits per sample with an 8000 Hz sampling rate (and the data must be binary, not ASCII text).

      I created the files myself with CoolEdit 96, a predecessor to Adobe Audition (it was a shareware program). And yes, I have been using that program since about 1996. You can search for CoolEdit on the web and you may find a version of it.

      • valdez Says:

        Thanks for your reply. Finally I’d like to use Matlab, but first I downloaded that CoolEdit software. To see how such a file looks like. There is one another option for ‘raw pcm’ files, bytes order. I suppose it should be MSB, LSB order. Am I right?

        In such a case, if I’d like to test your filter, the first approach is to process “delta” signal to get coefficients values in the output file. So that test file should look like: first value 32 768, next many zeros. Is it correct?

  2. valdez Says:

    Because I need fixed point, I started to testing your code from second part. 32767 written as 0xFF7F worked as impulse signal I after filtration I get coefficients.
    Now I’m trying to generate so coded signal in Matlab.
    Thank you for this site 🙂

    • Shawn Says:

      The MSB, LSB order depends on the type of CPU you are using (big endian or little endian). On an Intel machine, which is little endian, the value 32768 would look like 0xFF7F instead of 0x7FFF if viewed byte by byte, because the 8 least significant bytes are stored at the lowest address. Your idea of using an impulse is a good way to test; it’s one I often use myself.

  3. umair Says:

    Hi . I have another question related to fir filter. What if I want to implement a Cascade of Fir filters. For example cascade of 3 fir filters how would the code look like any idea or at least give some pseudo code. One more thing , why isn’t the memory of the FIR filter not taken into consideration in this code.

    • Shawn Says:

      Hi. Your two questions are related. To cascade more than one FIR filter together, you would need to modify the firFloat function to take a pointer to the filter history as an argument. In the example code, I have one fixed buffer to hold history (called insamp). And then the firFloatInit would also need to have a pointer to the filter history passed in. When running multiple filters, you would need to call firFloatInit for each filter. Then you would need to run firFloat multiple times, first with floatInput as input and floatOutput as output. Next with floatOutput as input and floatOutput2 as output, etc.

      To make things cleaner, you might want to consider defining a filter structure that includes the coefficients, the filter length and the filter memory, and then passing a pointer to that structure for the init and filter functions.

      • umair Says:

        Thank you for your answer. I can understand it but can you give some sample code i mean how to do it by structures? I am trying for 2 weeks but couldn’t figure out what to do. Your help would be highly appreciated 🙂

  4. Shawn Says:

    There is some sample code for how to use structures with an FIR filter in one of my other tutorials. Go to https://sestevenson.wordpress.com/tone-detection-using-the-teager-kaiser-energy-operator-part-2/ and look for the code listing for tone_detect.c. From line 168 to the bottom are init and FIR functions. Look at the tone_detect.h listing to see how the FIRClass structure is defined.

  5. Bibek Tiwari Says:

    This work is really commendable. I really appreciate it but i wonder about the coefficients of FIR filter. Would you mind sharing the idea about generating the coefficients of FIR filter, any useful sites or sample program. ?

    Sincerly,
    Bibek Tiwari

    • Shawn Says:

      I designed the bandpass filter in the sample code using a program called Scilab. It’s a math scripting program similar to Matlab, but definitely different than Matlab. I would suggest using Matlab and it’s filter design packages if you have them (or access to them). Other alternatives are Octave (quite similar to Matlab) and Scilab. Scilab has some detailed documentation on how to design filters using their software. Here is my Scilab script:

      // bandpass filter centred at 1000 Hz

      nf=63;
      bedge=[0,800/8000;900/8000,1100/8000;1200/8000,4000/8000];
      des=[0 1 0];
      wate=[1 1 1];
      [hn]=eqfir(nf,bedge,des,wate);
      [hm,fr]=frmag(hn,256);
      plot(fr,hm)

      print(‘coeff.txt’,hn)


Leave a Reply

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

WordPress.com Logo

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

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s


%d bloggers like this: