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

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

  1. Abdul Sattar Says:

    Thank you so much for posting very nice tutorial. I would like to implement your example into AVR Atmega16 controller. I have the ADC function ReadADC(0x00) which reads the 10bits value from the ADC Registers. where should i input the ADC value to the above example. I am new to controller program, I shall be thankful for your guides and suggestions. Please guide me where i should make necessary changes in above example.

    • Shawn Says:

      Abdul, thanks for the compliment. You should read your ADC samples into an array and then pass the address of that array to the firFixed function as the second argument. In the main program, I read samples from a file into the array called “input”. In your code you should read ADC samples into a similar array. In the example I have processed 80 samples at a time, but you should change this value to what is appropriate for you application. Processing one sample at a time would minimize delay through the filter, but would be the most costly in terms of cycles.

      • Abdul Sattar Says:

        Dear Shawn, thank you for your reply. i have read ADC value as input[SAMPLES]=ADCRead(0x00), and passed the address of the input[SAMPLES] to the function firStoreNewSamples as inp = firStoreNewSamples(&input[SAMPLES],size). i am confused about the size integer. as size is the length, would you please suggest me what should be the value of size. since FirFixed function has five arguments as firFixed( coeffs, inp, output, size, FILTER_LEN ); and inp is read from firStoreNewSamples . My program is listed below kindly chek it. I would be thank full to you, thank you.

        #include
        #include
        #include

        #define ADC_VREF_TYPE 0x20
        #define PWM1DCReg OCR1A

        //////////////////////////////////////////////////////////////
        // Filter Code Definitions
        //////////////////////////////////////////////////////////////

        // maximum number of inputs that can be handled
        // in one function call
        #define MAX_INPUT_LEN 25
        // maximum length of filter than can be handled
        #define MAX_FLT_LEN 16
        // 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 ];

        #define SAMPLES 25

        // 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 0x3fffffff ) {
        acc = 0x3fffffff;
        } else if ( acc > 15);
        }
        }

        //Low pass
        //cutoff frequency 1000Hz, Sampling frequency 8915Hz.
        #define FILTER_LEN 16
        int16_t coeffs[ FILTER_LEN ] =
        {
        0xFF45, 0xFE83, 0xFDB2, 0xFFC9, 0x088E, 0x18A4,
        0x2B29, 0x37AF, 0x37AF, 0x2B29, 0x18A4, 0x088E,
        0xFFC9, 0xFDB2, 0xFE83, 0xFF45
        };

        void InitPort(void)
        {

        PORTA = 0x00;
        DDRA = 0x00;
        PORTB = 0x00;
        DDRB = 0xff;
        PORTC = 0x00;
        DDRC = 0xff;
        PORTD = 0x00;
        DDRD = 0x20; //PIND.5 Filter Output

        }

        void InitADC( void)
        {
        ADMUX = ADC_VREF_TYPE & 0xff; // Select channel 0
        ADCSRA = 0xC4; // Enable ADC & start 1:st dummy conversion
        }

        int16_t ReadADC( unsigned char channel)
        {
        int16_t ADC_Value;

        ADMUX = channel | (ADC_VREF_TYPE & 0xff); // Select channel

        _delay_us(10); // Delay needed for the stabilization of the ADC input voltage

        ADCSRA |= 0x40; // Start conversion

        while (!(ADCSRA & 0x10)); // Check if conversion is ready

        ADCSRA |= 0x10; // Clear Conversion ready flag by setting the bit

        ADC_Value = ADCL; // Read 8 low bits first (important)
        ADC_Value |= (int)ADCH << 8; // Read 2 high bits and multiply with 256

        return ADC_Value;

        }

        void timer1 (void)
        {
        // Timer/Counter 1 initialization
        // Clock source: System Clock
        // Clock value: 8000.000 kHz
        // Mode: 8 bit Pulse Width Modulation
        // OC1A output: Non-Inv.
        // OC1B output: Discon.
        // PWM output frequency is doubled
        // Noise Canceler: Off91,09
        // Input Capture on Falling Edge

        TCCR1A = 0x91;
        TCCR1B = 0x09;
        TCNT1H = 0x00;
        TCNT1L = 0x00;
        ICR1H = 0x00;
        ICR1L = 0x00;
        OCR1AH = 0x00;
        OCR1AL = 0x00;
        OCR1BH = 0x00;
        OCR1BL = 0x00;

        }

        int main( void )
        {

        // Initialize all the devices
        cli(); // Clear all Interrupts
        InitPort();
        InitADC();
        timer1();

        int size;
        int16_t input[SAMPLES];
        int16_t output[SAMPLES];
        int16_t *inp;

        // initialize the filter
        firFixedInit();

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

        // apply each filter
        firFixed( coeffs, inp, output, size, FILTER_LEN );

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

        // External Interrupt(s) initialization
        // INT0: Off
        // INT1: Off
        MCUCR = 0x00;
        MCUCSR = 0x00;

        // Timer(s)/Counter(s) Interrupt(s) initialization
        TIMSK = 0x83;

        // Analog Comparator initialization
        // Analog Comparator: Off
        // Analog Comparator Input Capture by Timer/Counter 1: Off
        ACSR = 0x80;
        SFIOR = 0x00;

        sei(); //Enable Global interrupt

        PWM1DCReg = output; // output for filter

        PORTB = ~PORTB; // Toggle output port to allow measurement of ISR / Sampling Rate
        PORTC = ADCH; // Check the ADCH Status result

        while(1);

        }

      • Shawn Says:

        There are a number of things wrong with your code. It appears you are reading only one sample at a time, so it would make sense to change the #define SAMPLES to 1 (or just use a variable instead of an array). The “size” was originally the number of samples read from file in my example. You should just pass 1 instead if you are doing 1 sample at a time. The do loop should be eliminated completely since you are no longer reading from a file. You should be either polling the timer interrupt or writing an ISR so that you can read from the ADC periodically (8000 or 8915 times per second, judging by comments in your code). And you probably need to write to your output register each time you have a new output.

        I would also recommend starting with a simpler program where you simply read the ADC and then output what you read. Once you have that working, put in the filter code.

        I hope that helps.

  2. Abdul Sattar Says:

    Dear Shawn, Thank you for your guide and suggestion. I am new to controller programming i must start first very simpler example. I am very much interested to learn programming the controller if possible please suggest me how should i start. I have Beginning knowledge of AVR architecture ATmega series 8-bit microcontrollers. and i am using AVR Studio for programming the controller. Please guide me how should i start working on controller.

    • Shawn Says:

      Abdul. My advice would be to start with some example programs that come with your development kit (or tools, or elsewhere) and make sure that you can get them working properly. Then try experimenting by modifying these programs.

      Once you get bored with that, try making your own program. Start with something simple like setting up the ADC and reading from it. Then add more functionality one step at a time. Add some code to set up the timer and verify that it works. Add code to poll the timer interrupt. Then try coding an interrupt service routine (ISR). Then add code to the ISR to read the ADC. The trick is to write a small amount of code at a time and verify that it works each step of the way. That is usually easier than writing a whole bunch of code at a time, and then trying to debug it all at once (especially if you are doing something that is new to you).

      Good luck!

      • Abdul Sattar Says:

        Dear Shawn, Thank you so much for your valued suggestion. And i look forward for your help if i found any problem in programming in future. Thank you again.

  3. Joseph Says:

    Shawn, I am trying to make a graphic equalizer using FIR filters and this has been quite helpful so far. But I was wondering if you have any information on how to apply gain to the output of each band?

    • Shawn Says:

      Hi Joseph. For a graphic equalizer, you want to apply a gain to each band by multiplying the end result of the filtering by a number between 0 and 1. So each filter should have a maximum gain of 1 (at the centre frequency) and then you should have a set of gains for each band. It is best to express the gains in dB (deciBels) if making a graphical display or control (like a slider). The gain in dB is 20*log10(g) where g is the multiplication value (between just above 0 and 1). Then -6 dB corresponds to a gain of about 0.5. That is 20*log10(0.5) is approximately -6. The reverse formula is g = pow10(G/20) where G is the gain in dB. So if the gain is -24 dB, the multiplier value is g = pow10(-24/20) = 0.063. The reason for using decibels is because the perception of loudness follows a logarithmic scale.

      I hope that helps you.

  4. Slayer Says:

    Hi Shawn, Great tutorial! Your blog has much more practical stuff than most signal processing books. I have a couple of follow up questions on your tutorial and it would be great if you could point me in the right direction.

    Am trying to implement(in hardware/verilog) a 32 bit mac with 2 16 bit inputs. If i want to store back the 32 bit result in a 16 bit register, I understand that for fixed point signed-signed multiplication we need to just (without rounding) extract bit 30 through 15 (and not 31 becuse of the double sign) .

    1) Am a bit confused if I need different hardware for signed fractional multiply and pure integer multiply. Say I use two separate 16 bit registers to store the 32bit result of the multiplication. Do you think just using a conditional shifter between the 32bit multiplier result (immediately after the multplier hardware) and just before storing it into 2 separate registers would be sufficient? The conditional shifter shifting left by 1 for fractional multiply and shifting by 0 for integer multiply?

    Thanks a ton 🙂

    • Shawn Says:

      Hi Slayer,

      Thanks for the comment. Your second question seems to have been lost. To answer your first question, you should be able to use the same multiplier for both integer and fixed point multiplication, just as you have described. For the fixed point multiplier, there is one special case that you need to be aware of, which is multiplying -1 by -1 (where -1 is 0x8000 in fractional Q.15). When you multiply 0x8000 by 0x8000, the result will be 0x4000,0000. If you shift that left you end up with 0x8000,0000, so the end result is that -1 times -1 is equal to -1. Some ways to handle this are:

      1) Never multiply -1 by -1
      2) Saturate the result to 0x7FFF,FFFF
      3) Raise an overflow flag and handle the error somehow.

      I hope that helps.


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: