## Tone Detection Using the Teager Kaiser Energy Operator (Part 1)

The Teager Kaiser Energy Operator (TKEO) is a non-linear operator that can provide an estimate of the instantaneous frequency and amplitude of a signal with AM and FM modulation. Here I will show how it can be used to measure the frequency of a non-modulated sinusoid or tone. The TKEO is defined as follows for discrete-time:

$\Psi [x(n)] = {x(n)}^{2}-x(n+1)x(n-1)$

Applying the TKEO to a sinusoid, one can derive this equation:

$\Psi (x(n)) = {A}^{2}{sin(\Omega )}^{2}$

where $\Omega = \frac{2\pi f}{{F}_{S}}$ is the normalized angular frequency

The Discrete Energy Separation Algorithms known as DESA-1 and DESA-2 use the TKEO on a signal and its derivative, and are used to estimate the frequency of a sinusoid. For a sinusoid without noise, they actually give an exact measurement. The equations for each are as follows:

DESA-1: $\left(sin{\left(\frac{\Omega }{2}\right)}\right)^{2}=\frac{\Psi ({x}_{n}-{x}_{n-1})+\Psi ({x}_{n+1}-{x}_{n})}{8\Psi {x}_{n}}$

DESA-2: ${\left(sin(\Omega )\right)}^{2}=\frac{\Psi ({x}_{n+1}-{x}_{n-1})}{4\Psi {x}_{n}}$

It is straightforward to solve these equations for the frequency f. The advantage of these expressions over the TKEO on its own is that the frequency is separated out from the amplitude. There are also DESA-1 and DESA-2 equations (not shown here) for deriving the instantaneous amplitude of a sinusoid (useful for finding the amplitude of amplitude modulated signals).

I have written some example C code to demonstrate the performance of the DESA algorithms for estimating the frequency of tones. The code reads in a block of samples, calculates the frequency estimates for each new sample, and then finds the mean and variance of the estimates. The test program prints out the mean frequency, the variance and the standard deviation for every block of samples. For the input file, use a headerless 16 bit PCM file sampled at 8000 Hz.

The variables x1, x2, x3 are the delayed inputs, that is, x[n-1], x[n-2] and x[n-3]. I have shifted everything in time to make the equations causal. That is, there are no future samples such as x[n+1]. The variables diff0 to diff3 are the current and delayed differences, that is diff0 is x[n] – x[n-1] and diff1 is diff0 delayed by one sample.


/*
* desa.c
*
* Tone detection using the Teager Kaiser Energy Operator to
* measure frequency (Desa-1 and Desa-2 algorithms)
* (DESA = Discrete Energy Separation Algorithm)
*
*/

#include <stdio.h>
#include <stdint.h>
#include <memory.h>
#include <math.h>

#define BLOCK   40          // samples processed each invocation
#define SAMPLE_RATE 8000.0  // sample rate of input signal

/*
* Purpose: detect a tone using DESA-1 algorithm
*
* Parameters:
*      input       pointer to input samples
*      variance    the variance of the frequency estimates
*
* Return value: frequency estimate in Hz
*/
double desa1( double *input, double *variance )
{
// detector variables
static double diff0 = 0.0;  // delayed differences
static double diff1 = 0.0;
static double diff2 = 0.0;
static double diff3 = 0.0;
static double x1 = 0.0;     // delayed inputs
static double x2 = 0.0;
static double x3 = 0.0;

double num; // numerator
double den; // denominator
double freq[BLOCK]; // frequency estimates
int i;

// calculate the frequency estimate for each sample
for ( i = 0; i < BLOCK; i++ )
{
diff0 = input[i] - x1;
num = diff2 * diff2 - diff1 * diff3
+ diff1 * diff1 - diff0 * diff2;
den = x2 * x2 - x1 * x3;
freq[i] = SAMPLE_RATE * asin(sqrt(num/(8.0 * den))) / M_PI;
// handle errors - division by zero, square root of
// negative number or asin of number > 1 or < -1
if (isnan(freq[i]))
{
freq[i] = 0;
}

diff3 = diff2;
diff2 = diff1;
diff1 = diff0;
x3 = x2;
x2 = x1;
x1 = input[i];
}

// calculate mean frequency
double mean = 0.0;
for ( i = 0; i < BLOCK; i++ )
{
mean += freq[i];
}
mean /= (double)BLOCK;

// calculate the variance in the frequency estimates
*variance = 0.0;
for ( i = 0; i < BLOCK; i++ )
{
*variance += freq[i] * freq[i];
}
*variance /= (double)BLOCK;
*variance -= (mean * mean);

return mean;
}

/*
* Purpose: detect a tone using DESA-2 algorithm
*
* Parameters:
*      input       pointer to input samples
*      variance    the variance of the frequency estimates
*
* Return value: frequency estimate in Hz
*/
double desa2( double *input, double *variance )
{
// detector variables
static double diff0 = 0.0;  // delayed differences
static double diff1 = 0.0;
static double diff2 = 0.0;
static double x1 = 0.0;     // delayed inputs
static double x2 = 0.0;
static double x3 = 0.0;

double num; // numerator
double den; // denominator
double freq[BLOCK]; // frequency estimates
int i;

// calculate the frequency estimate for each sample
for ( i = 0; i < BLOCK; i++ )
{
diff0 = input[i] - x2;
num = diff1 * diff1 - diff0 * diff2;
den = x2 * x2 - x1 * x3;
freq[i] = SAMPLE_RATE * asin(sqrt(num/(4.0 * den)))
/ (2.0 * M_PI);
// handle errors - division by zero, square root of
// negative number or asin of number > 1 or < -1
if (isnan(freq[i]))
{
freq[i] = 0;
}

diff2 = diff1;
diff1 = diff0;
x3 = x2;
x2 = x1;
x1 = input[i];
}

// calculate mean frequency
double mean = 0.0;
for ( i = 0; i < BLOCK; i++ )
{
mean += freq[i];
}
mean /= (double)BLOCK;

// calculate the variance in the frequency estimates
*variance = 0.0;
for ( i = 0; i < BLOCK; i++ )
{
*variance += freq[i] * freq[i];
}
*variance /= (double)BLOCK;
*variance -= (mean * mean);

return mean;
}

// convert 16 bit ints to doubles
void intToFloat( int16_t *input, double *output, int length )
{
int i;

for ( i = 0; i < length; i++ ) {
output[i] = (double)input[i];
}
}

int main( int argc, char *argv[] )
{
int16_t intData[BLOCK];
double inputData[BLOCK];
double frequency;
double variance;
int numWords;
int sampleCount;
char *inFileName;
FILE *inFile;

inFileName = NULL;

if ( argc == 2 )
{
inFileName = argv[1];
printf("input file name = %s\n",inFileName);
}
else
{
printf("Incorrect arguments, usage: desa filename\n");
return(1);
}

inFile = fopen(inFileName,"rb");
if ( !inFile )
{
printf("Exiting.Cannot open input file %s\n",inFileName);
fclose( inFile );
return(1);
}

// start counting frames
sampleCount = 0;

numWords = fread(intData, sizeof(int16_t), BLOCK, inFile );

// until end of file
while( numWords == BLOCK )
{
intToFloat( intData, inputData, numWords );
printf("sample = %d\n",sampleCount);

// get the frequency estimates
frequency = desa1( inputData, &variance );
printf("Desa1: Mean freq = %f, var = %f, std dev = %f\n",
frequency, variance, sqrt(variance));

frequency = desa2( inputData, &variance );
printf("Desa2: Mean freq = %f, var = %f, std dev = %f\n",
frequency, variance, sqrt(variance));

sampleCount += BLOCK;
numWords = fread(intData, sizeof(int16_t), BLOCK, inFile );
}

printf("\nFinished. sampleCount = %d\n",sampleCount);

fclose( inFile );
return 0;
}



For noise-free tones, the frequency estimation works well between about 200 to 3700 Hz for DESA-1 and about 200 to 1950 Hz for DESA-2. Note that the DESA-2 algorithm has the advantage of less computations per sample than DESA-1, but it won’t work beyond one quarter of the sampling frequency. If you are concerned with minimizing execution time, it is probably better to use DESA-1, because the input signal can be down-sampled by a factor of two more, which should likely make up for the fact that there are a few more operations for finding the frequency estimate.

A drawback of the DESA algorithms is that they perform poorly in noise. With additive white noise at 20 dB signal to noise ratio, and a 1000 Hz tone, the program above shows a large standard deviation of about 60 to 200 Hz.

While the above program is a neat demonstration of how to get the frequency of a tone without a lot of code, it is not practical for real world situations in which there is noise. In the next lesson I will show how to make a useful tone detector using the DESA algorithms.

If you are interested in finding the instantaneous frequency and amplitude of a general AM and FM modulated signal, I recommend taking a look at the thesis paper here (http://folk.uio.no/eivindkv/ek-thesis-2003-05-12-final-2.pdf ). It has a lot of analysis of the TKEO, and shows how to use a bank of bandpass filters to deal with noise. The following IEEE papers also explain the use of the energy operator:

Detection of Multi-Tone Signals Based on Energy Operators, Edgar F. Velez, Proceedings of the IEEE-SP International Symposium on Time-Frequency and Time-Scale Analysis 1994

Energy Separation in Signal Modulations with Applications to Speech Analysis, Petros Maragos, James F. Kaiser, Thomas F. Quatieri, IEEE Transactions on Signal Processing,Vol. 41, No. 10, October 1993