FIR Filter

Table of Contents

    FIR Filter Design with Python

    In this paragraph we see how to design a FIR filter using SciPy, NumPy and MatplotLib.

    You can install SciPy, NumPy and MatplotLib as normal Python libraries, otherwise you can simply install the Enthought Python distribution.

    The Enthought Python distribution is packaging for you a Python interpreter and all the above libraries in a single convenient installer: http://enthought.com/products/epd_free.php

    Python FIR Filter Design

    from numpy import sin, arange, pi
    from scipy.signal import lfilter, firwin
    from pylab import figure, plot, grid, show
    
    #------------------------------------------------
    # Create a signal for demonstration.
    #------------------------------------------------
    # 320 samples of (1000Hz + 15000 Hz) at 48 kHz
    sample_rate = 48000.
    nsamples = 320
    
    F_1KHz = 1000.
    A_1KHz = 1.0
    
    F_15KHz = 15000.
    A_15KHz = 0.5
    
    t = arange(nsamples) / sample_rate
    signal = A_1KHz * sin(2*pi*F_1KHz*t) + A_15KHz*sin(2*pi*F_15KHz*t)
    
    #------------------------------------------------
    # Create a FIR filter and apply it to signal.
    #------------------------------------------------
    # The Nyquist rate of the signal.
    nyq_rate = sample_rate / 2.
    
    # The cutoff frequency of the filter: 6KHz
    cutoff_hz = 6000.0
    
    # Length of the filter (number of coefficients, i.e. the filter order + 1)
    numtaps = 29
    
    # Use firwin to create a lowpass FIR filter
    fir_coeff = firwin(numtaps, cutoff_hz/nyq_rate)
    
    # Use lfilter to filter the signal with the FIR filter
    filtered_signal = lfilter(fir_coeff, 1.0, signal)
    
    #------------------------------------------------
    # Plot the original and filtered signals.
    #------------------------------------------------
    
    # The first N-1 samples are "corrupted" by the initial conditions
    warmup = numtaps - 1
    
    # The phase delay of the filtered signal
    delay = (warmup / 2) / sample_rate
    
    figure(1)
    # Plot the original signal
    plot(t, signal)
    
    # Plot the filtered signal, shifted to compensate for the phase delay
    plot(t-delay, filtered_signal, 'r-')
    
    # Plot just the "good" part of the filtered signal.  The first N-1
    # samples are "corrupted" by the initial conditions.
    plot(t[warmup:]-delay, filtered_signal[warmup:], 'g', linewidth=4)
    
    grid(True)
    
    show()
    
    #------------------------------------------------
    # Print values
    #------------------------------------------------
    def print_values(label, values):
        var = "float32_t %s[%d]" % (label, len(values))
        print "%-30s = {%s}" % (var, ', '.join(["%+.10f" % x for x in values]))
    
    print_values('signal', signal)
    print_values('fir_coeff', fir_coeff)
    print_values('filtered_signal', filtered_signal)
    

    coefficients

    float32_t fir_coeff[29]        = {-0.0018225230, -0.0015879294, +0.0000000000, +0.0036977508, +0.0080754303, +0.0085302217, -0.0000000000, -0.0173976984, -0.0341458607, -0.0333591565, +0.0000000000, +0.0676308395, +0.1522061835, +0.2229246956, +0.2504960933, +0.2229246956, +0.1522061835, +0.0676308395, +0.0000000000, -0.0333591565, -0.0341458607, -0.0173976984, -0.0000000000, +0.0085302217, +0.0080754303, +0.0036977508, +0.0000000000, -0.0015879294, -0.0018225230}
    

    /media/uploads/emilmont/dsp_fir.png

    mbed-dsp library

    Import librarymbed-dsp

    CMSIS DSP library

    CMSIS DSP API

    Import program

    00001 #include "arm_math.h" 
    00002 #include "math_helper.h" 
    00003 #include <stdio.h>
    00004 
    00005 #define BLOCK_SIZE          32
    00006 #define NUM_BLOCKS          10
    00007 
    00008 #define TEST_LENGTH_SAMPLES (BLOCK_SIZE * NUM_BLOCKS)
    00009 
    00010 #define SNR_THRESHOLD_F32   140.0f
    00011 #define NUM_TAPS            29
    00012 
    00013 /* -------------------------------------------------------------------
    00014  * The input signal and reference output (computed with MATLAB)
    00015  * are defined externally in arm_fir_lpf_data.c.
    00016  * ------------------------------------------------------------------- */
    00017 extern float32_t testInput_f32_1kHz_15kHz[TEST_LENGTH_SAMPLES]; 
    00018 extern float32_t refOutput[TEST_LENGTH_SAMPLES]; 
    00019 
    00020 /* ------------------------------------------------------------------- 
    00021  * Declare State buffer of size (numTaps + blockSize - 1)
    00022  * ------------------------------------------------------------------- */
    00023 static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1]; 
    00024 
    00025 /* ---------------------------------------------------------------------- 
    00026  * FIR Coefficients buffer generated using fir1() MATLAB function. 
    00027  * fir1(28, 6/24)
    00028  * ------------------------------------------------------------------- */ 
    00029 const float32_t firCoeffs32[NUM_TAPS] = { 
    00030     -0.0018225230f, -0.0015879294f, +0.0000000000f, +0.0036977508f, +0.0080754303f,
    00031     +0.0085302217f, -0.0000000000f, -0.0173976984f, -0.0341458607f, -0.0333591565f,
    00032     +0.0000000000f, +0.0676308395f, +0.1522061835f, +0.2229246956f, +0.2504960933f,
    00033     +0.2229246956f, +0.1522061835f, +0.0676308395f, +0.0000000000f, -0.0333591565f,
    00034     -0.0341458607f, -0.0173976984f, -0.0000000000f, +0.0085302217f, +0.0080754303f,
    00035     +0.0036977508f, +0.0000000000f, -0.0015879294f, -0.0018225230f
    00036 };
    00037 
    00038 /* ---------------------------------------------------------------------- 
    00039  * FIR LPF Example 
    00040  * ------------------------------------------------------------------- */
    00041 int32_t main(void) { 
    00042     /* Call FIR init function to initialize the instance structure. */
    00043     arm_fir_instance_f32 S;
    00044     arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32[0], &firStateF32[0], BLOCK_SIZE); 
    00045     
    00046     /* ---------------------------------------------------------------------- 
    00047      * Call the FIR process function for every blockSize samples  
    00048      * ------------------------------------------------------------------- */ 
    00049     for (uint32_t i=0; i < NUM_BLOCKS; i++) {
    00050         float32_t* signal =  testInput_f32_1kHz_15kHz + (i * BLOCK_SIZE);
    00051         arm_fir_f32(&S, signal, signal, BLOCK_SIZE);
    00052     } 
    00053     
    00054     /* ---------------------------------------------------------------------- 
    00055      * Compare the generated output against the reference output computed
    00056      * in MATLAB.
    00057      * ------------------------------------------------------------------- */
    00058     float32_t snr = arm_snr_f32(refOutput, testInput_f32_1kHz_15kHz, TEST_LENGTH_SAMPLES);
    00059     printf("snr: %f\n\r", snr);
    00060     if (snr < SNR_THRESHOLD_F32) {
    00061         printf("Failed\n\r");
    00062     } else {
    00063         printf("Success\n\r");
    00064     }
    00065 } 
    

    mbed DSP API

    Import program

    00001 #include "mbed.h"
    00002 #include "dsp.h"
    00003 
    00004 #define BLOCK_SIZE              (32)
    00005 #define NUM_BLOCKS              (10)
    00006 #define TEST_LENGTH_SAMPLES     (BLOCK_SIZE * NUM_BLOCKS)
    00007 
    00008 #define SAMPLE_RATE             (48000)
    00009 
    00010 #define SNR_THRESHOLD_F32       (50.0f)
    00011 
    00012 float32_t expected_output[TEST_LENGTH_SAMPLES];
    00013 float32_t          output[TEST_LENGTH_SAMPLES];
    00014 
    00015 /* FIR Coefficients buffer generated using fir1() MATLAB function: fir1(28, 6/24) */
    00016 #define NUM_TAPS            29
    00017 const float32_t firCoeffs32[NUM_TAPS] = { 
    00018     -0.0018225230f, -0.0015879294f, +0.0000000000f, +0.0036977508f, +0.0080754303f,
    00019     +0.0085302217f, -0.0000000000f, -0.0173976984f, -0.0341458607f, -0.0333591565f,
    00020     +0.0000000000f, +0.0676308395f, +0.1522061835f, +0.2229246956f, +0.2504960933f,
    00021     +0.2229246956f, +0.1522061835f, +0.0676308395f, +0.0000000000f, -0.0333591565f,
    00022     -0.0341458607f, -0.0173976984f, -0.0000000000f, +0.0085302217f, +0.0080754303f,
    00023     +0.0036977508f, +0.0000000000f, -0.0015879294f, -0.0018225230f
    00024 };
    00025 #define WARMUP    (NUM_TAPS-1)
    00026 #define DELAY     (WARMUP/2)
    00027 
    00028 int main() {
    00029     Sine_f32 sine_1KHz(  1000, SAMPLE_RATE, 1.0);
    00030     Sine_f32 sine_15KHz(15000, SAMPLE_RATE, 0.5);
    00031     FIR_f32<NUM_TAPS> fir(firCoeffs32);
    00032     
    00033     float32_t buffer_a[BLOCK_SIZE];
    00034     float32_t buffer_b[BLOCK_SIZE];
    00035     for (float32_t *sgn=output; sgn<(output+TEST_LENGTH_SAMPLES); sgn += BLOCK_SIZE) {
    00036         sine_1KHz.generate(buffer_a);           // Generate a 1KHz sine wave
    00037         sine_15KHz.process(buffer_a, buffer_b); // Add a 15KHz sine wave
    00038         fir.process(buffer_b, sgn);             // FIR low pass filter: 6KHz cutoff
    00039     }
    00040     
    00041     sine_1KHz.reset();
    00042     for (float32_t *sgn=expected_output; sgn<(expected_output+TEST_LENGTH_SAMPLES); sgn += BLOCK_SIZE) {
    00043         sine_1KHz.generate(sgn);        // Generate a 1KHz sine wave
    00044     }
    00045     
    00046     float snr = arm_snr_f32(&expected_output[DELAY-1], &output[WARMUP-1], TEST_LENGTH_SAMPLES-WARMUP);
    00047     printf("snr: %f\n\r", snr);
    00048     if (snr < SNR_THRESHOLD_F32) {
    00049         printf("Failed\n\r");
    00050     } else {
    00051         printf("Success\n\r");
    00052     }
    00053 }