tool/mbed/mbed-sdk/workspace_tools/dev/dsp_fir.py

   1 """
   2 mbed SDK
   3 Copyright (c) 2011-2013 ARM Limited
   4
   5 Licensed under the Apache License, Version 2.0 (the "License");
   6 you may not use this file except in compliance with the License.
   7 You may obtain a copy of the License at
   8
   9     http://www.apache.org/licenses/LICENSE-2.0
  10
  11 Unless required by applicable law or agreed to in writing, software
  12 distributed under the License is distributed on an "AS IS" BASIS,
  13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14 See the License for the specific language governing permissions and
  15 limitations under the License.
  16 """
  17 from numpy import sin, arange, pi
  18 from scipy.signal import lfilter, firwin
  19 from pylab import figure, plot, grid, show
  20
  21 #------------------------------------------------
  22 # Create a signal for demonstration.
  23 #------------------------------------------------
  24 # 320 samples of (1000Hz + 15000 Hz) at 48 kHz
  25 sample_rate = 48000.
  26 nsamples = 320
  27
  28 F_1KHz = 1000.
  29 A_1KHz = 1.0
  30
  31 F_15KHz = 15000.
  32 A_15KHz = 0.5
  33
  34 t = arange(nsamples) / sample_rate
  35 signal = A_1KHz * sin(2*pi*F_1KHz*t) + A_15KHz*sin(2*pi*F_15KHz*t)
  36
  37 #------------------------------------------------
  38 # Create a FIR filter and apply it to signal.
  39 #------------------------------------------------
  40 # The Nyquist rate of the signal.
  41 nyq_rate = sample_rate / 2.
  42
  43 # The cutoff frequency of the filter: 6KHz
  44 cutoff_hz = 6000.0
  45
  46 # Length of the filter (number of coefficients, i.e. the filter order + 1)
  47 numtaps = 29
  48
  49 # Use firwin to create a lowpass FIR filter
  50 fir_coeff = firwin(numtaps, cutoff_hz/nyq_rate)
  51
  52 # Use lfilter to filter the signal with the FIR filter
  53 filtered_signal = lfilter(fir_coeff, 1.0, signal)
  54
  55 #------------------------------------------------
  56 # Plot the original and filtered signals.
  57 #------------------------------------------------
  58
  59 # The first N-1 samples are "corrupted" by the initial conditions
  60 warmup = numtaps - 1
  61
  62 # The phase delay of the filtered signal
  63 delay = (warmup / 2) / sample_rate
  64
  65 figure(1)
  66 # Plot the original signal
  67 plot(t, signal)
  68
  69 # Plot the filtered signal, shifted to compensate for the phase delay
  70 plot(t-delay, filtered_signal, 'r-')
  71
  72 # Plot just the "good" part of the filtered signal.  The first N-1
  73 # samples are "corrupted" by the initial conditions.
  74 plot(t[warmup:]-delay, filtered_signal[warmup:], 'g', linewidth=4)
  75
  76 grid(True)
  77
  78 show()
  79
  80 #------------------------------------------------
  81 # Print values
  82 #------------------------------------------------
  83 def print_values(label, values):
  84     var = "float32_t %s[%d]" % (label, len(values))
  85     print "%-30s = {%s}" % (var, ', '.join(["%+.10f" % x for x in values]))
  86
  87 print_values('signal', signal)
  88 print_values('fir_coeff', fir_coeff)
  89 print_values('filtered_signal', filtered_signal)