Improvement of Goertzel filter. Using Hamming window now. Add test routine.

.gitignore

@@ -84,6 +84,7 @@ src/test/test_filter
 src/test/test_sendevolumenregler
 src/test/test_compandor
 src/test/test_emphasis
+src/test/test_goertzel
 src/test/test_dtmf
 src/test/test_dms
 src/test/test_sms

src/amps/dsp.c

@@ -114,7 +114,7 @@
 #define TACS_SAT_DEVIATION	(1700.0 / TACS_SPEECH_DEVIATION)	/* no emphasis (panasonic / TI) */
 #define TACS_MAX_DISPLAY	(8000.0 / TACS_SPEECH_DEVIATION)	/* no emphasis */
 #define TACS_BITRATE		8000
-#define SAT_DURATION		0.05	/* duration of SAT signal measurement */
+#define SAT_BANDWIDTH		30.0	/* distance between two SAT tones, also bandwidth for goertzel filter */
 #define SAT_QUALITY		0.85	/* quality needed to detect SAT signal */
 #define SAT_PRINT		10	/* print sat measurement every 0.5 seconds */
 #define DTX_LEVEL		0.50	/* SAT level needed to mute/unmute */
@@ -133,7 +133,7 @@ static double sat_freq[4] = {
 	5970.0,
 	6000.0,
 	6030.0,
-	5800.0, /* noise level to check against */
+	5790.0, /* noise level to check against */
 };
 
 static sample_t dsp_sine_sat[65536];
@@ -243,8 +243,11 @@ int dsp_init_sender(amps_t *amps, int tolerant)
 	amps->fsk_deviation = (!tacs) ? AMPS_FSK_DEVIATION : TACS_FSK_DEVIATION;
 	dsp_init_ramp(amps);
 
-	/* allocate ring buffer for SAT signal detection */
+	/* allocate ring buffer for SAT signal detection
-	amps->sat_samples = (int)((double)amps->sender.samplerate * SAT_DURATION + 0.5);
+	 * the bandwidth of the Goertzel filter is the reciprocal of the duration
+	 * we half our bandwidth, so that other supervisory signals will be canceled out completely by goertzel filter
+	 */
+	amps->sat_samples = (int)((double)amps->sender.samplerate * (1.0 / (SAT_BANDWIDTH / 2.0)) + 0.5);
 	spl = calloc(sizeof(*spl), amps->sat_samples);
 	if (!spl) {
 		PDEBUG(DDSP, DEBUG_ERROR, "No memory!\n");
@@ -732,7 +735,7 @@ static void sender_receive_frame(amps_t *amps, sample_t *samples, int length)
 
 
 /* decode SAT and signaling tone */
-/* compare supervisory signal against noise floor on 5800 Hz */
+/* compare supervisory signal against noise floor at 5790 Hz */
 static void sat_decode(amps_t *amps, sample_t *samples, int length)
 {
 	double result[3], sat_quality, sig_quality, sat_level, sig_level;
@@ -742,8 +745,8 @@ static void sat_decode(amps_t *amps, sample_t *samples, int length)
 	audio_goertzel(&amps->sat_goertzel[4], samples, length, 0, &result[2], 1);
 
 	/* normalize sat level and signaling tone level */
-	sat_level = result[0] / ((!tacs) ? AMPS_SAT_DEVIATION : TACS_SAT_DEVIATION) / 0.63662;
+	sat_level = result[0] / ((!tacs) ? AMPS_SAT_DEVIATION : TACS_SAT_DEVIATION);
-	sig_level = result[2] / ((!tacs) ? AMPS_FSK_DEVIATION : TACS_FSK_DEVIATION) / 0.63662;
+	sig_level = result[2] / ((!tacs) ? AMPS_FSK_DEVIATION : TACS_FSK_DEVIATION);
 
 	/* get normalized quality of SAT and signaling tone */
 	sat_quality = (result[0] - result[1]) / result[0];

src/anetz/dsp.c

@@ -41,7 +41,7 @@
 #define TX_PEAK_TONE		(10500.0 / SPEECH_DEVIATION)	/* 10.5 kHz, no emphasis */
 #define TX_PEAK_PAGE		(15000.0 / SPEECH_DEVIATION)	/* 15 kHz, no emphasis */
 #define MAX_DISPLAY		(15000.0 / SPEECH_DEVIATION)	/* 15 kHz, no emphasis */
-#define CHUNK_DURATION		0.010	/* 10 ms */
+#define CHUNK_DURATION		0.010	/* 10 m = 100 Hz bandwidth (-7.6 DB @ +-100 Hz) */
 #define TONE_THRESHOLD		0.05
 #define QUAL_THRESHOLD		0.5
 
@@ -158,15 +158,15 @@ static void fsk_decode_chunk(anetz_t *anetz, sample_t *spl, int max)
 {
 	double level, result[2], quality[2];
 
-	level = audio_level(spl, max);
+	level = audio_mean_level(spl, max);
+	/* convert mean (if level comes from a sine curve) to peak value */
+	level = level * M_PI / 2.0 / TX_PEAK_TONE;
 
 	audio_goertzel(anetz->fsk_tone_goertzel, spl, max, 0, result, 2);
 
-	/* normalize quality of tones and level */
+	/* calculate quality of tones */
 	quality[0] = result[0] / level;
 	quality[1] = result[1] / level;
-	/* adjust level, so we get peak of sine curve */
-	level = level / 0.63662 / TX_PEAK_TONE;
 	/* show tones */
 	display_measurements_update(anetz->dmp_tone_level, level * 100.0, 0.0);
 	display_measurements_update(anetz->dmp_tone_quality, quality[1] * 100.0, 0.0);

src/libgoertzel/goertzel.c

@@ -35,7 +35,7 @@
  * For a perfect rectangualr wave, the result would equal the peak level.
  * For a sine wave the result would be factor (2 / PI) below peak level.
  */
-double audio_level(sample_t *samples, int length)
+double audio_mean_level(sample_t *samples, int length)
 {
 	double level, sk;
 	int n;
@@ -54,8 +54,23 @@ double audio_level(sample_t *samples, int length)
 	return level;
 }
 
+/* use hamming window */
+double window[256];
+int window_generated = 0;
+
+static void gen_window(void)
+{
+	int i;
+
+	for (i = 0; i < 256; i++)
+		window[i] = 0.54 - 0.46 * cos(2.0 * M_PI * (double)i / 256.0);
+	window_generated = 1;
+}
+
 void audio_goertzel_init(goertzel_t *goertzel, double freq, int samplerate)
 {
+	if (!window_generated)
+		gen_window();
 	memset(goertzel, 0, sizeof(*goertzel));
 	goertzel->coeff = 2.0 * cos(2.0 * M_PI * freq / (double)samplerate);
 }
@@ -67,10 +82,11 @@ void audio_goertzel_init(goertzel_t *goertzel, double freq, int samplerate)
 /* filter frequencies and return their levels
  *
  * samples: pointer to sample buffer
- * length: length of buffer
+ * length: length of buffer: -7.4 dB @ +-(1 / duration) Hz and -INF @ +-(1 / duration * 2) Hz
+ *  -> if duration is 10 ms, we got -7.4 dB @ +-100 Hz and -INF at +-200 Hz
  * offset: for ring buffer, start here and wrap around to 0 when length has been hit
  * coeff: array of coefficients (coeff << 15)
- * result: array of result levels (average value of the sine, that is 1 / (PI/2) of the sine's peak)
+ * result: array of result levels (peak value of the target frequency)
  * k: number of frequencies to check
  */
 void audio_goertzel(goertzel_t *goertzel, sample_t *samples, int length, int offset, double *result, int k)
@@ -87,7 +103,7 @@ void audio_goertzel(goertzel_t *goertzel, sample_t *samples, int length, int off
 		cos2pik = goertzel[i].coeff;
 		/* note: after 'length' cycles, offset is restored to its initial value */
 		for (n = 0; n < length; n++) {
-			sk = (cos2pik * sk1) - sk2 + samples[offset++];
+			sk = (cos2pik * sk1) - sk2 + samples[offset++] * window[n * 256 / length];
 			sk2 = sk1;
 			sk1 = sk;
 			if (offset == length)
@@ -98,7 +114,7 @@ void audio_goertzel(goertzel_t *goertzel, sample_t *samples, int length, int off
 			(sk * sk) -
 			(cos2pik * sk * sk2) +
 			(sk2 * sk2)
-				) / (double)length * 2.0 * 0.63662; /* 1 / (PI/2) */
+				) / (double)length * 4 / 1.08;
 	}
 }
 

src/libgoertzel/goertzel.h

@@ -1,5 +1,5 @@
 
-double audio_level(sample_t *samples, int length);
+double audio_mean_level(sample_t *samples, int length);
 
 typedef struct goertzel {
 	double coeff;

src/nmt/dsp.c

@@ -64,7 +64,7 @@
 #define MAX_DISPLAY		1.4	/* something above speech level */
 #define DIALTONE_HZ		425.0	/* dial tone frequency */
 #define TX_PEAK_DIALTONE	1.0	/* dial tone peak FIXME: Not found in the specs! */
-#define SUPER_DURATION		0.25	/* duration of supervisory signal measurement */
+#define SUPER_BANDWIDTH		30.0	/* distance between two SAT tones, also bandwidth for goertzel filter */
 #define SUPER_PRINT		2	/* print supervisory signal measurement every 0.5 seconds */
 #define SUPER_LOST_COUNT	4	/* number of measures to loose supervisory signal */
 #define SUPER_DETECT_COUNT	6	/* number of measures to detect supervisory signal */
@@ -76,7 +76,7 @@ static double super_freq[5] = {
 	3985.0, /* 0-Signal 2 */
 	4015.0, /* 0-Signal 3 */
 	4045.0, /* 0-Signal 4 */
-	3900.0, /* noise level to check against */
+	3895.0, /* noise level to check against */
 };
 
 /* table for fast sine generation */
@@ -129,8 +129,11 @@ int dsp_init_sender(nmt_t *nmt, double deviation_factor)
 		return -EINVAL;
 	}
 
-	/* allocate ring buffer for supervisory signal detection */
+	/* allocate ring buffer for SAT signal detection
-	nmt->super_samples = (int)((double)nmt->sender.samplerate * SUPER_DURATION + 0.5);
+	 * the bandwidth of the Goertzel filter is the reciprocal of the duration
+	 * we half our bandwidth, so that other supervisory signals will be canceled out completely by goertzel filter
+	 */
+	nmt->super_samples = (int)((double)nmt->sender.samplerate * (1.0 / (SUPER_BANDWIDTH / 2)) + 0.5);
 	spl = calloc(1, nmt->super_samples * sizeof(*spl));
 	if (!spl) {
 		PDEBUG(DDSP, DEBUG_ERROR, "No memory!\n");
@@ -260,7 +263,7 @@ static void fsk_receive_bit(void *inst, int bit, double quality, double level)
 	nmt_receive_frame(nmt, nmt->rx_frame, quality, level, frames_elapsed);
 }
 
-/* compare supervisory signal against noise floor on 3900 Hz */
+/* compare supervisory signal against noise floor around 3895 Hz */
 static void super_decode(nmt_t *nmt, sample_t *samples, int length)
 {
 	double result[2], level, quality;
@@ -269,7 +272,7 @@ static void super_decode(nmt_t *nmt, sample_t *samples, int length)
 	audio_goertzel(&nmt->super_goertzel[4], samples, length, 0, &result[1], 1); /* noise floor detection */
 
 	/* normalize supervisory level */
-	level = result[0] / 0.63662 / TX_PEAK_SUPER;
+	level = result[0] / TX_PEAK_SUPER;
 
 	quality = (result[0] - result[1]) / result[0];
 	if (quality < 0)

src/test/Makefile.am

@@ -5,6 +5,7 @@ noinst_PROGRAMS = \
 	test_sendevolumenregler \
 	test_compandor \
 	test_emphasis \
+	test_goertzel \
 	test_dtmf \
 	test_dms \
 	test_sms \
@@ -150,3 +151,12 @@ test_v27scrambler_LDADD = \
 	$(top_builddir)/src/libv27/libv27.a \
 	-lm
 
+test_goertzel_SOURCES = test_goertzel.c dummy.c
+
+test_goertzel_LDADD = \
+	$(COMMON_LA) \
+	$(top_builddir)/src/libdebug/libdebug.a \
+	$(top_builddir)/src/libgoertzel/libgoertzel.a \
+	$(top_builddir)/src/liboptions/liboptions.a \
+	-lm
+

src/test/test_goertzel.c

@@ -0,0 +1,51 @@
+#include <stdio.h>
+#include <stdint.h>
+#include <math.h>
+#include <string.h>
+#include "../libsample/sample.h"
+#include "../libgoertzel/goertzel.h"
+#include "../libdebug/debug.h"
+
+#define level2db(level)		(20 * log10(level))
+#define db2level(db)		pow(10, (double)db / 20.0)
+
+#define SAMPLERATE	48000
+
+static void gen_samples(sample_t *samples, double freq)
+{
+	int i;
+	double value;
+
+	for (i = 0; i < SAMPLERATE; i++) {
+		value = cos(2.0 * M_PI * freq / (double)SAMPLERATE * (double)i);
+		samples[i] = value;
+	}
+}
+
+int num_kanal;
+
+int main(void)
+{
+	goertzel_t goertzel;
+	sample_t samples[SAMPLERATE];
+	double frequency = 1000;
+	double duration = 1.0/100.0;
+	double level;
+	double i;
+
+	printf("testing goertzel with frequency %.1f and duration 1 / %.0f\n", frequency, 1.0 / duration);
+
+	for (i = 700; i < 1301; i = i + 10) {
+		gen_samples(samples, (double)i);
+		audio_goertzel_init(&goertzel, frequency, SAMPLERATE);
+		audio_goertzel(&goertzel, samples, SAMPLERATE * duration, 0, &level, 1);
+		printf("%s%.0f Hz: %.1f dB", debug_db(level), i, level2db(level));
+		if ((int)round(i) == (int)round(frequency))
+			printf(" level=%.6f\n", level);
+		else
+			printf("\n");
+	}
+
+	return 0;
+}
+