Fixing compander

6 files changed, 231 insertions, 57 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 7721c5b..c636965 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,3 +1,3 @@
 AUTOMAKE_OPTIONS = foreign
-SUBDIRS = common anetz bnetz nmt
+SUBDIRS = common anetz bnetz nmt test
 
diff --git a/src/common/compander.c b/src/common/compander.c
index d42d6de..4413bb5 100644
--- a/src/common/compander.c
+++ b/src/common/compander.c
@@ -17,42 +17,47 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <stdio.h>
 #include <stdint.h>
 #include <string.h>
 #include <math.h>
 #include "compander.h"
 
-/* this is the 0 DB level that stays 0 DB after compression / espansion */
-#define ZERO_DB_LEVEL	16384.0
+//#define db2level(db)			pow(10, (double)db / 20.0)
 
-/* what factor shall the gain raise and fall after given attack/recovery time */
-#define ATTACK_FACTOR	1.5
-#define RECOVERY_FACTOR	0.75
+/* factor is the gain (raise and fall) after given attack/recovery time */
+#define COMPRESS_ATTACK_FACTOR		1.83		/* about 1.5 after 12 dB step up */
+#define COMPRESS_RECOVERY_FACTOR	0.44		/* about 0.75 after 12 dB step down */
+#define EXPAND_ATTACK_FACTOR		1.145		/* about 0.57 after 6 dB step up */
+#define EXPAND_RECOVERY_FACTOR		0.753		/* about 1.51 after 6 dB step down */
 
 /* Minimum level value to keep state */
-#define ENVELOP_MIN	0.001
+#define ENVELOPE_MIN	0.001
 
 static double sqrt_tab[10000];
 
 /*
- * Generate companding tables according to NMT specification
+ * Init compander according to ITU-T G.162 specification
  *
  * Hopefully this is correct
  *
  */
-void init_compander(compander_t *state, int samplerate, double attack_ms, double recovery_ms)
+void init_compander(compander_t *state, int samplerate, double attack_ms, double recovery_ms, int unaffected_level)
 {
 	int i;
 
 	memset(state, 0, sizeof(*state));
 
-	state->envelop_e = 1.0;
-	state->envelop_c = 1.0;
-	/* ITU-T G.162: 1.5 times the steady state after attack_ms */
-	state->step_up = pow(ATTACK_FACTOR, 1000.0 / attack_ms / (double)samplerate);
-
-	/* ITU-T G.162: 0.75 times the steady state after recovery_ms */
-	state->step_down = pow(RECOVERY_FACTOR, 1000.0 / recovery_ms / (double)samplerate);
+	state->c.peak = 1.0;
+	state->c.envelope = 1.0;
+	state->e.peak = 1.0;
+	state->e.envelope = 1.0;
+	state->c.step_up = pow(COMPRESS_ATTACK_FACTOR, 1000.0 / attack_ms / (double)samplerate);
+	state->c.step_down = pow(COMPRESS_RECOVERY_FACTOR, 1000.0 / recovery_ms / (double)samplerate);
+	state->e.step_up = pow(EXPAND_ATTACK_FACTOR, 1000.0 / attack_ms / (double)samplerate);
+	state->e.step_down = pow(EXPAND_RECOVERY_FACTOR, 1000.0 / recovery_ms / (double)samplerate);
+	state->c.unaffected = unaffected_level;
+	state->e.unaffected = unaffected_level;
 
 	// FIXME: make global, not at instance
 	for (i = 0; i < 10000; i++)
@@ -62,64 +67,85 @@ void init_compander(compander_t *state, int samplerate, double attack_ms, double
 void compress_audio(compander_t *state, int16_t *samples, int num)
 {
 	int32_t sample;
-	double value, envelop, step_up, step_down;
+	double value, peak, envelope, step_up, step_down, unaffected;
 	int i;
 
-	step_up = state->step_up;
-	step_down = state->step_down;
-	envelop = state->envelop_c;
+	step_up = state->c.step_up;
+	step_down = state->c.step_down;
+	peak = state->c.peak;
+	envelope = state->c.envelope;
+	unaffected = state->c.unaffected;
 
-//	printf("envelop=%.4f\n", envelop);
+//	printf("envelope=%.4f\n", envelope);
 	for (i = 0; i < num; i++) {
-		/* normalize sample value to 0 DB level */
-		value = (double)(*samples) / ZERO_DB_LEVEL;
+		/* normalize sample value to unaffected level */
+		value = (double)(*samples) / unaffected;
 
-		if (fabs(value) > envelop)
-			envelop *= step_up;
+		/* 'peak' is the level that raises directly with the signal
+		 * level, but falls with specified recovery rate. */
+		if (fabs(value) > peak)
+			peak = fabs(value);
 		else
-			envelop *= step_down;
-		if (envelop < ENVELOP_MIN)
-			envelop = ENVELOP_MIN;
+			peak *= step_down;
+		/* 'evelope' is the level that raises with the specified attack
+		 * rate to 'peak', but falls with specified recovery rate. */
+		if (peak > envelope)
+			envelope *= step_up;
+		else
+			envelope = peak;
+		if (envelope < ENVELOPE_MIN)
+			envelope = ENVELOPE_MIN;
 
-		value = value / sqrt_tab[(int)(envelop / 0.001)];
+		value = value / sqrt_tab[(int)(envelope / 0.001)];
+//if (i > 47000.0 && i < 48144)
+//printf("time=%.4f envelope=%.4fdb, value=%.4f\n", (double)i/48000.0, 20*log10(envelope), value);
 
 		/* convert back from 0 DB level to sample value */
-		sample = (int)(value * ZERO_DB_LEVEL);
+		sample = (int)(value * unaffected);
 		if (sample > 32767)
 			sample = 32767;
 		else if (sample < -32768)
 			sample = -32768;
 		*samples++ = sample;
 	}
+//exit(0);
 
-	state->envelop_c = envelop;
+	state->c.envelope = envelope;
+	state->c.peak = peak;
 }
 
 void expand_audio(compander_t *state, int16_t *samples, int num)
 {
 	int32_t sample;
-	double value, envelop, step_up, step_down;
+	double value, peak, envelope, step_up, step_down, unaffected;
 	int i;
 
-	step_up = state->step_up;
-	step_down = state->step_down;
-	envelop = state->envelop_e;
+	step_up = state->e.step_up;
+	step_down = state->e.step_down;
+	peak = state->e.peak;
+	envelope = state->e.envelope;
+	unaffected = state->e.unaffected;
 
 	for (i = 0; i < num; i++) {
 		/* normalize sample value to 0 DB level */
-		value = (double)(*samples) / ZERO_DB_LEVEL;
+		value = (double)(*samples) / unaffected;
 
-		if (fabs(value) > envelop)
-			envelop *= step_up;
+		/* for comments: see compress_audio() */
+		if (fabs(value) > peak)
+			peak = fabs(value);
+		else
+			peak *= step_down;
+		if (peak > envelope)
+			envelope *= step_up;
 		else
-			envelop *= step_down;
-		if (envelop < ENVELOP_MIN)
-			envelop = ENVELOP_MIN;
+			envelope = peak;
+		if (envelope < ENVELOPE_MIN)
+			envelope = ENVELOPE_MIN;
 
-		value = value * envelop;
+		value = value * envelope;
 
 		/* convert back from 0 DB level to sample value */
-		sample = (int)(value * ZERO_DB_LEVEL);
+		sample = (int)(value * unaffected);
 		if (sample > 32767)
 			sample = 32767;
 		else if (sample < -32768)
@@ -127,6 +153,7 @@ void expand_audio(compander_t *state, int16_t *samples, int num)
 		*samples++ = sample;
 	}
 
-	state->envelop_e = envelop;
+	state->e.envelope = envelope;
+	state->e.peak = peak;
 }
 
diff --git a/src/common/compander.h b/src/common/compander.h
index d48e284..53d1967 100644
--- a/src/common/compander.h
+++ b/src/common/compander.h
@@ -1,11 +1,21 @@
 typedef struct compander {
-	double	step_up;
-	double	step_down;
-	double	envelop_c;
-	double	envelop_e;
+	struct {
+		double	unaffected;
+		double	step_up;
+		double	step_down;
+		double	peak;
+		double	envelope;
+	} c;
+	struct {
+		double	unaffected;
+		double	step_up;
+		double	step_down;
+		double	peak;
+		double	envelope;
+	} e;
 } compander_t;
 
-void init_compander(compander_t *state, int samplerate, double attack_ms, double recovery_ms);
+void init_compander(compander_t *state, int samplerate, double attack_ms, double recovery_ms, int unaffected_level);
 void compress_audio(compander_t *state, int16_t *samples, int num);
 void expand_audio(compander_t *state, int16_t *samples, int num);
 
diff --git a/src/nmt/dsp.c b/src/nmt/dsp.c
index 79403b9..e4e1bd7 100644
--- a/src/nmt/dsp.c
+++ b/src/nmt/dsp.c
@@ -33,8 +33,9 @@
 #define PI			M_PI
 
 /* signalling */
-#define TX_PEAK_FSK		16384	/* peak amplitude of signalling FSK */
-#define TX_PEAK_SUPER		1638	/* peak amplitude of supervisory signal */
+#define TX_AUDIO_0dBm0		16384	/* works quite well */
+#define TX_PEAK_FSK		16384.0	/* peak amplitude of signalling FSK */
+#define TX_PEAK_SUPER		1638.0	/* peak amplitude of supervisory signal */
 #define BIT_RATE		1200	/* baud rate */
 #define STEPS_PER_BIT		10	/* step every 1/12000 sec */
 #define DIALTONE_HZ		425.0	/* dial tone frequency */
@@ -84,7 +85,7 @@ int dsp_init_sender(nmt_t *nmt)
 	int i;
 
 	/* attack (3ms) and recovery time (13.5ms) according to NMT specs */
-	init_compander(&nmt->cstate, 8000, 3.0, 13.5);
+	init_compander(&nmt->cstate, 8000, 3.0, 13.5, TX_AUDIO_0dBm0);
 
 	if ((nmt->sender.samplerate % (BIT_RATE * STEPS_PER_BIT))) {
 		PDEBUG(DDSP, DEBUG_ERROR, "Sample rate must be a multiple of %d bits per second.\n", BIT_RATE * STEPS_PER_BIT);
@@ -238,7 +239,8 @@ static void fsk_receive_bit(nmt_t *nmt, int bit, double quality, double level)
 
 	/* send telegramm */
 	frames_elapsed = (double)(nmt->rx_sample_count_current - nmt->rx_sample_count_last) / (double)(nmt->samples_per_bit * 166);
-	nmt_receive_frame(nmt, nmt->fsk_filter_frame, nmt->fsk_filter_qualitysum / 140.0, nmt->fsk_filter_levelsum / 140.0, frames_elapsed);
+	/* convert level so that received level at TX_PEAK_FSK results in 1.0 (100%) */
+	nmt_receive_frame(nmt, nmt->fsk_filter_frame, nmt->fsk_filter_qualitysum / 140.0, nmt->fsk_filter_levelsum / 140.0 * 32768.0 / TX_PEAK_FSK, frames_elapsed);
 }
 
 char *show_level(int value)
@@ -339,8 +341,8 @@ static void super_decode(nmt_t *nmt, int16_t *samples, int length)
 
 #if 0
 	/* normalize levels */
-	result[0] *= 32768.0 / (double)TX_PEAK_SUPER / 0.63662;
-	result[1] *= 32768.0 / (double)TX_PEAK_SUPER / 0.63662;
+	result[0] *= 32768.0 / TX_PEAK_SUPER / 0.63662;
+	result[1] *= 32768.0 / TX_PEAK_SUPER / 0.63662;
 	printf("signal=%.4f noise=%.4f\n", result[0], result[1]);
 #endif
 
@@ -412,7 +414,7 @@ void sender_receive(sender_t *sender, int16_t *samples, int length)
 	spl = nmt->fsk_filter_spl;
 	for (i = 0; i < length; i++) {
 #ifdef DEBUG_MODULATOR
-		printf("|%s|\n", show_level((int)((samples[i] / (double)TX_PEAK_FSK) * 50)+50));
+		printf("|%s|\n", show_level((int)((samples[i] / TX_PEAK_FSK) * 50)+50));
 #endif
 		spl[pos++] = samples[i];
 		if (nmt->fsk_filter_mute) {
@@ -441,7 +443,8 @@ void sender_receive(sender_t *sender, int16_t *samples, int length)
 		spl = nmt->sender.rxbuf;
 		pos = nmt->sender.rxbuf_pos;
 		for (i = 0; i < count; i++) {
-			spl[pos++] = down[i];
+#warning hacking: remove after preemphasis implementation
+			spl[pos++] = down[i] / 2;
 			if (pos == 160) {
 				call_tx_audio(nmt->sender.callref, spl, 160);
 				pos = 0;
diff --git a/src/test/Makefile.am b/src/test/Makefile.am
new file mode 100644
index 0000000..2375047
--- /dev/null
+++ b/src/test/Makefile.am
@@ -0,0 +1,12 @@
+AM_CPPFLAGS = -Wall -g $(all_includes)
+
+noinst_PROGRAMS = \
+	test_compander
+
+test_compander_SOURCES = test_compander.c
+
+test_compander_LDADD = \
+	$(COMMON_LA) \
+	$(top_builddir)/src/common/libcommon.a \
+	-lm
+
diff --git a/src/test/test_compander.c b/src/test/test_compander.c
new file mode 100644
index 0000000..d1661b4
--- /dev/null
+++ b/src/test/test_compander.c
@@ -0,0 +1,122 @@
+#include <stdio.h>
+#include <stdint.h>
+#include <math.h>
+#include <string.h>
+#include "../common/compander.h"
+
+#define level2db(level)		(20 * log10(level))
+#define db2level(db)		pow(10, (double)db / 20.0)
+
+#define SAMPLERATE		48000
+#define ATTACK_MS		15.0
+#define RECOVERY_MS		15.0
+#define UNAFFECTED		10000.0
+
+static double test_frequency[3] = { 2000.0, 4000.0, 1000.0 };
+
+static int16_t samples_4db[SAMPLERATE];
+static int16_t samples_16db[SAMPLERATE];
+static int16_t samples_2db[SAMPLERATE];
+static int16_t samples_8db[SAMPLERATE];
+static int16_t samples_0db[SAMPLERATE];
+
+/* generate 2 samples: one with -4 dB, the other with -16 dB */
+static void generate_test_sample(double test_frequency)
+{
+	int i;
+	double value;
+
+	for (i = 0; i < SAMPLERATE; i++) {
+		value = cos(2.0 * M_PI * test_frequency / (double)SAMPLERATE * i);
+		samples_4db[i] = (int)(UNAFFECTED * value * db2level(-4));
+		samples_16db[i] = (int)(UNAFFECTED * value * db2level(-16));
+		samples_2db[i] = (int)(UNAFFECTED * value * db2level(-2));
+		samples_8db[i] = (int)(UNAFFECTED * value * db2level(-8));
+		samples_0db[i] = (int)(UNAFFECTED * value);
+	}
+}
+
+static void check_level(int16_t *samples, double duration, const char *desc, double target_db)
+{
+	int i;
+	int last = 0, envelop = 0;
+	int up = 0;
+	double factor;
+
+	int when = (int)((double)SAMPLERATE + (double)SAMPLERATE * duration / 1000.0 + 0.5);
+
+	for (i = 0; i < when; i++) {
+		if (last < samples[i]) {
+			up = 1;
+		} else if (last > samples[i]) {
+			if (up) {
+				envelop = last;
+			}
+			up = 0;
+		}
+#if 0
+		if (i >= (SAMPLERATE-(SAMPLERATE/100)) && (i % (SAMPLERATE/(SAMPLERATE/10))) == 0)
+			printf("%s: envelop = %.4f (when=%d)\n", desc, level2db((double)envelop / UNAFFECTED), i);
+#endif
+		last = samples[i];
+	}
+	factor = (envelop / UNAFFECTED) / db2level(target_db);
+	printf("%s: envelop after the instance of %.1f ms is %.4f db factor =%.4f\n", desc, duration, level2db((double)envelop / UNAFFECTED), factor);
+}
+
+int main(void)
+{
+	compander_t cstate;
+	int16_t samples[SAMPLERATE * 2];
+	int f;
+
+	init_compander(&cstate, SAMPLERATE, ATTACK_MS, RECOVERY_MS, UNAFFECTED);
+
+	for (f = 0; f < 3; f++) {
+		/* -16 and -4 dB */
+		generate_test_sample(test_frequency[f]);
+
+		/* low to high transition */
+		memcpy(samples, samples_16db, SAMPLERATE * 2);
+		memcpy(samples + SAMPLERATE, samples_4db, SAMPLERATE * 2);
+
+		compress_audio(&cstate, samples, SAMPLERATE * 2);
+
+		check_level(samples, ATTACK_MS, "compressor attack", -2.0);
+
+		/* high to low transition */
+		memcpy(samples, samples_4db, SAMPLERATE * 2);
+		memcpy(samples + SAMPLERATE, samples_16db, SAMPLERATE * 2);
+
+		compress_audio(&cstate, samples, SAMPLERATE * 2);
+
+		check_level(samples, RECOVERY_MS, "compressor recovery", -8.0);
+
+		/* low to high transition */
+		memcpy(samples, samples_8db, SAMPLERATE * 2);
+		memcpy(samples + SAMPLERATE, samples_2db, SAMPLERATE * 2);
+
+		expand_audio(&cstate, samples, SAMPLERATE * 2);
+
+		check_level(samples, ATTACK_MS, "expander attack", -4.0);
+
+		/* high to low transition */
+		memcpy(samples, samples_2db, SAMPLERATE * 2);
+		memcpy(samples + SAMPLERATE, samples_8db, SAMPLERATE * 2);
+
+		expand_audio(&cstate, samples, SAMPLERATE * 2);
+
+		check_level(samples, RECOVERY_MS, "expander recovery", -16.0);
+
+		/* 0 DB */
+		memcpy(samples, samples_0db, SAMPLERATE * 2);
+		memcpy(samples + SAMPLERATE, samples_0db, SAMPLERATE * 2);
+
+		check_level(samples, RECOVERY_MS, "unaffected level", 0.0);
+
+		puts("");
+	}
+
+	return 0;
+}
+