commit 54f81882da338cfe10b71aa924dd0df1d73741a2
parent 9eee0b23d1b0004f99990ba8e58f2eeb8e91c32c
Author: Egor Achkasov <eaachkasov@edu.hse.ru>
Date: Thu, 9 Jan 2025 17:04:47 +0100
Merge branch 'master' into ara
Diffstat:
7 files changed, 71 insertions(+), 80 deletions(-)
diff --git a/GNUmakefile b/GNUmakefile
@@ -9,12 +9,12 @@ OBJ += $(OBJDIR)/main.o
TARGET = $(BUILDDIR)/mmse
+all: $(TARGET)
+
# Debug target
dbg: CFLAGS += -g
dbg: $(TARGET)
-all: $(TARGET)
-
# Compile the elf
$(TARGET): $(OBJ)
@mkdir -p $(BUILDDIR)
diff --git a/include/mmserv.h b/include/mmserv.h
@@ -11,8 +11,8 @@
*/
typedef struct {
- data_t re;
- data_t im;
+ data_t re;
+ data_t im;
} complex;
/** Calculate MMSE estimation of x in y = H*x + n
diff --git a/script/gen_data.py b/script/gen_data.py
@@ -1,6 +1,6 @@
#!/usr/bin/python
-from util import read_defines, interleave
+from util import read_defines
import numpy as np
from numpy.random import random, normal
@@ -35,16 +35,15 @@ n = normal(0, NOISE_STD_DEVIATION, (NUM_RX_ANT, NUM_SC)) \
+ 1.j * normal(0, NOISE_STD_DEVIATION, (NUM_RX_ANT, NUM_SC))
# Received signal
y = np.einsum("ijk,jk->ik", H, x) + n
-
# Noise covariance matrix
R = np.eye(NUM_TX_ANT, NUM_TX_ANT, dtype=np.complex64) * NOISE_STD_DEVIATION**2
-# Cast the complex data to int16
-x_data = interleave(x)
-H_data = interleave(H)
-R_data = interleave(R)
-y_data = interleave(y)
-data_tuple = (x_data, H_data, R_data, y_data)
+data_tuple = (
+ x.real.astype(np.float32), x.imag.astype(np.float32),
+ H.real.astype(np.float32), H.imag.astype(np.float32),
+ R.real.astype(np.float32), R.imag.astype(np.float32),
+ y.real.astype(np.float32), y.imag.astype(np.float32)
+)
class Section:
@@ -58,10 +57,14 @@ class Section:
sections = [
- Section("x", "data/x.txt", "3", 32, x.size * 2),
- Section("H", "data/H.txt", "3", 32, H.size * 2),
- Section("R", "data/R.txt", "3", 32, R.size * 2),
- Section("y", "data/y.txt", "3", 32, y.size * 2),
+ Section("x_re", "data/x_re.txt", "3", 32, x.size),
+ Section("x_im", "data/x_im.txt", "3", 32, x.size),
+ Section("H_re", "data/H_re.txt", "3", 32, H.size),
+ Section("H_im", "data/H_im.txt", "3", 32, H.size),
+ Section("R_re", "data/R_re.txt", "3", 32, R.size),
+ Section("R_im", "data/R_im.txt", "3", 32, R.size),
+ Section("y_re", "data/y_re.txt", "3", 32, y.size),
+ Section("y_im", "data/y_im.txt", "3", 32, y.size),
]
# Create "data" directory if it does not exist
@@ -93,4 +96,3 @@ if WRITE_DATA_S:
for n in range(4):
s += "%02x" % bs[i+3-n]
print(" .word 0x%s" % s)
-
diff --git a/script/mmse.py b/script/mmse.py
@@ -1,41 +1,35 @@
#!/usr/bin/python
-from util import read_defines, deinterleave, interleave
+""" Python implementation of the MMSE receiver with Cholesky decomposition """
+
+from util import read_defines, load_xHRy
import numpy as np
from os import path, makedirs
-NUM_RX_ANT, NUM_TX_ANT, NUM_SC = read_defines()
-
# Read the data from the files
-x = deinterleave(np.fromfile("data/x.bin", dtype=np.float32), (NUM_TX_ANT, NUM_SC))
-H = deinterleave(np.fromfile("data/H.bin", dtype=np.float32), (NUM_RX_ANT, NUM_TX_ANT, NUM_SC))
-R = deinterleave(np.fromfile("data/R.bin", dtype=np.float32), (NUM_TX_ANT, NUM_TX_ANT, NUM_SC))
-y = deinterleave(np.fromfile("data/y.bin", dtype=np.float32), (NUM_RX_ANT, NUM_SC))
-
+NUM_RX_ANT, NUM_TX_ANT, NUM_SC = read_defines()
+x, H, R, y = load_xHRy(NUM_RX_ANT, NUM_TX_ANT, NUM_SC)
x_mmse = np.empty((NUM_TX_ANT, NUM_SC), np.complex64)
for sc in range(NUM_SC):
HH = H[..., sc].conj().T
L = np.linalg.cholesky(HH @ H[..., sc] + R[..., sc])
HHy = np.einsum("ij,j->i", HH, y[:, sc])
- print("\nHHy:\n", HHy)
- print("\nL:\n", L)
z = np.empty((NUM_TX_ANT,), dtype=np.complex64)
for i in range(NUM_TX_ANT):
z[i] = (HHy[i] - np.sum([L[i, j]*z[j] for j in range(i)])) / L[i, i]
- print("\nz:\n", z)
LH = L.conj().T
for i in range(NUM_TX_ANT-1, -1, -1):
x_mmse[i, sc] = (z[i] - np.sum([LH[i, j]*x_mmse[j, sc] for j in range(i+1, NUM_TX_ANT)], 0)) / LH[i, i]
- print("\nx_mmse:\n", x_mmse)
# Output the data
-x_mmse_data = interleave(x_mmse)
out_dir = path.join(path.dirname(__file__), "..", "out")
if not path.exists(out_dir):
makedirs(out_dir)
-with open(path.join(out_dir, "x_mmse_python.bin"), "wb") as f:
- f.write(x_mmse_data.tobytes())
+with open(path.join(out_dir, "x_mmse_python_re.bin"), "wb") as f:
+ f.write(x_mmse.real.astype(np.float32).tobytes())
+with open(path.join(out_dir, "x_mmse_python_im.bin"), "wb") as f:
+ f.write(x_mmse.imag.astype(np.float32).tobytes())
diff --git a/script/mmse_nosqrt.py b/script/mmse_nosqrt.py
@@ -1,21 +1,17 @@
#!/usr/bin/python
-from util import read_defines, deinterleave, interleave
+from util import read_defines, load_xHRy
import numpy as np
#from scipy.linalg import ldl
from os import path, makedirs
-NUM_RX_ANT, NUM_TX_ANT, NUM_SC = read_defines()
-
# Read the data from the files
-x = deinterleave(np.fromfile("data/x.bin", dtype=np.float32), (NUM_TX_ANT, NUM_SC))
-H = deinterleave(np.fromfile("data/H.bin", dtype=np.float32), (NUM_RX_ANT, NUM_TX_ANT, NUM_SC))
-R = deinterleave(np.fromfile("data/R.bin", dtype=np.float32), (NUM_TX_ANT, NUM_TX_ANT, NUM_SC))
-y = deinterleave(np.fromfile("data/y.bin", dtype=np.float32), (NUM_RX_ANT, NUM_SC))
-
+NUM_RX_ANT, NUM_TX_ANT, NUM_SC = read_defines()
+x, H, R, y = load_xHRy(NUM_RX_ANT, NUM_TX_ANT, NUM_SC)
x_mmse = np.empty((NUM_TX_ANT, NUM_SC), np.complex64)
+
def ldl(A):
n = A.shape[0]
D = np.zeros_like(A)
@@ -37,24 +33,20 @@ for sc in range(NUM_SC):
HH = H[..., sc].conj().T
L, D = ldl(HH @ H[..., sc] + R[..., sc])
HHy = np.einsum("ij,j->i", HH, y[:, sc])
- print("\nHHy:\n", HHy)
- print("\nL:\n", L)
- print("\nD:\n", D)
z = np.empty((NUM_TX_ANT,), dtype=np.complex64)
for i in range(NUM_TX_ANT):
z[i] = (HHy[i] - np.sum([L[i, j]*z[j] for j in range(i)]))
- print("\nz:\n", z)
LH = L.conj().T
for i in range(NUM_TX_ANT-1, -1, -1):
x_mmse[i, sc] = z[i] / D[i, i] - np.sum([L[i, j] * x[j, sc] for j in range(i+1, NUM_TX_ANT)])
- print("\nx_mmse:\n", x_mmse)
# Output the data
-x_mmse_data = interleave(x_mmse)
out_dir = path.join(path.dirname(__file__), "..", "out")
if not path.exists(out_dir):
makedirs(out_dir)
-with open(path.join(out_dir, "x_mmse_python.bin"), "wb") as f:
- f.write(x_mmse_data.tobytes())
+with open(path.join(out_dir, "x_mmse_nosqrt_python_re.bin"), "wb") as f:
+ f.write(x_mmse.real.astype(np.float32).tobytes())
+with open(path.join(out_dir, "x_mmse_nosqrt_python_im.bin"), "wb") as f:
+ f.write(x_mmse.imag.astype(np.float32).tobytes())
diff --git a/script/util.py b/script/util.py
@@ -28,32 +28,26 @@ def read_defines():
return NUM_RX_ANT, NUM_TX_ANT, NUM_SC
-def interleave(data: np.ndarray) -> np.ndarray:
- """Cast a np.complex64 array to np.float32 array.
-
- Args:
- data (np.ndarray): The complex data to be casted. Dtype: np.complex64
-
- Returns:
- np.ndarray: The interleaved data. Dtype: np.float32
- """
- res = np.empty(2 * data.size, dtype=np.float32)
- res[0::2] = data.flatten().real.astype(np.float32)
- res[1::2] = data.flatten().imag.astype(np.float32)
- return res
-
-
-def deinterleave(data: np.ndarray, shape: tuple) -> np.ndarray:
- """Cast a np.float32 array to np.complex64 array.
-
- Args:
- data (np.ndarray): The interleaved data to be casted. Dtype: np.float32
- shape (tuple): The shape of the deinterleaved data.
+def load_xHRy(NUM_RX_ANT, NUM_TX_ANT, NUM_SC) -> tuple:
+ """Load x, H, R and y from the data files.
+ Assumes that the following files are present in the data directory:
+ x_re.bin, x_im.bin, H_re.bin, H_im.bin, R_re.bin, R_im.bin, y_re.bin, y_im.bin
+ and they contain float32 data of the appropriate sizes.
Returns:
- np.ndarray: The deinterleaved data. Dtype: np.complex64
+ tuple: x, H, R, y
"""
- res = np.empty(shape, np.complex64)
- res.real = data[0::2].reshape(shape)
- res.imag = data[1::2].reshape(shape)
- return res
+ x = np.fromfile("data/x_re.bin", dtype=np.float32)
+ x = x + 1j*np.fromfile("data/x_im.bin", dtype=np.float32)
+ x = x.reshape((NUM_TX_ANT, NUM_SC))
+ H = np.fromfile("data/H_re.bin", dtype=np.float32)
+ H = H + 1j*np.fromfile("data/H_im.bin", dtype=np.float32)
+ H = H.reshape((NUM_RX_ANT, NUM_TX_ANT, NUM_SC))
+ R = np.fromfile("data/R_re.bin", dtype=np.float32)
+ R = R + 1j*np.fromfile("data/R_im.bin", dtype=np.float32)
+ R = R.reshape((NUM_TX_ANT, NUM_TX_ANT, NUM_SC))
+ y = np.fromfile("data/y_re.bin", dtype=np.float32)
+ y = y + 1j*np.fromfile("data/y_im.bin", dtype=np.float32)
+ y = y.reshape((NUM_RX_ANT, NUM_SC))
+
+ return x, H, R, y
diff --git a/script/viz_compare_x.py b/script/viz_compare_x.py
@@ -1,6 +1,6 @@
#!/usr/bin/env python3
-from util import read_defines, deinterleave
+from util import read_defines
import numpy as np
from matplotlib import pyplot as plt
@@ -8,11 +8,15 @@ from matplotlib import pyplot as plt
# read `out/x_mmse.bin` and plot the complex signal data samples
NUM_RX_ANT, NUM_TX_ANT, NUM_SC = read_defines()
-
-with open("out/x_mmse.bin", "rb") as f:
- x_mmse = deinterleave(np.fromfile(f, dtype=np.float32), (NUM_TX_ANT, NUM_SC))
-with open("data/x.bin", "rb") as f:
- x = deinterleave(np.fromfile(f, dtype=np.float32), (NUM_TX_ANT, NUM_SC))
+x = np.fromfile("data/x_re.bin", dtype=np.float32) \
+ + 1j * np.fromfile("data/x_im.bin", dtype=np.float32)
+x = x.reshape((NUM_TX_ANT, NUM_SC))
+x_mmse = np.fromfile("out/x_mmse_re.bin", dtype=np.float32) \
+ + 1j * np.fromfile("out/x_mmse_im.bin", dtype=np.float32)
+x_mmse = x_mmse.reshape((NUM_TX_ANT, NUM_SC))
+x_mmse_python = np.fromfile("out/x_mmse_python_re.bin", dtype=np.float32) \
+ + 1j * np.fromfile("out/x_mmse_python_im.bin", dtype=np.float32)
+x_mmse_python = x_mmse_python.reshape((NUM_TX_ANT, NUM_SC))
for tx in range(NUM_TX_ANT):
plt.scatter(x_mmse[tx].real,
@@ -25,6 +29,11 @@ for tx in range(NUM_TX_ANT):
label=f"Orig Tx {tx}",
marker="o",
color=plt.cm.hsv(tx / NUM_TX_ANT))
+ plt.scatter(x_mmse_python[tx].real,
+ x_mmse_python[tx].imag,
+ label=f"MMSE Python Tx {tx}",
+ marker="^",
+ color=plt.cm.hsv(tx / NUM_TX_ANT))
plt.axhline(0, color='black')
plt.axvline(0, color='black')
plt.xlim(-1.1, 1.1)
