Implementation of FIR Filtering in C (Part 3)

Part 2 showed an example of a FIR filter in C using fixed point. This tutorial on FIR filtering shows how to apply several different FIR filters to the same input data. The examples for this part are also in fixed point. The example is a single C file with the FIR filter code at the top, and a small test program at the bottom. In an actual implementation, you would likely want to split the code into several files.

Click the following link for a PDF version of the code example:

Code Example PDF

The code example is shown below:

#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
int16_t insamp[ BUFFER_LEN ];

// FIR init
void firFixedInit( void )
{
    memset( insamp, 0, sizeof( insamp ) );
}

// store new input samples
int16_t *firStoreNewSamples( int16_t *inp, int length )
{
    // put the new samples at the high end of the buffer
    memcpy( &insamp[MAX_FLT_LEN - 1], inp,
            length * sizeof(int16_t) );
    // return the location at which to apply the filtering
    return &insamp[MAX_FLT_LEN - 1];
}

// move processed samples
void firMoveProcSamples( int length )
{
    // shift input samples back in time for next time
    memmove( &insamp[0], &insamp[length],
            (MAX_FLT_LEN - 1) * sizeof(int16_t) );
}

// 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;

    // apply the filter to each input sample
    for ( n = 0; n < length; n++ ) {
        // calculate output n
        coeffp = coeffs;
        inputp = &input[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);
    }
}

//////////////////////////////////////////////////////////////
//  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
};

// Moving average (lowpass) filter of length 8
// There is a null in the spectrum at 1000 Hz

#define FILTER_LEN_MA   8
int16_t coeffsMa[ FILTER_LEN_MA ] =
{
    32768/8, 32768/8, 32768/8, 32768/8,
    32768/8, 32768/8, 32768/8, 32768/8
};

// number of samples to read per loop
#define SAMPLES   80

int main( void )
{
    int size;
    int16_t input[SAMPLES];
    int16_t output[SAMPLES];
    int16_t *inp;
    FILE   *in_fid;
    FILE   *out_fid;
    FILE   *out_fid2;

    // 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 files
    out_fid = fopen( "outputFixed.pcm", "wb" );
    if ( out_fid == 0 ) {
        printf("couldn't open outputFixed.pcm");
        return;
    }
    out_fid2 = fopen( "outputFixedMa.pcm", "wb" );
    if ( out_fid == 0 ) {
        printf("couldn't open outputFixedMa.pcm");
        return;
    }

    // initialize the filter
    firFixedInit();

    // process all of the samples
    do {
        // read samples from file
        size = fread( input, sizeof(int16_t), SAMPLES, in_fid );
        // store new samples in working array
        inp = firStoreNewSamples( input, size );

        // apply each filter
        firFixed( coeffs, inp, output, size, FILTER_LEN );
        fwrite( output, sizeof(int16_t), size, out_fid );

        firFixed( coeffsMa, inp, output, size, FILTER_LEN_MA );
        fwrite( output, sizeof(int16_t), size, out_fid2 );

        // move processed samples
        firMoveProcSamples( size );
    } while ( size != 0 );

    fclose( in_fid );
    fclose( out_fid );
    fclose( out_fid2 );

    return 0;
}

There are a few differences from the code example of Part 2. First, I have created a function to store the input samples to the input sample array (firStoreNewSamples). This function is called once for every block of input samples that are processed. The calling function passes in a pointer to the new input samples, and the number of new samples to copy. The function returns the address at which to apply the FIR filter.

Second, I have added a function to move the samples after processing a block of samples (firMoveProcSamples). Again, this function is called once per block of samples, not once per FIR filter applied.

The FIR filtering function (firFixed) has the same argument list as in the Part 2 example, but the “input” argument is a bit different in this case. The input pointer passed in should be the address returned from the firStoreNewSamples function, rather than a pointer to the input sample buffer.

The test program shows an example where two different FIR filters are applied to the same output data. First one input file is opened (for input samples) and two output files are opened (one for each filter). In the sample processing loop, a block of up to 80 samples is read and stored into the working array for the filters. Next the 63 tap bandpass filter is applied by calling firFixed, and the block of output samples is written to file. Afterwards, the 8 tap moving average filter is applied, and the output samples are written to a different file. Finally, the sample buffer is shifted to prepare for the next block of input samples.

The code I have shown works for however many filters that you want to implement. Remember to keep track of the maximum filter tap length, and input sample block size, and change the #define statements appropriately. That concludes my tutorial on basic FIR filters.

Advertisements
Explore posts in the same categories: FIR filters, fixed point

12 Comments on “Implementation of FIR Filtering in C (Part 3)”

  1. Dat Says:

    Hi, Thanks for the great post on fixed point fir.
    I am coding for a pic microchip and can really use
    the fixed point implementation. Can you please tell me how I can calculate the coefficients for use with your filter? I will need to use it for a bandpass and notch filter. Thanks.

    • Shawn Says:

      I would recommend you to use Matlab to generate your filter coefficients if you have access to it. I calculated the coefficients for the example using a program called Scilab. (Scilab is a free alternative to Matlab.) Another program you might try is called Octave. If you go with Scilab, do a search for the document called “Signal Processing with Scilab”. It has a bunch of filter design examples in it.

      All of these programs will include algorithms/routines that will allow you to generate filter coefficients. Make sure the maximum gain of the filter is not too much larger than unity, and that none of the coefficients is greater than or equal to 1 (or less than -1). Once you have your coefficients in floating point, multiply each one of them by 32768 and round to an integer.

      Below is some Scilab code to generate a bandpass filter at 1000 Hz.

      Good luck with your filters!

      // 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);
      hn=hn*32768;
      hn=round(hn);
      [hm,fr]=frmag(hn,256);
      plot(fr,hm)

      print(‘coeff.txt’,hn)


  2. […] tutorials to get you started: Introduction to FIR Digital Filters The simple FIR filter equation Implementation of FIR Filtering in C FIR Filter Implementation An FIR Filter Example Hope these links help. Reply With […]

  3. Lee Says:

    Hey Shawn. GREAT post! I’m trying to implement FIR on a wav format file using your coding, but the output wav seems to be unreadable. My windows media player gave me this error message:

    Windows Media Player cannot play the file. The Player might not support the file type or might not support the codec that was used to compress the file.

    Do you have a clue on that? Thank you.

    • Shawn Says:

      Hi Lee. The code works on headerless files. The headerless files are basically the same as a mono, uncompressed .wav file, but without the file header. My code doesn’t read wav file headers or write them for the output file. So you have a couple of options:

      The first option is to write code to read the wav file header of the input file (in which you can verify that its mono, uncompressed, 16 bits). And write a wav file header before writing out the output samples. Here’s a good link about the wav file headers: http://www.sonicspot.com/guide/wavefiles.html

      The second option is to use headerless (sometimes called raw) files. There are a number of programs that will play these files, convert to wav and back, etc. The program I use at work is Adobe Audition. There is also a free program called Audacity (I haven’t tried it myself). There is also an old shareware program called CoolEdit (which eventually became Adobe Audition) and you may be able to find a version of that on the web for free. You won’t be able to register it, but the unregistered versions can be quite usable (for example letting you pick 2 out of 6 different sets of options each time you start).

      Hope that helps.

  4. Lee Says:

    Thank you Shawn! Thank you so much for replying 🙂

    I’ve tried out the first option you suggested, it goes all well for me to read the wav header of input but not at writing a wav header for the output. (p/s: I’m just not too good at programming or writing codes =P)

    Do you have any tutorial on handling wav header files?
    Thanks again!

    • Shawn Says:

      For this particular case, you can write the header data to the output file as you are reading it from the input file. In other words you can copy the input file header to the output file. It should work okay since the output file will have the same number of samples as the input file for the code I have written above.

      Another option I didn’t mention before is to use a library that someone else has written. There is one called “libsndfile” that works with all kinds of audio file formats (and is written in C).

  5. Brandon Says:

    This is a really great blog that’s helping me come up to speed on some DSP applications I’m assigned to at work. Thank you, and I hope you continue it.

  6. Steve Says:

    Nice tutorial. How did you generate the input test sample files?

    • Shawn Says:

      Read some of the previous comments (March 13, 2012), and you’ll see a description of how I generated the test sample files. I used a program called CoolEdit 95, which works fine with Windows 7, other than the fact that the help doesn’t work at all.


  7. i can’t run this code using turboc

    • Shawn Says:

      I compiled and linked the code with GCC when I wrote it. But there isn’t anything GCC specific in the code as far as I know. Does the code compile and link but not run? The header stdint.h was introduced with the C99 standard. Does your version of Turbo C support C99 or include stdint.h?


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: