dtmf.py

#!/usr/bin/python

# this python file was written by luke trusheim
# distribution is unrestricted provided it is without charge and includes attribution

import sys
import time

import base64
import random as random

import datetime
import time
import math

import matplotlib.pyplot as plt
import numpy as np

from cpe367_wav import cpe367_wav
from cpe367_sig_analyzer import cpe367_sig_analyzer
from my_fifo import my_fifo


bp_out_g = 0
abs_out_g = 0
lp_out_g = 0
detected_g = 0

scale_factor = 128


############################################
############################################
# define routine for detecting DTMF tones
def process_wav(fpath_sig_in):
	
		
	###############################
	# define list of signals to be displayed by and analyzer
	#  note that the signal analyzer already includes: 'symbol_val','symbol_det','error'
	
	current_analyzed_signal = '852'
	bpass = current_analyzed_signal + '_bpass'
	abs = current_analyzed_signal + '_abs'
	lp = current_analyzed_signal + '_lp'
	detected = current_analyzed_signal + '_detected'

	more_sig_list = ['sig_1', bpass, abs, lp, detected, '1336_lp', '1336_detected', '770_lp', '770_detected']
	
	# sample rate is 4kHz
	fs = 4000
	
	# instantiate signal analyzer and load data
	s2 = cpe367_sig_analyzer(more_sig_list,fs)
	s2.load(fpath_sig_in)
	s2.print_desc()
	
	########################
	# students: setup filters

	keypad = [
    ['1', '2', '3', 'A'],
    ['4', '5', '6', 'B'],
    ['7', '8', '9', 'C'],
    ['*', '0', '#', 'D']
	]

	# gains = [1/14.295, 1/13.59, 1/13.097, 1/12.81, 1/13.435, 1/14.74, 1/17.409, 1/23.38]
	bp_ak = {"697" : [], "770": [], "852" : [], "941": [], "1209" : [], "1336": [], "1477" : [], "1633": []}
	bp_bk = {"697" : [], "770": [], "852" : [], "941": [], "1209" : [], "1336": [], "1477" : [], "1633": []}
	bp_ak["697"], bp_bk["697"] = bp_filter_coeff(697, fs, 0.92)
	bp_ak["770"], bp_bk["770"] = bp_filter_coeff(770, fs, 0.90)
	bp_ak["852"], bp_bk["852"] = bp_filter_coeff(852, fs, 0.91)
	bp_ak["941"], bp_bk["941"] = bp_filter_coeff(941, fs, 0.92)
	bp_ak["1209"], bp_bk["1209"] = bp_filter_coeff(1209, fs, 0.92)
	bp_ak["1336"], bp_bk["1336"] = bp_filter_coeff(1336, fs, 0.86)
	bp_ak["1477"], bp_bk["1477"] = bp_filter_coeff(1477, fs, 0.87)
	bp_ak["1633"], bp_bk["1633"] = bp_filter_coeff(1633, fs, 0.87)

	lp_bk = {"697" : [], "770": [], "852" : [], "941": [], "1209" : [], "1336": [], "1477" : [], "1633": []}
	lp_ak = {"697" : [], "770": [], "852" : [], "941": [], "1209" : [], "1336": [], "1477" : [], "1633": []}
	lp_ak["697"], lp_bk["697"] = lp_filter_coeff(697, fs)
	lp_ak["770"], lp_bk["770"] = lp_filter_coeff(770, fs)
	lp_ak["852"], lp_bk["852"] = lp_filter_coeff(852, fs)
	lp_ak["941"], lp_bk["941"] = lp_filter_coeff(941, fs)
	lp_ak["1209"], lp_bk["1209"] = lp_filter_coeff(1209, fs)
	lp_ak["1336"], lp_bk["1336"] = lp_filter_coeff(1336, fs)
	lp_ak["1477"], lp_bk["1477"] = lp_filter_coeff(1477, fs)
	lp_ak["1633"], lp_bk["1633"] = lp_filter_coeff(1633, fs)

	thresholds = {"697" : 28, "770": 30, "852" : 28, "941": 27, "1209" : 29, "1336": 28, "1477" : 28, "1633": 28}

	M = 3

	fifos = {}
	
	# process input	
	xin = 0
	for n_curr in range(s2.get_len()):
	
		# read next input sample from the signal analyzer
		xin = s2.get('xin',n_curr)

		detected_697_Hz = detect_freq(xin, M, thresholds, "697", bp_bk, bp_ak, lp_bk, lp_ak, fifos)
		
		

		# detected_770_Hz = detected_852_Hz = detected_941_Hz = detected_1209_Hz = detected_1336_Hz = detected_1477_Hz = detected_1633_Hz = 0
		detected_770_Hz = detect_freq(xin, M, thresholds, "770", bp_bk, bp_ak, lp_bk, lp_ak, fifos)

		s2.set('770_lp', n_curr, lp_out_g)
		s2.set('770_detected', n_curr, detected_g)

		detected_852_Hz = detect_freq(xin, M, thresholds, "852", bp_bk, bp_ak, lp_bk, lp_ak, fifos)

		# s2.set(bpass, n_curr, bp_out_g)
		# s2.set(abs, n_curr, abs_out_g)
		# s2.set(lp, n_curr, lp_out_g)
		# s2.set(detected, n_curr, detected_g)

		detected_941_Hz = detect_freq(xin, M, thresholds, "941", bp_bk, bp_ak, lp_bk, lp_ak, fifos)

		detected_1209_Hz = detect_freq(xin, M, thresholds, "1209", bp_bk, bp_ak, lp_bk, lp_ak, fifos)

		detected_1336_Hz = detect_freq(xin, M, thresholds, "1336", bp_bk, bp_ak, lp_bk, lp_ak, fifos)

		s2.set('1336_lp', n_curr, lp_out_g)
		s2.set('1336_detected', n_curr, detected_g)
		
		detected_1477_Hz = detect_freq(xin, M, thresholds, "1477", bp_bk, bp_ak, lp_bk, lp_ak, fifos)

		detected_1633_Hz = detect_freq(xin, M, thresholds, "1633", bp_bk, bp_ak, lp_bk, lp_ak, fifos)

		detected_symbol = decode(detected_697_Hz, detected_770_Hz, detected_852_Hz, detected_941_Hz, detected_1209_Hz, detected_1336_Hz, detected_1477_Hz, detected_1633_Hz, keypad)
		
		########################
		# students: evaluate each filter and implement other processing blocks

		
		########################
		# students: combine results from filtering stages
		#  and find (best guess of) symbol that is present at this sample time

		try:
			symbol_val_det = int(detected_symbol)
		except:
			symbol_val_det = 10

		# save intermediate signals as needed, for plotting
		#  add signals, as desired!
		s2.set('sig_1',n_curr,xin)

		# save detected symbol
		s2.set('symbol_det',n_curr,symbol_val_det)

		# get correct symbol (provided within the signal analyzer)
		symbol_val = s2.get('symbol_val',n_curr)

		# compare detected signal to correct signal
		symbol_val_err = 0
		if symbol_val != symbol_val_det: symbol_val_err = 1
		
		# save error signal
		s2.set('error',n_curr,symbol_val_err)
		
	
	# display mean of error signal
	err_mean = s2.get_mean('error')
	print('mean error = '+str( round(100 * err_mean,1) )+'%')
		
	# define which signals should be plotted
	# plot_sig_list = ['sig_1', bpass, abs, lp, detected, 'symbol_val','symbol_det','error']
	plot_sig_list = ['1336_lp', '1336_detected', '770_lp', '770_detected', 'symbol_val','symbol_det', 'error']
	
	# plot results
	s2.plot(plot_sig_list)
	
	return True

def detect_freq(xin, M, thresholds, target_freq, bp_bk, bp_ak, lp_bk, lp_ak, fifos):
	global bp_out_g
	global abs_out_g
	global lp_out_g
	global detected_g

	if target_freq not in fifos:
		fifos[target_freq] = {
            'inputs': my_fifo(M),
            'bp_outputs': my_fifo(M),
            'abs_outputs': my_fifo(M),
            'lp_outputs': my_fifo(M),
        }

	inputs 		= fifos[target_freq]["inputs"]
	bp_outputs 	= fifos[target_freq]["bp_outputs"]
	abs_outputs = fifos[target_freq]["abs_outputs"]
	lp_outputs 	= fifos[target_freq]["lp_outputs"]

	inputs.update(xin)
		
	bp_out = bp_bk[target_freq][0] * inputs.get(0) - bp_ak[target_freq][1] * bp_outputs.get(0) - bp_ak[target_freq][2] * bp_outputs.get(1)
	bp_out /= scale_factor
	bp_outputs.update(bp_out)
	bp_out_g = bp_out

	abs_out = abs(bp_outputs.get(0))
	abs_outputs.update(abs_out)
	abs_out_g = abs_out

	lp_out = (lp_bk[target_freq][0] * abs_outputs.get(0) + lp_bk[target_freq][1] * abs_outputs.get(1) + lp_bk[target_freq][2] * abs_outputs.get(2) 
		   - lp_ak[target_freq][1] * lp_outputs.get(0) - lp_ak[target_freq][2] * lp_outputs.get(1))
	lp_out /= scale_factor
	lp_outputs.update(lp_out)
	lp_out_g = lp_out

	detected = 1 if lp_outputs.get(0) > thresholds[target_freq] else 0
	detected_g = detected

	return detected



def bp_filter_coeff(f0, fs, pole_radius):
	F = f0/fs
	w0 = 2 * math.pi * F

	zeros = [0, 0]
	bk = np.poly(zeros)

	poles = np.array([pole_radius * math.e ** (1j * w0), pole_radius * math.e ** (-1j * w0)])
	ak = np.poly(poles)

	z = (math.e ** (1j * w0))		
	G = abs((z ** 2) / (ak[0] * (z ** 2) + ak[1] * z + ak[2])) # calculate gain
	bk = 1/G * bk

	bk *= scale_factor
	ak *= scale_factor

	if (scale_factor != 1):
		bk = np.round(bk, 0)
		ak = np.round(ak, 0)

	return ak, bk
       

def lp_filter_coeff(fc, fs):
	Fc = fc/fs
	w0 = 2 * math.pi * Fc

	zeros = np.array([math.e ** (1j * w0), math.e ** (-1j * w0)])
	bk = np.poly(zeros)
	pole_radius = 0.63

	poles = [pole_radius, pole_radius]
	ak = np.poly(poles)

	z = 1		
	G = abs((bk[0] * (z ** 2) + bk[1] * z + bk[2]) / (ak[0] * (z ** 2) + ak[1] * z + ak[2])) # calculate gain
	bk = 1/G * bk

	bk *= scale_factor
	ak *= scale_factor

	if (scale_factor != 1):
		bk = np.round(bk, 0)
		ak = np.round(ak, 0)

	return ak, bk   

def decode(detected_697_Hz, detected_770_Hz, detected_852_Hz, detected_941_Hz, detected_1209_Hz, detected_1336_Hz, detected_1477_Hz, detected_1633_Hz, keypad):
	row_index = column_index = 0
	
	if (detected_697_Hz):
		row_index = 0
	elif (detected_770_Hz):
		row_index = 1
	elif (detected_852_Hz):
		row_index = 2
	elif (detected_941_Hz):
		row_index = 3

	if (detected_1209_Hz):
		column_index = 0
	elif (detected_1336_Hz):
		column_index = 1
	elif (detected_1477_Hz):
		column_index = 2
	elif (detected_1633_Hz):
		column_index = 3

	return keypad[row_index][column_index]

	
############################################
############################################
# define main program
def main():

	# check python version!
	major_version = int(sys.version[0])
	if major_version < 3:
		print('Sorry! must be run using python3.')
		print('Current version: ')
		print(sys.version)
		return False
		
	# assign file name
	# fpath_sig_in = 'dtmf_signals_slow.txt'
	fpath_sig_in = 'dtmf_signals_fast.txt'
	
	
	# let's do it!
	return process_wav(fpath_sig_in)

	
############################################
############################################
# call main function
if __name__ == '__main__':
	
	main()
	quit()