Diff: dev/dsp_fir.py

Revision:

0:4a2e5f0422d6

diff -r 000000000000 -r 4a2e5f0422d6 dev/dsp_fir.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dev/dsp_fir.py	Wed Mar 30 16:23:44 2016 +0000
@@ -0,0 +1,89 @@
+"""
+mbed SDK
+Copyright (c) 2011-2013 ARM Limited
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+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)