detect.cpp

GNU ccAudio2 is a stand-alone portable C++ class framework for manipulating audio data. It has exist

//
//   DTMF Receiver module, part of:
//      BSD Telephony Of Mexico "Zapata" Telecom Library, version 1.10  12/9/01
//
//   Part of the "Zapata" Computer Telephony Technology.
//
//   See http://www.bsdtelephony.com.mx
//
//  The technologies, software, hardware, designs, drawings, scheumatics, board
//  layouts and/or artwork, concepts, methodologies (including the use of all
//  of these, and that which is derived from the use of all of these), all other
//  intellectual properties contained herein, and all intellectual property
//  rights have been and shall continue to be expressly for the benefit of all
//  mankind, and are perpetually placed in the public domain, and may be used,
//  copied, and/or modified by anyone, in any manner, for any legal purpose,
//  without restriction.
//
// tone_detect.c - General telephony tone detection, and specific
//                      detection of DTMF.
//
//        Copyright (C) 2001  Steve Underwood <[e-mail]>
//
//        Despite my general liking of the GPL, I place this code in the
//        public domain for the benefit of all mankind - even the slimy
//        ones who might try to proprietize my work and use it to my
//        detriment.
//

#include "private.h"
#include <cstring>
#include <cstdio>
#include <cctype>
#include <cmath>
#include <ctime>

#ifdef  HAVE_STDINT_H
#include <stdint.h>
#endif

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#include "audio2.h"


/* Basic DTMF specs:
 *
 * Minimum tone on = 40ms
 * Minimum tone off = 50ms
 * Maximum digit rate = 10 per second
 * Normal twist <= 8dB accepted
 * Reverse twist <= 4dB accepted
 * S/N >= 15dB will detect OK
 * Attenuation <= 26dB will detect OK
 * Frequency tolerance +- 1.5% will detect, +-3.5% will reject
 */

using namespace ost;

#define SAMPLE_RATE		8000.0

#define DTMF_THRESHOLD		8.0e7
#define FAX_THRESHOLD		8.0e7
#define FAX_2ND_HARMONIC	2.0     /* 4dB */
#define DTMF_NORMAL_TWIST	8.0     /* 8dB */
#define DTMF_REVERSE_TWIST	4.0     /* 4dB normal */
#define DTMF_RELATIVE_PEAK_ROW	6.3     /* 8dB */
#define DTMF_RELATIVE_PEAK_COL	6.3     /* 8dB */
#define DTMF_2ND_HARMONIC_ROW	2.5     /* 4dB normal */
#define DTMF_2ND_HARMONIC_COL	63.1    /* 18dB */

DTMFDetect::DTMFDetect()
{
	int i;
	float theta;
	static float dtmf_row[] = { 697.0, 770.0, 852.0, 941.0 };
	static float dtmf_col[] = { 1209.0, 1336.0, 1477.0, 1633.0 };
	static float fax_freq = 1100.0;

	state = (dtmf_detect_state_t *)malloc(sizeof(dtmf_detect_state_t));
	memset(state, 0, sizeof(state));

	for(i = 0; i < 4; i++)
	{
		theta = (float)(2.0 * M_PI * (dtmf_row[i] / SAMPLE_RATE));
		dtmf_detect_row[i].fac = (float)(2.0 * cos(theta));

		theta = (float)(2.0 * M_PI * (dtmf_col[i] / SAMPLE_RATE));
		dtmf_detect_col[i].fac = (float)(2.0 * cos(theta));

		theta = (float)(2.0 * M_PI * (dtmf_row[i] * 2.0 / SAMPLE_RATE));
		dtmf_detect_row_2nd[i].fac = (float)(2.0 * cos(theta));

		theta = (float)(2.0 * M_PI * (dtmf_col[i] * 2.0 / SAMPLE_RATE));
		dtmf_detect_col_2nd[i].fac = (float)(2.0 * cos(theta));

		goertzelInit(&state->row_out[i], &dtmf_detect_row[i]);
		goertzelInit(&state->col_out[i], &dtmf_detect_col[i]);
		goertzelInit(&state->row_out2nd[i], &dtmf_detect_row_2nd[i]);
		goertzelInit(&state->col_out2nd[i], &dtmf_detect_col_2nd[i]);

		state->energy = 0.0;
	}

	// Same for the fax detector
	theta = (float)(2.0 * M_PI * (fax_freq / SAMPLE_RATE));
	fax_detect.fac = (float)(2.0 * cos(theta));
	goertzelInit(&state->fax_tone, &fax_detect);

	// Same for the fax detector 2nd harmonic
	theta = (float)(2.0 * M_PI * (fax_freq / SAMPLE_RATE));
	fax_detect_2nd.fac = (float)(2.0 * cos(theta));
	goertzelInit(&state->fax_tone2nd, &fax_detect_2nd);

	state->current_digits = 0;
	state->current_sample = 0;
	state->detected_digits = 0;
	state->lost_digits = 0;
	state->digits[0] = '\0';
	state->mhit = 0;
}

DTMFDetect::~DTMFDetect()
{
	if(state)
	{
		free(state);
		state = NULL;
	}
	return;
}

void DTMFDetect::goertzelInit(goertzel_state_t *s, tone_detection_descriptor_t *t)
{
	s->v2 = s->v3 = 0.0;
	s->fac = t->fac;
}

void DTMFDetect::goertzelUpdate(goertzel_state_t *s,
                     Sample x[],
                     int samples)
{
	int i;
	float v1;

	for (i = 0;  i < samples;  i++)
	{
		v1 = s->v2;
		s->v2 = s->v3;
		s->v3 = s->fac*s->v2 - v1 + x[i];
	}
}

float DTMFDetect::goertzelResult (goertzel_state_t *s)
{
	return s->v3 * s->v3 + s->v2 * s->v2 - s->v2 *s->v3 *s->fac;
}

int DTMFDetect::putSamples(Linear amp, int samples)
{
	static char dtmf_positions[] = "123A" "456B" "789C" "*0#D";
	float row_energy[4];
	float col_energy[4];
	float fax_energy;
	float fax_energy_2nd;
	float famp;
	float v1;
	int i;
	int j;
	int sample;
	int best_row;
	int best_col;
	int hit;
	int limit;

	hit = 0;
	for (sample = 0;  sample < samples;  sample = limit)
	{
		// 102 is optimised to meet the DTMF specs.
		if ((samples - sample) >= (102 - state->current_sample))
			limit = sample + (102 - state->current_sample);
		else
			limit = samples;

		// The following unrolled loop takes only 35% (rough estimate) of the 
		// time of a rolled loop on the machine on which it was developed 
		for(j = sample;  j < limit;  j++)
		{
			famp = amp[j];
	  		state->energy += famp*famp;
	    
			// With GCC 2.95, the following unrolled code seems to take about 35%
			// (rough estimate) as long as a neat little 0-3 loop
			v1 = state->row_out[0].v2;
			state->row_out[0].v2 = state->row_out[0].v3;
			state->row_out[0].v3 = state->row_out[0].fac*state->row_out[0].v2 - v1 + famp;

			v1 = state->col_out[0].v2;
			state->col_out[0].v2 = state->col_out[0].v3;
			state->col_out[0].v3 = state->col_out[0].fac*state->col_out[0].v2 - v1 + famp;

			v1 = state->row_out[1].v2;
			state->row_out[1].v2 = state->row_out[1].v3;
			state->row_out[1].v3 = state->row_out[1].fac*state->row_out[1].v2 - v1 + famp;

			v1 = state->col_out[1].v2;
			state->col_out[1].v2 = state->col_out[1].v3;
			state->col_out[1].v3 = state->col_out[1].fac*state->col_out[1].v2 - v1 + famp;

			v1 = state->row_out[2].v2;
			state->row_out[2].v2 = state->row_out[2].v3;
			state->row_out[2].v3 = state->row_out[2].fac*state->row_out[2].v2 - v1 + famp;

			v1 = state->col_out[2].v2;
			state->col_out[2].v2 = state->col_out[2].v3;
			state->col_out[2].v3 = state->col_out[2].fac*state->col_out[2].v2 - v1 + famp;

			v1 = state->row_out[3].v2;
			state->row_out[3].v2 = state->row_out[3].v3;
			state->row_out[3].v3 = state->row_out[3].fac*state->row_out[3].v2 - v1 + famp;

			v1 = state->col_out[3].v2;
			state->col_out[3].v2 = state->col_out[3].v3;
			state->col_out[3].v3 = state->col_out[3].fac*state->col_out[3].v2 - v1 + famp;

			v1 = state->col_out2nd[0].v2;
			state->col_out2nd[0].v2 = state->col_out2nd[0].v3;
			state->col_out2nd[0].v3 = state->col_out2nd[0].fac*state->col_out2nd[0].v2 - v1 + famp;

			v1 = state->row_out2nd[0].v2;
			state->row_out2nd[0].v2 = state->row_out2nd[0].v3;
			state->row_out2nd[0].v3 = state->row_out2nd[0].fac*state->row_out2nd[0].v2 - v1 + famp;

			v1 = state->col_out2nd[1].v2;
			state->col_out2nd[1].v2 = state->col_out2nd[1].v3;
			state->col_out2nd[1].v3 = state->col_out2nd[1].fac*state->col_out2nd[1].v2 - v1 + famp;

			v1 = state->row_out2nd[1].v2;
			state->row_out2nd[1].v2 = state->row_out2nd[1].v3;
			state->row_out2nd[1].v3 = state->row_out2nd[1].fac*state->row_out2nd[1].v2 - v1 + famp;

			v1 = state->col_out2nd[2].v2;
			state->col_out2nd[2].v2 = state->col_out2nd[2].v3;
			state->col_out2nd[2].v3 = state->col_out2nd[2].fac*state->col_out2nd[2].v2 - v1 + famp;

			v1 = state->row_out2nd[2].v2;
			state->row_out2nd[2].v2 = state->row_out2nd[2].v3;
			state->row_out2nd[2].v3 = state->row_out2nd[2].fac*state->row_out2nd[2].v2 - v1 + famp;

			v1 = state->col_out2nd[3].v2;
			state->col_out2nd[3].v2 = state->col_out2nd[3].v3;
			state->col_out2nd[3].v3 = state->col_out2nd[3].fac*state->col_out2nd[3].v2 - v1 + famp;

			v1 = state->row_out2nd[3].v2;
			state->row_out2nd[3].v2 = state->row_out2nd[3].v3;
			state->row_out2nd[3].v3 = state->row_out2nd[3].fac*state->row_out2nd[3].v2 - v1 + famp;

			v1 = state->fax_tone.v2;
			state->fax_tone.v2 = state->fax_tone.v3;
			state->fax_tone.v3 = state->fax_tone.fac*state->fax_tone.v2 - v1 + famp;

			v1 = state->fax_tone.v2;
			state->fax_tone2nd.v2 = state->fax_tone2nd.v3;
			state->fax_tone2nd.v3 = state->fax_tone2nd.fac*state->fax_tone2nd.v2 - v1 + famp;
		}
		state->current_sample += (limit - sample);
		if(state->current_sample < 102)
			continue;

		fax_energy = goertzelResult(&state->fax_tone);

		// We are at the end of a DTMF detection block
		// Find the peak row and the peak column
		row_energy[0] = goertzelResult (&state->row_out[0]);
		col_energy[0] = goertzelResult (&state->col_out[0]);

		for(best_row = best_col = 0, i = 1;  i < 4;  i++)
		{
			row_energy[i] = goertzelResult (&state->row_out[i]);
			if(row_energy[i] > row_energy[best_row])
				 best_row = i;
			col_energy[i] = goertzelResult (&state->col_out[i]);
			if(col_energy[i] > col_energy[best_col])
				best_col = i;
		}
		hit = 0;

		// Basic signal level test and the twist test
		if(row_energy[best_row] >= DTMF_THRESHOLD &&
			col_energy[best_col] >= DTMF_THRESHOLD &&
			col_energy[best_col] < row_energy[best_row] * DTMF_REVERSE_TWIST &&
			col_energy[best_col] * DTMF_NORMAL_TWIST > row_energy[best_row])
		{
			// Relative peak test
			for(i = 0;  i < 4;  i++)
			{
				if ((i != best_col && 
					col_energy[i]*DTMF_RELATIVE_PEAK_COL > col_energy[best_col]) ||
					(i != best_row && row_energy[i]*DTMF_RELATIVE_PEAK_ROW > row_energy[best_row]))
					break;
			}
			// ... and second harmonic test
			if(i >= 4 &&
				(row_energy[best_row] + col_energy[best_col]) > 42.0*state->energy &&
				goertzelResult (&state->col_out2nd[best_col])*DTMF_2ND_HARMONIC_COL < col_energy[best_col] &&
				goertzelResult (&state->row_out2nd[best_row])*DTMF_2ND_HARMONIC_ROW < row_energy[best_row])
			{
				hit = dtmf_positions[(best_row << 2) + best_col];
				// Look for two successive similar results
				// The logic in the next test is:
				//   We need two successive identical clean detects, with
				//   something different preceeding it. This can work with
				//   back to back differing digits. More importantly, it
				//   can work with nasty phones that give a very wobbly start
				//   to a digit.
				if (hit == state->hit3  &&  state->hit3 != state->hit2)
				{
					state->mhit = hit;
					state->digit_hits[(best_row << 2) + best_col]++;
					state->detected_digits++;
					if (state->current_digits < 128)
					{
						state->digits[state->current_digits++] = hit;
						state->digits[state->current_digits] = '\0';
					}
					else
					{
						state->lost_digits++;
					}
				}
			}
		}

		if (!hit && (fax_energy >= FAX_THRESHOLD) && (fax_energy > state->energy * 21.0))
		{
			fax_energy_2nd = goertzelResult(&state->fax_tone2nd);
			if (fax_energy_2nd * FAX_2ND_HARMONIC < fax_energy)
			{
				// XXX Probably need better checking than just this the energy
				hit = 'f';
				state->fax_hits++;
			} /* Don't reset fax hits counter */
		} else {
			if (state->fax_hits > 5)
			{
				state->mhit = 'f';
				state->detected_digits++;
				if (state->current_digits < 128)
				{
					state->digits[state->current_digits++] = hit;
					state->digits[state->current_digits] = '\0';
				}
				else
					state->lost_digits++;
			}
			state->fax_hits = 0;
		}
		state->hit1 = state->hit2;
		state->hit2 = state->hit3;
		state->hit3 = hit;
		// Reinitialise the detector for the next block
		for (i = 0;  i < 4;  i++)
		{
			goertzelInit (&state->row_out[i], &dtmf_detect_row[i]);
			goertzelInit (&state->col_out[i], &dtmf_detect_col[i]);
			goertzelInit (&state->row_out2nd[i], &dtmf_detect_row_2nd[i]);
			goertzelInit (&state->col_out2nd[i], &dtmf_detect_col_2nd[i]);
		}
		goertzelInit (&state->fax_tone, &fax_detect);
		goertzelInit (&state->fax_tone2nd, &fax_detect_2nd);
		state->energy = 0.0;
		state->current_sample = 0;
	}
	if ((!state->mhit) || (state->mhit != hit))
	{
		state->mhit = 0;
		return(0);
	}
	return (hit);
}

int DTMFDetect::getResult(char *buf, int max)
{
	if (max > state->current_digits)
		max = state->current_digits;
	if (max > 0)
	{
		memcpy (buf, state->digits, max);
		memmove (state->digits, state->digits + max, state->current_digits - max);
		state->current_digits -= max;
	}
	buf[max] = '\0';
	return  max;
}