SciPy FIR Filter
This content relates to a deprecated version of Mbed
Mbed 2 is now deprecated. For the latest version please see the Mbed OS documentation.
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
FIR Filter Design¶
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}
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 }