redesigned ddc, replaced delay lines with common buffer

Makefile

@@ -28,7 +28,7 @@ clean:
 
 encode: encode.o mmap_file.o pcm.o wav.o alsa.o yuv.o img.o ppm.o sdl.o
 
-decode: decode.o mmap_file.o pcm.o wav.o alsa.o window.o ddc.o delay.o yuv.o img.o ppm.o sdl.o
+decode: decode.o mmap_file.o pcm.o wav.o alsa.o window.o ddc.o buffer.o yuv.o img.o ppm.o sdl.o
 
-debug: debug.o mmap_file.o pcm.o wav.o alsa.o window.o ddc.o delay.o yuv.o img.o ppm.o sdl.o
+debug: debug.o mmap_file.o pcm.o wav.o alsa.o window.o ddc.o buffer.o yuv.o img.o ppm.o sdl.o
 

buffer.c

@@ -5,31 +5,34 @@ To the extent possible under law, the author(s) have dedicated all copyright and
 You should have received a copy of the CC0 Public Domain Dedication along with this software. If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.
 */
 
-
+#include "buffer.h"
-#include "delay.h"
 #include <stdlib.h>
 
-float do_delay(struct delay *d, float input)
+float *do_buffer(struct buffer *d, float input)
 {
-	d->s[d->last] = input;
+	d->s[d->last0] = input;
-	d->last = (d->last + 1) < d->len ? d->last + 1 : 0;
+	d->s[d->last1] = input;
-	return d->s[d->last];
+	d->last0 = (d->last0 - 1) < 0 ? d->len : d->last0 - 1;
+	d->last1 = (d->last1 - 1) < 0 ? d->len : d->last1 - 1;
+	int last = d->last0 < d->last1 ? d->last0 : d->last1;
+	return d->s + last;
 }
 
-struct delay *alloc_delay(int samples)
+struct buffer *alloc_buffer(int samples)
 {
-	int len = samples + 1;
+	int len = 2 * samples;
-	struct delay *d = malloc(sizeof(struct delay));
+	struct buffer *d = malloc(sizeof(struct buffer));
 	d->s = malloc(sizeof(float) * len);
-	d->last = 0;
+	d->last0 = 0;
+	d->last1 = samples;
 	d->len = len;
 	for (int i = 0; i < len; i++)
 		d->s[i] = 0.0;
 	return d;
 }
-void free_delay(struct delay *delay)
+void free_buffer(struct buffer *buffer)
 {
-	free(delay->s);
+	free(buffer->s);
-	free(delay);
+	free(buffer);
 }
 

buffer.h

@@ -5,17 +5,17 @@ To the extent possible under law, the author(s) have dedicated all copyright and
 You should have received a copy of the CC0 Public Domain Dedication along with this software. If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.
 */
 
-
+#ifndef BUFFER_H
-#ifndef DELAY_H
+#define BUFFER_H
-#define DELAY_H
+struct buffer {
-struct delay {
 	float *s;
-	int last;
+	int last0;
+	int last1;
 	int len;
 };
 
-float do_delay(struct delay *, float);
+float *do_buffer(struct buffer *d, float input);
-struct delay *alloc_delay(int);
+struct buffer *alloc_buffer(int samples);
-void free_delay(struct delay *);
+void free_buffer(struct buffer *buffer);
 #endif
 

ddc.c

@@ -5,78 +5,70 @@ To the extent possible under law, the author(s) have dedicated all copyright and
 You should have received a copy of the CC0 Public Domain Dedication along with this software. If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.
 */
 
-
+#include <complex.h> 
 #include <math.h>
 #include <stdlib.h>
-#include <complex.h> 
 #include "window.h"
 #include "ddc.h"
 
-void do_ddc(struct ddc *ddc, float *input, float complex *output)
+void do_ddc(struct ddc *ddc, float *input, complex float *output)
 {
-	// this works only for L <= M
+	int N = ddc->N;
-	for (int k = 0, last = ddc->last, in = 0; k < ddc->L; k++) {
+	int M = ddc->M;
-		while (ddc->skip < ddc->M) {
+	int L = ddc->L;
-			ddc->s[ddc->last] = input[in++];
+	int offset = 0;
-			last = ddc->last;
+	for (int k = 0; k < L; k++) {
-			ddc->last = (ddc->last + 1) < ddc->samples ? ddc->last + 1 : 0;
+		float *x = input + (((L * M - 1) - offset) / L);
-			ddc->skip += ddc->L;
+		complex float *b = ddc->b + (N * (offset % L));
-		}
+		offset += M;
 
-		ddc->skip %= ddc->M;
+		complex float sum = 0.0;
-
+		for (int i = 0; i < N; i++)
-		float complex sum = 0.0;
+			sum += b[i] * x[i];
-		for (int i = ddc->offset; i < ddc->taps; i += ddc->L) {
-			sum += ddc->b[i] * ddc->s[last];
-			last += last ? - 1 : ddc->samples - 1;
-		}
-
-		ddc->offset = (ddc->offset + ddc->M) % ddc->L;
 
 		output[k] = ddc->osc * sum;
 		ddc->osc *= ddc->d;
-//		ddc->osc /= cabsf(ddc->osc); // not really needed
+		// ddc->osc /= cabsf(ddc->osc); // not really needed
 	}
 }
-struct ddc *alloc_ddc(float freq, float bw, float step, int taps, int L, int M, float (*window)(float, float, float), float a)
+struct ddc *alloc_ddc(int L, int M, float carrier, float bw, float rate, int taps, float (*window)(float, float, float), float a)
 {
-	float lstep = step / (float)L;
+	float lstep = 1.0 / ((float)L * rate);
-	float ostep = step * (float)M / (float)L;
+	float ostep = (float)M * lstep;
 	struct ddc *ddc = malloc(sizeof(struct ddc));
-	ddc->taps = taps;
+	ddc->N = (taps + L - 1) / L;
-	ddc->samples = (taps + L - 1) / L;
+	ddc->b = malloc(sizeof(complex float) * ddc->N * L);
-	ddc->b = malloc(sizeof(float complex) * ddc->taps);
-	ddc->s = malloc(sizeof(float) * ddc->samples);
 	ddc->osc = I;
-	ddc->d = cexpf(-I * 2.0 * M_PI * freq * ostep);
+	ddc->d = cexpf(-I * 2.0 * M_PI * carrier * ostep);
-	ddc->offset = (M - 1) % L;
-	ddc->last = 0;
-	ddc->skip = 0;
 	ddc->L = L;
 	ddc->M = M;
-	for (int i = 0; i < ddc->samples; i++)
+	complex float *b = malloc(sizeof(complex float) * taps);
-		ddc->s[i] = 0.0;
 	float sum = 0.0;
-	for (int i = 0; i < ddc->taps; i++) {
+	for (int i = 0; i < taps; i++) {
-		float N = (float)ddc->taps;
+		float N = taps;
-		float n = (float)i;
+		float n = i;
 		float x = n - (N - 1.0) / 2.0;
-		float l = 2.0 * M_PI * bw * lstep;
+		float l = 2.0 * bw * lstep;
-		float w = window(n, ddc->taps, a);
+		float h = l * sinc(l * x);
-		float h = 0.0 == x ? l / M_PI : sinf(l * x) / (x * M_PI);
+		float w = window(n, N, a);
-		float b = w * h;
+		float t = w * h;
-		sum += b;
+		sum += t;
-		complex float o = cexpf(I * 2.0 * M_PI * freq * lstep * n);
+		complex float o = cexpf(I * 2.0 * M_PI * carrier * lstep * n);
-		ddc->b[i] = b * o * (float)L;
+		b[i] = t * o;
 	}
-	for (int i = 0; i < ddc->taps; i++)
+	for (int i = 0; i < taps; i++)
-		ddc->b[i] /= sum;
+		b[i] /= sum;
+	for (int i = 0; i < ddc->N * L; i++)
+		ddc->b[i] = 0.0;
+	for (int i = 0; i < L; i++)
+		for (int j = i, k = 0; j < taps; j += L, k++)
+			ddc->b[i * ddc->N + k] = b[j];
+	free(b);
 	return ddc;
 }
 void free_ddc(struct ddc *ddc)
 {
 	free(ddc->b);
-	free(ddc->s);
 	free(ddc);
 }
 

ddc.h

@@ -5,28 +5,22 @@ To the extent possible under law, the author(s) have dedicated all copyright and
 You should have received a copy of the CC0 Public Domain Dedication along with this software. If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.
 */
 
-
 #ifndef DDC_H
 #define DDC_H
 #include "window.h"
 
 struct ddc {
-	float complex *b;
+	complex float *b;
-	float *s;
+	complex float osc;
-	float complex osc;
+	complex float d;
-	float complex d;
+	int N;
-	int offset;
-	int skip;
-	int last;
-	int taps;
-	int samples;
 	int L;
 	int M;
 };
 
-void do_ddc(struct ddc *, float *, float complex *);
+void do_ddc(struct ddc *ddc, float *input, complex float *output);
-struct ddc *alloc_ddc(float, float, float, int, int, int, float (*)(float, float, float), float);
+struct ddc *alloc_ddc(int L, int M, float carrier, float bw, float rate, int taps, float (*window)(float, float, float), float a);
-void free_ddc(struct ddc *);
+void free_ddc(struct ddc *ddc);
 
 #endif
 

debug.c

@@ -15,7 +15,7 @@ You should have received a copy of the CC0 Public Domain Dedication along with t
 #include "mmap_file.h"
 #include "pcm.h"
 #include "ddc.h"
-#include "delay.h"
+#include "buffer.h"
 #include "yuv.h"
 #include "utils.h"
 #include "img.h"
@@ -47,7 +47,6 @@ int main(int argc, char **argv)
 	if (channels > 1)
 		fprintf(stderr, "using first of %d channels\n", channels);
 
-	const float step = 1.0 / rate;
 	float complex cnt_last = -I;
 	float complex dat_last = -I;
 
@@ -105,18 +104,19 @@ int main(int argc, char **argv)
 	float drate = rate * (float)factor_L / (float)factor_M;
 	float dstep = 1.0 / drate;
 	fprintf(stderr, "using factor of %ld/%ld, working at %.2fhz\n", factor_L, factor_M, drate);
-	float *cnt_amp = malloc(sizeof(float) * factor_M);
-	float *dat_amp = malloc(sizeof(float) * factor_M);
 	float complex *cnt_q = malloc(sizeof(float complex) * factor_L);
 	float complex *dat_q = malloc(sizeof(float complex) * factor_L);
 	// same factor to keep life simple and have accurate horizontal sync
-	struct ddc *cnt_ddc = alloc_ddc(1200.0, 200.0, step, cnt_taps, factor_L, factor_M, kaiser, 2.0);
+	struct ddc *cnt_ddc = alloc_ddc(factor_L, factor_M, 1200.0, 200.0, rate, cnt_taps, kaiser, 2.0);
-	struct ddc *dat_ddc = alloc_ddc(1900.0, 800.0, step, dat_taps, factor_L, factor_M, kaiser, 2.0);
+	struct ddc *dat_ddc = alloc_ddc(factor_L, factor_M, 1900.0, 800.0, rate, dat_taps, kaiser, 2.0);
 	// delay input by phase shift of other filter to synchronize outputs
-	struct delay *cnt_delay = alloc_delay((dat_taps - 1) / (2 * factor_L));
+	int cnt_delay = (dat_taps - 1) / (2 * factor_L);
-	struct delay *dat_delay = alloc_delay((cnt_taps - 1) / (2 * factor_L));
+	int dat_delay = (cnt_taps - 1) / (2 * factor_L);
 
-	short *buff = (short *)malloc(sizeof(short) * channels * factor_M);
+	short *pcm_buff = (short *)malloc(sizeof(short) * channels * factor_M);
+
+	// 0.1 second history + enough room for delay and taps
+	struct buffer *buffer = alloc_buffer(0.1 * rate + 2 * fmaxf(cnt_delay, dat_delay));
 
 	const float sync_porch_len = 0.003;
 	const float porch_len = 0.0015; (void)porch_len;
@@ -143,15 +143,14 @@ int main(int argc, char **argv)
 	for (int out = factor_L;; out++, hor_ticks++, cal_ticks++, vis_ticks++) {
 		if (out >= factor_L) {
 			out = 0;
-			if (!read_pcm(pcm, buff, factor_M))
+			if (!read_pcm(pcm, pcm_buff, factor_M))
 				break;
-			for (int j = 0; j < factor_M; j++) {
+			float *buff = 0;
-				float amp = (float)buff[j * channels] / 32767.0;
+			for (int j = 0; j < factor_M; j++)
-				cnt_amp[j] = do_delay(cnt_delay, amp);
+				buff = do_buffer(buffer, (float)pcm_buff[j * channels] / 32767.0);
-				dat_amp[j] = do_delay(dat_delay, amp);
+
-			}
+			do_ddc(cnt_ddc, buff + cnt_delay, cnt_q);
-			do_ddc(cnt_ddc, cnt_amp, cnt_q);
+			do_ddc(dat_ddc, buff + dat_delay, dat_q);
-			do_ddc(dat_ddc, dat_amp, dat_q);
 		}
 
 		float cnt_freq = fclampf(1200.0 + cargf(cnt_q[out] * conjf(cnt_last)) / (2.0 * M_PI * dstep), 1100.0, 1300.0);
@@ -417,8 +416,8 @@ int main(int argc, char **argv)
 
 	free_ddc(cnt_ddc);
 	free_ddc(dat_ddc);
-	free(cnt_amp);
+	free_buffer(buffer);
-	free(dat_amp);
+	free(pcm_buff);
 
 	return 0;
 }

decode.c

@@ -14,7 +14,7 @@ You should have received a copy of the CC0 Public Domain Dedication along with t
 #include <time.h>
 #include "pcm.h"
 #include "ddc.h"
-#include "delay.h"
+#include "buffer.h"
 #include "yuv.h"
 #include "utils.h"
 #include "img.h"
@@ -289,23 +289,20 @@ int demodulate(struct pcm *pcm, float *cnt_freq, float *dat_freq, float *drate)
 	static float dstep;
 	static float complex cnt_last = -I;
 	static float complex dat_last = -I;
-	static float *cnt_amp;
-	static float *dat_amp;
 	static float complex *cnt_q;
 	static float complex *dat_q;
 	static struct ddc *cnt_ddc;
 	static struct ddc *dat_ddc;
-	static struct delay *cnt_delay;
+	static struct buffer *buffer;
-	static struct delay *dat_delay;
+	static int cnt_delay;
-	static short *buff;
+	static int dat_delay;
+	static short *pcm_buff;
 
 	static int init = 0;
 	if (!init) {
 		init = 1;
 		rate = rate_pcm(pcm);
 		channels = channels_pcm(pcm);
-		const float step = 1.0 / rate;
-
 #if DN && UP
 		// 320 / 0.088 = 160 / 0.044 = 40000 / 11 = 3636.(36)~ pixels per second for Y, U and V
 		factor_L = 40000;
@@ -331,18 +328,19 @@ int demodulate(struct pcm *pcm, float *cnt_freq, float *dat_freq, float *drate)
 		*drate = rate * (float)factor_L / (float)factor_M;
 		dstep = 1.0 / *drate;
 		fprintf(stderr, "using factor of %ld/%ld, working at %.2fhz\n", factor_L, factor_M, *drate);
-		cnt_amp = malloc(sizeof(float) * factor_M);
-		dat_amp = malloc(sizeof(float) * factor_M);
 		cnt_q = malloc(sizeof(float complex) * factor_L);
 		dat_q = malloc(sizeof(float complex) * factor_L);
 		// same factor to keep life simple and have accurate horizontal sync
-		cnt_ddc = alloc_ddc(1200.0, 200.0, step, cnt_taps, factor_L, factor_M, kaiser, 2.0);
+		cnt_ddc = alloc_ddc(factor_L, factor_M, 1200.0, 200.0, rate, cnt_taps, kaiser, 2.0);
-		dat_ddc = alloc_ddc(1900.0, 800.0, step, dat_taps, factor_L, factor_M, kaiser, 2.0);
+		dat_ddc = alloc_ddc(factor_L, factor_M, 1900.0, 800.0, rate, dat_taps, kaiser, 2.0);
 		// delay input by phase shift of other filter to synchronize outputs
-		cnt_delay = alloc_delay((dat_taps - 1) / (2 * factor_L));
+		cnt_delay = (dat_taps - 1) / (2 * factor_L);
-		dat_delay = alloc_delay((cnt_taps - 1) / (2 * factor_L));
+		dat_delay = (cnt_taps - 1) / (2 * factor_L);
 
-		buff = (short *)malloc(sizeof(short) * channels * factor_M);
+		pcm_buff = (short *)malloc(sizeof(short) * channels * factor_M);
+
+		// 0.1 second history + enough room for delay and taps
+		buffer = alloc_buffer(0.1 * rate + 2 * fmaxf(cnt_delay, dat_delay));
 
 		// start immediately below
 		out = factor_L;
@@ -350,22 +348,20 @@ int demodulate(struct pcm *pcm, float *cnt_freq, float *dat_freq, float *drate)
 
 	if (out >= factor_L) {
 		out = 0;
-		if (!read_pcm(pcm, buff, factor_M)) {
+		if (!read_pcm(pcm, pcm_buff, factor_M)) {
 			init = 0;
-			free(buff);
+			free(pcm_buff);
 			free_ddc(cnt_ddc);
 			free_ddc(dat_ddc);
-			free(cnt_amp);
+			free_buffer(buffer);
-			free(dat_amp);
 			return 0;
 		}
-		for (int j = 0; j < factor_M; j++) {
+		float *buff = 0;
-			float amp = (float)buff[j * channels] / 32767.0;
+		for (int j = 0; j < factor_M; j++)
-			cnt_amp[j] = do_delay(cnt_delay, amp);
+			buff = do_buffer(buffer, (float)pcm_buff[j * channels] / 32767.0);
-			dat_amp[j] = do_delay(dat_delay, amp);
+
-		}
+		do_ddc(cnt_ddc, buff + cnt_delay, cnt_q);
-		do_ddc(cnt_ddc, cnt_amp, cnt_q);
+		do_ddc(dat_ddc, buff + dat_delay, dat_q);
-		do_ddc(dat_ddc, dat_amp, dat_q);
 	}
 
 	*cnt_freq = fclampf(1200.0 + cargf(cnt_q[out] * conjf(cnt_last)) / (2.0 * M_PI * dstep), 1100.0, 1300.0);