Fixedpoint FIR filter in C?

Question

There is an algorithm for a FIR filter but it's floatingpoint: FIR filter implementation in C programming

If I want a fixedpoint algorithm with this spec, how would I do it?

the FIR-filter receives and sends 8-bit fixed-point numbers in the Q7-format via the standard input and output. Remember to output the measured time (number of ticks) also in hex format. Following the guidelines presented in the previous section, your program should call getchar() to read a Q7-value. should call putchar() to write a Q7-value.

The coefficients are

c0=0.0299 c1=0.4701 c2=0.4701 c3=0.299

And for a fixedpoint algorithm I would need to implement my own multiplication for fixedpoint number, right?

Should I store a fixepdpoint number like a struct i.e.

struct point
{
    int integer;
    int fraction;
};

Should I use shifts to implement the numbering and specifically how?

The number are 32-bit so could I write the shifts like below?

#define SHIFT_AMOUNT 16 // 2^16 = 65536
#define SHIFT_MASK ((1 << SHIFT_AMOUNT) - 1)

So I think that I must implement one multiplication algorithm and then the FIR algorithm itself? Is that correct? Can you help me?

Update

I compiled and ran a program like in the anser but it's giving me unexpected output.

#include <stdio.h>
#include "system.h"

#define FBITS 16 /* number of fraction bits */
const int c0 = (( 299<<FBITS) + 5000) / 10000; /* (int)(0.0299*(1<<FBITS) + 0.5) */
const int c1 = ((4701<<FBITS) + 5000) / 10000; /* (int)(0.4701*(1<<FBITS) + 0.5) */
/* Ditto for C3 and C2 */
const int c2 = (( 4701<<FBITS) + 5000) / 10000; /* (int)(0.4701 *(1<<FBITS) + 0.5) */
const int c3 = ((299<<FBITS) + 5000) / 10000; /* (int)(0.299*(1<<FBITS) + 0.5) */

#define HALF (1 << (FBITS) >> 1) /* Half adjust for rounding = (int)(0.5 * (1<<FBITS)) */
signed char input[4]; /* The 4 most recent input values */
int output = 0;

void firFixed()
{
 signed char sum = c0*input[0] + c1*input[1] + c2*input[2] + c3*input[3];
 output = (signed char)((sum + HALF) >> FBITS);
 printf("output: %d\n", output);
}

int main( void )
{   
    int i=0;
    signed char inVal;
    while (scanf("%c", &inVal) > 0)
    {    
     if (i>3)
     {
       i=0;
     }       
     input[i]=inVal;           
     firFixed();
     i++;      
    }
    return 0;
}

Why is output not computed correctly andd why is output written several times after one input?

Update

I tried writing the fixedpoint FIR filter, the algorithm might not be 100 % correct:

    #include <stdio.h>
    #include "system.h"

    #define FBITS 16 /* number of fraction bits */
    const int c0 = (( 299<<FBITS) + 5000) / 10000; /* (int)(0.0299*(1<<FBITS) + 0.5) */
    const int c1 = ((4701<<FBITS) + 5000) / 10000; /* (int)(0.4701*(1<<FBITS) + 0.5) */
    /* Ditto for C3 and C2 */
    const int c2 = (( 4701<<FBITS) + 5000) / 10000; /* (int)(0.4701 *(1<<FBITS) + 0.5) */
    const int c3 = ((299<<FBITS) + 5000) / 10000; /* (int)(0.299*(1<<FBITS) + 0.5) */

    #define HALF (1 << (FBITS) >> 1) /* Half adjust for rounding = (int)(0.5 * (1<<FBITS)) */
    signed char input[4]; /* The 4 most recent input values */

    char get_q7( void );
    void put_q7( char );

    void firFixed()
    {
     int sum = c0*input[0] + c1*input[1] + c2*input[2] + c3*input[3];
     signed char output = (signed char)((sum + HALF) >> FBITS);
     put_q7(output);
    }

    int main( void )
    {   
        int i=0;
        while(1)
        {    
         if (i>3)
         {
           i=0;
         }       
         input[i]=get_q7();           
         firFixed();
         i++;      
        } 
        return 0;
    }

#include <sys/alt_stdio.h>

char get_q7( void );

char prompt[] = "Enter Q7 (in hex-code): ";
char error1[] = "Illegal hex-code - character ";
char error2[] = " is not allowed";
char error3[] = "Number too big";
char error4[] = "Line too long";
char error5[] = "Line too short";

char get_q7( void )
{
    int c; /* Current character */
    int i; /* Loop counter */
    int num;
    int ok = 0; /* Flag: 1 means input is accepted */

    while( ok == 0 )
    {
        num = 0;
        for( i = 0; prompt[i]; i += 1 )
            alt_putchar( prompt[i] );

        i = 0; /* Number of accepted characters */
        while( ok == 0 )
        {
            c = alt_getchar();
            if( c == (char)26/*EOF*/ ) return( -1 );
            if( (c >= '0') && (c <= '9') )
            {
                num = num << 4;
                num = num | (c & 0xf);
                i = i + 1;
            }
            else if( (c >= 'A') && (c <= 'F') )
            {
                num = num << 4;
                num = num | (c + 10 - 'A');
                i = i + 1;
            }
            else if( (c >= 'a') && (c <= 'f') )
            {
                num = num << 4;
                num = num | (c + 10 - 'a');
                i = i + 1;
            }
            else if( c == 10 ) /* LF finishes line */
            {
                if( i > 0 ) ok = 1;
                else
                {    /* Line too short */
                    for( i = 0; error5[i]; i += 1 )
                        alt_putchar( error5[i] );
                    alt_putchar( '\n' );
                    break; /* Ask for a new number */
                }
            }
            else if( (c & 0x20) == 'X' || (c < 0x20) )
            {
                /* Ignored - do nothing special */
            }
            else
            {   /* Illegal hex-code */
                for( i = 0; error1[i]; i += 1 )
                    alt_putchar( error1[i] );
                alt_putchar( c );
                for( i = 0; error2[i]; i += 1 )
                    alt_putchar( error2[i] );
                alt_putchar( '\n' );
                break; /* Ask for a new number */
            }
            if( ok )
            {
                if( i > 10 )
                {
                    alt_putchar( '\n' );
                    for( i = 0; error4[i]; i += 1 )
                        alt_putchar( error4[i] );
                    alt_putchar( '\n' );
                    ok = 0;
                    break; /* Ask for a new number */
                }
                if( num >= 0 && num <= 255 )
                    return( num );
                for( i = 0; error3[i]; i += 1 )
                    alt_putchar( error3[i] );
                alt_putchar( '\n' );
                ok = 0;
                break; /* Ask for a new number */
            }
        }
    }
    return( 0 ); /* Dead code, or the compiler complains */
}


#include <sys/alt_stdio.h>

void put_q7( char );    /* prototype */

char prom[] = "Calculated FIR-value in Q7 (in hex-code): 0x";

char hexasc (char in)   /* help function */
{
    in = in & 0xf;
    if (in <=9 ) return (in + 0x30);
    if (in > 9 ) return (in - 0x0A + 0x41);
    return (-1);
}

void put_q7( char inval)
{
    int i; /* Loop counter */   
        for( i = 0; prom[i]; i += 1 )
            alt_putchar( prom[i] );
    alt_putchar (hexasc ((inval & 0xF0) >> 4));
    alt_putchar (hexasc (inval & 0x0F));
    alt_putchar ('\n');     
}

Answer 1

Each point in the FIR filter result is just a weighted sum of values from the unfiltered data. You shouldn't need anything other than plain multiplication and addition if you have 8-bit input data and 32-bit arithmetic.

A quick visit to Wikipedia tells me that Q7 is essentially an 8-bit 2's complement integer, so if the target platform uses 2's complement, then simply describing the byte received as (signed char) will give it the correct numerical value when promoted to an int. If you premultiply the coefficients by a power of two, then the weighted sum will be multiplied by that same power of 2. Rounded division is then simply adding a half-adjust value followed by a signed right shift. With 16-bit fractions, the premultiplied constants are:

#define FBITS 16 /* number of fraction bits */
const int C0 = (( 299<<FBITS) + 5000) / 10000; /* (int)(0.0299*(1<<FBITS) + 0.5) */
const int C1 = ((4701<<FBITS) + 5000) / 10000; /* (int)(0.4701*(1<<FBITS) + 0.5) */
/* Ditto for C3 and C2 */
#define HALF (1 << (FBITS) >> 1) /* Half adjust for rounding = (int)(0.5 * (1<<FBITS)) */

The reason for that oddness it to get the significant bits you want without depending on any floating point rounding. Now, if:

signed char input[4];

...contain the 4 most recent input values, your output value is:

sum = c0*input[0] + c1*input[1] + c2*input[2] + c3*input[3];
output = (signed char)((sum + HALF) >> FBITS);

Since all of your coefficients are positive and sum to 1.0, there's no possibility of overflow.

There are a number of optimizations you can try after you get a simple version working. One possible minor glitch with other coefficients is for that rounding of the C0-C3 constants to produce values that don't exactly add up to 1<<FBITS. I tested and it doesn't happen with these values (you'd need the c0*(1<<LBITS) to have a fraction part of exactly 0.5; meaning that all the other scaled coefficients would also have 0.5 as their fraction parts. They'd all round up and the sum would be too large by 2. That could add a very small unintended gain to your filter.

That can't occur with the coefficients you gave.

Edit: I forgot. Both the integer part and fraction part are in the same 32-bit int during the sum calculation. With 8 bits of input(7+sign) and 16 bits of fraction, you can have up to 2^(32 - 16 - 8) = 2^8 = 256 points in the filter (at this point, you will obviously have an array of coefficients, and a multiply-add loop to compute the sum. Should the (input size) + (fraction bits) + log2(filter size) exceed 32, then you can try to expand the sum field to a C99 long long or int64_t value, if that's available, or write extended-precision add and shift logic if not. Extended precision in hardware is far better to use, when available.