You are viewing an older revision! See the latest version

SciPy FIR Filter

Table of Contents

    FIR Filter Design with Python

    """
320 samples of (1000Hz + 15000 Hz) at 48 kHz
"""
from numpy import sin, arange, pi
from scipy.signal import lfilter, firwin
from pylab import figure, plot, grid, show

#------------------------------------------------
# Create a signal for demonstration.
#------------------------------------------------
sample_rate = 48000.
nsamples = 320

F_1KHz = 1000.
A_1KHz = 1.0

F_15KHz = 15000.
A_15KHz = 0.5

t = arange(nsamples) 
signal = A_1KHz * sin(2*pi*(F_1KHz/sample_rate)*t) + A_15KHz*sin(2*pi*(F_15KHz/sample_rate)*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)

figure(1)
# Plot the original signal
plot(t[:-warmup], signal[delay:-delay])

# Plot the filtered signal, shifted to compensate for the phase delay
plot(t[:-warmup], filtered_signal[warmup:], 'r-')

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)

    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