dsp_fir.py

00001 """
00002 mbed SDK
00003 Copyright (c) 2011-2013 ARM Limited
00004 
00005 Licensed under the Apache License, Version 2.0 (the "License");
00006 you may not use this file except in compliance with the License.
00007 You may obtain a copy of the License at
00008 
00009     http://www.apache.org/licenses/LICENSE-2.0
00010 
00011 Unless required by applicable law or agreed to in writing, software
00012 distributed under the License is distributed on an "AS IS" BASIS,
00013 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014 See the License for the specific language governing permissions and
00015 limitations under the License.
00016 """
00017 from numpy import sin, arange, pi
00018 from scipy.signal import lfilter, firwin
00019 from pylab import figure, plot, grid, show
00020 
00021 #------------------------------------------------
00022 # Create a signal for demonstration.
00023 #------------------------------------------------
00024 # 320 samples of (1000Hz + 15000 Hz) at 48 kHz
00025 sample_rate = 48000.
00026 nsamples = 320
00027 
00028 F_1KHz = 1000.
00029 A_1KHz = 1.0
00030 
00031 F_15KHz = 15000.
00032 A_15KHz = 0.5
00033 
00034 t = arange(nsamples) / sample_rate
00035 signal = A_1KHz * sin(2*pi*F_1KHz*t) + A_15KHz*sin(2*pi*F_15KHz*t)
00036 
00037 #------------------------------------------------
00038 # Create a FIR filter and apply it to signal.
00039 #------------------------------------------------
00040 # The Nyquist rate of the signal.
00041 nyq_rate = sample_rate / 2.
00042 
00043 # The cutoff frequency of the filter: 6KHz
00044 cutoff_hz = 6000.0
00045 
00046 # Length of the filter (number of coefficients, i.e. the filter order + 1)
00047 numtaps = 29
00048 
00049 # Use firwin to create a lowpass FIR filter
00050 fir_coeff = firwin(numtaps, cutoff_hz/nyq_rate)
00051 
00052 # Use lfilter to filter the signal with the FIR filter
00053 filtered_signal = lfilter(fir_coeff, 1.0, signal)
00054 
00055 #------------------------------------------------
00056 # Plot the original and filtered signals.
00057 #------------------------------------------------
00058 
00059 # The first N-1 samples are "corrupted" by the initial conditions
00060 warmup = numtaps - 1
00061 
00062 # The phase delay of the filtered signal
00063 delay = (warmup / 2) / sample_rate
00064 
00065 figure(1)
00066 # Plot the original signal
00067 plot(t, signal)
00068 
00069 # Plot the filtered signal, shifted to compensate for the phase delay
00070 plot(t-delay, filtered_signal, 'r-')
00071 
00072 # Plot just the "good" part of the filtered signal.  The first N-1
00073 # samples are "corrupted" by the initial conditions.
00074 plot(t[warmup:]-delay, filtered_signal[warmup:], 'g', linewidth=4)
00075 
00076 grid(True)
00077 
00078 show()
00079 
00080 #------------------------------------------------
00081 # Print values
00082 #------------------------------------------------
00083 def print_values(label, values):
00084     var = "float32_t %s[%d]" % (label, len(values))
00085     print "%-30s = {%s}" % (var, ', '.join(["%+.10f" % x for x in values]))
00086 
00087 print_values('signal', signal)
00088 print_values('fir_coeff', fir_coeff)
00089 print_values('filtered_signal', filtered_signal)