Daiki Kato / mbed-os-lychee

mbed-os for GR-LYCHEE

Dependents:   mbed-os-example-blinky-gr-lychee GR-Boads_Camera_sample GR-Boards_Audio_Recoder GR-Boads_Camera_DisplayApp ... more

tools/dev/dsp_fir.py@0:f782d9c66c49, 2018-02-02 (annotated)

Committer:
dkato
Date:
Fri Feb 02 05:42:23 2018 +0000
Revision:
0:f782d9c66c49
mbed-os for GR-LYCHEE

Who changed what in which revision?

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