Writing a Basic FIR Filter in C (1)

Introduction

A FIR (Finite Impulse Response) filter is mathematically identical to a convolution between the input array and a set of taps which represent the impulse response of the filter. This can be implemented using a series of dot products between the tap vector and a vector consisting of the current input sample and those which came immediately before it.

\[y_n = \sum_{N=0}^{N-1} x_{n-N} * b_N\]

where \(y\) is the output array, \(x\) is the input array, \(b\) is the array of taps and \(N\) is the number of taps.

When implementing this in code we need to decide how to treat the first \(N\) samples as we have no indication of what data should come before them. We can take one of two methods, either we can assume all previous samples are 0 or we can start our calculation at the \(N^{th}\) sample. For this example we will choose the zero assumption method as this gives us the same number of output samples as input.

Basic FIR Filter

We will first create an unoptimised FIR filter which cab be called multiple times as data comes in, with memory of the previous samples. This will require 4 functions, one to initialize the required memory, one to apply the FIR filter, one to reset the memory in the case of a discontinuity in the received data, and one to free all the filter memory. To pass pass around this memory we will first declare the following structure:

typedef struct FIR_Stream_Mem {
    size_t tap_len;
    double *taps;
    double *buffer;
} FIR_Stream_Mem;

The tap_len field contains the number of taps stored in the taps array. buffer will contain the previous samples and will also have size tap_len. We will allocate and populate this structure using the following function:

FIR_Stream_Mem *init_fir_stream(double* taps, size_t tap_len){
    /* Create the output object in memory */
    FIR_Stream_Mem *output = (FIR_Stream_Mem*)(malloc(sizeof(FIR_Stream_Mem)));

    /* Copy the taps across */
    output->tap_len = tap_len;
    output->taps = (double*)(malloc(sizeof(double) * tap_len));
    memcpy(output->taps, taps, tap_len * sizeof(double));

    /* Set the buffer to all zeros */
    output->buffer = (double*)(malloc(sizeof(double) * tap_len));
    memset(output->buffer, 0x00, tap_len * sizeof(double));

    return output;
}

This first allocates the memory to store the FIR_Stream_Mem structure itself, before allocating memory for the buffer and taps. The taps are copied over from the input to allow for the structure to be passed around more easily, and the buffer is set to all zeros as discussed above.

We are now ready to start using the FIR filter. To do this, we need to implement the algorithm so it can be run on our data. We will use the following function:

void apply_fir(FIR_Stream_Mem *mem, const double *input, double *output, size_t size){
    /* Loop through all input samples */
    for(size_t i = 0; i < size; i++){
        /* Calculate first output sample */
        double sample_output = input[i] * mem->taps[0];

        /* Loop through all the taps */
        for(size_t j = 1; j < mem->tap_len; j++){
            sample_output += mem->buffer[j] * mem->taps[j];
        }

        /* Update buffer */
        for(size_t j = mem->tap_len - 2; j > 0; j--){
            mem->buffer[j + 1] = mem->buffer[j];
        }
        mem->buffer[1] = input[i];
        output[i] = sample_output;
    }
}

This has the following arguments:

  • mem : A pointer to the FIR_Stream_Mem struct for the filter to be applied.

  • input : A pointer to an array of input samples.

  • output : A pointer to an array of output samples. This must be pre-allocated with a size greater than or equal to the size of the input array.

  • size : The number of samples to process. Note this must be less than or equal to the size of the input array

It loops through all the input samples and calculates the dot product as discussed above. It then shifts all the samples in the buffer to the right and adds the current sample to the front before calculating the next sample.

This function can now be applied on chunks of our input data as they come in, with the memory from the last chunk being kept in the FIR_Stream_Mem struct.

However, if a discontinuity is detected in the input stream of samples the history stored in buffer will be incorrect and we will have no way of determining what it should be. In this case it is best to reset the filter by setting the filter’s buffer to all zeros, which can be done with the following function:

void reset_fir(FIR_Stream_Mem *mem){
    memset(mem->buffer, 0x00, mem->tap_len * sizeof(double));
}

Finally, when we are done with the filter we should free all associated memory. This can be done with the following function:

void deinit_fir(FIR_Stream_Mem *mem){
    free(mem->buffer);
    free(mem->taps);
    free(mem);
}

Conclusion

We have created a FIR filter in C which can be applied to a stream of samples coming into our program. However, there are some optimizations which can be made to our code - these will be covered in part 2.