From dde4113e611dd0764e42406057589601d597da7f Mon Sep 17 00:00:00 2001 From: Andreas Eversberg Date: Sun, 28 Nov 2021 13:11:19 +0100 Subject: [PATCH] Improvement of Goertzel filter. Using Hamming window now. Add test routine. --- .gitignore | 1 + src/amps/dsp.c | 17 +++++++------ src/anetz/dsp.c | 10 ++++---- src/libgoertzel/goertzel.c | 26 +++++++++++++++---- src/libgoertzel/goertzel.h | 2 +- src/nmt/dsp.c | 15 ++++++----- src/test/Makefile.am | 10 ++++++++ src/test/test_goertzel.c | 51 ++++++++++++++++++++++++++++++++++++++ 8 files changed, 108 insertions(+), 24 deletions(-) create mode 100644 src/test/test_goertzel.c diff --git a/.gitignore b/.gitignore index 14e0170..383bec8 100644 --- a/.gitignore +++ b/.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 diff --git a/src/amps/dsp.c b/src/amps/dsp.c index 21af267..efd06c0 100644 --- a/src/amps/dsp.c +++ b/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 */ - amps->sat_samples = (int)((double)amps->sender.samplerate * SAT_DURATION + 0.5); + /* allocate ring buffer for SAT signal detection + * 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(&s->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; - sig_level = result[2] / ((!tacs) ? AMPS_FSK_DEVIATION : TACS_FSK_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); /* get normalized quality of SAT and signaling tone */ sat_quality = (result[0] - result[1]) / result[0]; diff --git a/src/anetz/dsp.c b/src/anetz/dsp.c index 40144a5..1bc20a1 100644 --- a/src/anetz/dsp.c +++ b/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); diff --git a/src/libgoertzel/goertzel.c b/src/libgoertzel/goertzel.c index cb06f99..d928175 100644 --- a/src/libgoertzel/goertzel.c +++ b/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; } } diff --git a/src/libgoertzel/goertzel.h b/src/libgoertzel/goertzel.h index 58fcbb8..0a8e0b5 100644 --- a/src/libgoertzel/goertzel.h +++ b/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; diff --git a/src/nmt/dsp.c b/src/nmt/dsp.c index d942c03..627c120 100644 --- a/src/nmt/dsp.c +++ b/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 */ - nmt->super_samples = (int)((double)nmt->sender.samplerate * SUPER_DURATION + 0.5); + /* allocate ring buffer for SAT signal detection + * 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) diff --git a/src/test/Makefile.am b/src/test/Makefile.am index e53579c..4ae1e3e 100644 --- a/src/test/Makefile.am +++ b/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 + diff --git a/src/test/test_goertzel.c b/src/test/test_goertzel.c new file mode 100644 index 0000000..ee1ba0f --- /dev/null +++ b/src/test/test_goertzel.c @@ -0,0 +1,51 @@ +#include +#include +#include +#include +#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; +} +