FIR Filter

Cookbook » FIR Filter

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 library mbed-dsp

CMSIS DSP library

CMSIS DSP API¶

Import program
cmsis_dsp_fir - main.cpp

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
dsp_fir - main.cpp

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 }

All wikipages

FIR Filter

Table of Contents

FIR Filter Design with Python¶

Python FIR Filter Design

coefficients

mbed-dsp library¶

Import library mbed-dsp

CMSIS DSP API¶

Import program
cmsis_dsp_fir - main.cpp

mbed DSP API¶

Import program
dsp_fir - main.cpp

Recent changes

Important Information for this Arm website

FIR Filter

Table of Contents

FIR Filter Design with Python¶

Python FIR Filter Design

coefficients

mbed-dsp library¶

Import librarymbed-dsp

CMSIS DSP API¶

Import programcmsis_dsp_fir - main.cpp

mbed DSP API¶

Import programdsp_fir - main.cpp

Recent changes

Important Information for this Arm website

Access Warning

Import library mbed-dsp

Import program
cmsis_dsp_fir - main.cpp

Import program
dsp_fir - main.cpp