ogl_beamforming

Commit: bce389f0d683d84cf92632e9e2f314505b928d8d
Parent: dca23e416ee968820bb82665c120e15000143278
Author: Randy Palamar
Date:   Tue, 29 Oct 2024 09:46:14 -0600

support inverse hadamard matrices with dimension 2^n * 12

we need this for arrays with 48 columns

Diffstat:
M
beamformer.c
 | 
3
++-
M
util.c
 | 
55
+++++++++++++++++++++++++++++++++++++++++++++++++++++--

2 files changed, 55 insertions(+), 3 deletions(-)
diff --git a/beamformer.c b/beamformer.c
@@ -119,7 +119,8 @@ alloc_shader_storage(BeamformerCtx *ctx, Arena a)
 	cs->hadamard_dim       = (uv2){.x = dec_data_dim.z, .y = dec_data_dim.z};
 	size hadamard_elements = dec_data_dim.z * dec_data_dim.z;
 	i32  *hadamard         = alloc(&a, i32, hadamard_elements);
-	fill_hadamard(hadamard, dec_data_dim.z);
+	i32  *tmp              = alloc(&a, i32, hadamard_elements);
+	fill_hadamard_transpose(hadamard, tmp, dec_data_dim.z);
 	glDeleteBuffers(1, &cs->hadamard_ssbo);
 	glCreateBuffers(1, &cs->hadamard_ssbo);
 	glNamedBufferStorage(cs->hadamard_ssbo, hadamard_elements * sizeof(i32), hadamard, 0);
diff --git a/util.c b/util.c
@@ -1,4 +1,19 @@
 /* See LICENSE for license details. */
+static i32 hadamard_12_12_transpose[] = {
+	1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
+	1, -1, -1,  1, -1, -1, -1,  1,  1,  1, -1,  1,
+	1,  1, -1, -1,  1, -1, -1, -1,  1,  1,  1, -1,
+	1, -1,  1, -1, -1,  1, -1, -1, -1,  1,  1,  1,
+	1,  1, -1,  1, -1, -1,  1, -1, -1, -1,  1,  1,
+	1,  1,  1, -1,  1, -1, -1,  1, -1, -1, -1,  1,
+	1,  1,  1,  1, -1,  1, -1, -1,  1, -1, -1, -1,
+	1, -1,  1,  1,  1, -1,  1, -1, -1,  1, -1, -1,
+	1, -1, -1,  1,  1,  1, -1,  1, -1, -1,  1, -1,
+	1, -1, -1, -1,  1,  1,  1, -1,  1, -1, -1,  1,
+	1,  1, -1, -1, -1,  1,  1,  1, -1,  1, -1, -1,
+	1, -1,  1, -1, -1, -1,  1,  1,  1, -1,  1, -1,
+};
+
 static void *
 mem_clear(u8 *p, u8 c, size len)
 {
@@ -258,9 +273,41 @@ parse_f64(s8 s)
 }
 
 static void
-fill_hadamard(i32 *m, u32 dim)
+fill_kronecker_sub_matrix(i32 *out, i32 out_stride, i32 scale, i32 *b, uv2 b_dim)
+{
+	for (u32 i = 0; i < b_dim.y; i++)
+		for (u32 j = 0; j < b_dim.x; j++)
+			out[i * out_stride + j] = scale * b[i * b_dim.x + j];
+}
+
+/* NOTE: this won't check for valid space/etc and assumes row major order */
+static void
+kronecker_product(i32 *out, i32 *a, uv2 a_dim, i32 *b, uv2 b_dim)
 {
-	ASSERT(dim && ISPOWEROF2(dim));
+	uv2 out_dim = {.x = a_dim.x * b_dim.x, .y = a_dim.y * b_dim.y};
+	for (u32 i = 0; i < a_dim.y; i++) {
+		for (u32 j = 0; j < a_dim.x; j++) {
+			fill_kronecker_sub_matrix(out + j * b_dim.y + i * out_dim.y * b_dim.x,
+			                          out_dim.y, a[i * a_dim.x + j], b, b_dim);
+		}
+	}
+}
+
+/* NOTE/TODO: to support even more hadamard sizes use the Paley construction */
+static void
+fill_hadamard_transpose(i32 *out, i32 *tmp, u32 dim)
+{
+	ASSERT(dim);
+	b32 power_of_2 = ISPOWEROF2(dim);
+
+	if (!power_of_2) {
+		ASSERT(dim % 12 == 0);
+		dim /= 12;
+	}
+
+	i32 *m;
+	if (power_of_2) m = out;
+	else            m = tmp;
 
 	#define IND(i, j) ((i) * dim + (j))
 	m[0] = 1;
@@ -275,4 +322,8 @@ fill_hadamard(i32 *m, u32 dim)
 		}
 	}
 	#undef IND
+
+	if (!power_of_2)
+		kronecker_product(out, tmp, (uv2){.x = dim, .y = dim}, hadamard_12_12_transpose,
+		                  (uv2){.x = 12, .y = 12});
 }