# Shawn's DSP Tutorials

tutorials on practical aspects of digital signal processing

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

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 155 156 157 158 159 160 161 162 163 164 | `#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.

### Like this:

Like Loading...

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

1. Abdul Sattar Says:
   April 7, 2012 at 6:12 am

   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.

   Reply
   - Shawn Says:
     April 11, 2012 at 8:48 am

     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.

     Reply
     - Abdul Sattar Says:
       April 12, 2012 at 3:59 am

       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

       #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:
       April 12, 2012 at 9:23 pm

       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:
   April 12, 2012 at 10:47 pm

   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.

   Reply
   - Shawn Says:
     April 18, 2012 at 2:40 pm

     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!

     Reply
     - Abdul Sattar Says:
       April 19, 2012 at 8:25 am

       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:
   January 29, 2014 at 6:57 am

   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?

   Reply
   - Shawn Says:
     February 1, 2014 at 9:36 pm

     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.

     Reply
4. Slayer Says:
   March 14, 2014 at 9:06 am

   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

   Reply
   - Shawn Says:
     March 14, 2014 at 9:57 am

     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.

     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: