Improve MUSIC implementation

pyproject.toml

@@ -4,3 +4,6 @@ line-length = 88
 [tool.pyright]
 typeCheckingMode = "strict"
 reportMissingTypeStubs = "warning"
+
+[tool.pytest.ini_options]
+python_files = "*.py"

src/aoa.py

@@ -50,6 +50,8 @@ class AoA:
         H_sm = self.smooth(data)
 
         auto_corr = np.matmul(H_sm, np.conj(H_sm).T)
+        # This matrix is by definition Hermitian.
+        # Therefore, all of its eigenvectors are orthogonal.
 
         if self.historical_autocorr.size == 0:
             self.historical_autocorr = np.expand_dims(auto_corr, 0)
@@ -62,15 +64,14 @@ class AoA:
         if self.historical_autocorr.shape[0] > WINDOW_SIZE:
             self.historical_autocorr = self.historical_autocorr[-WINDOW_SIZE:]
 
+        # Is the moving average also Hermitian?
         R = np.mean(self.historical_autocorr, axis=0)
-        RR_h = np.matmul(R, np.conj(R).T)
-        # This matrix is by definition Hermitian.
-        # Therefore, all of its eigenvectors are orthogonal.
 
         # The smallest eigenvectors span the noise subspace,
         # and the largest span the signal subspace.
-        eigvals, eigvecs = np.linalg.eig(RR_h)
-        self.E_n = eigvecs[:, eigvals < config.EIGVAL_THRESHOLD]
+        eigvals, eigvecs = np.linalg.eigh(R)
+        print(eigvals)
+        self.E_n = eigvecs[:, np.abs(eigvals) < config.EIGVAL_THRESHOLD]
 
     omega_base = np.exp(-2j * np.pi * config.DELTA_F)
     phi_base = np.exp(
@@ -89,9 +90,8 @@ class AoA:
         antenna_v = omega_t ** np.arange(self.N_subcarriers // 2)
         phis = phi_theta ** np.arange(self.N_rx)
         antenna_v = np.expand_dims(antenna_v, axis=-1)
-        phis = np.expand_dims(phis, axis=-2)
         steering = antenna_v * phis
-        return steering.reshape(-1)
+        return steering.T.reshape(-1)
 
     def evaluate(self, theta: float, tof: float):
         try:
@@ -101,9 +101,9 @@ class AoA:
             logger.exception(e)
             return 0
         E_n = self.E_n
-        E_n_H = np.conj(self.E_n).T
+        E_n_H = np.conj(E_n).T
         c = 1 / (0.001 + (steering_h @ E_n @ E_n_H @ steering))
-        return c.real
+        return np.abs(c.real)
 
 
 def test_smoothing():
@@ -118,3 +118,15 @@ def test_smoothing():
     print(expected)
     assert np.allclose(smoothed, expected)
     pass
+
+
+def test_steering_vector():
+    aoa = AoA()
+    aoa.N_subcarriers = 10
+    aoa.N_rx = 2
+    print(aoa.omega_base)
+    print(aoa.phi_base)
+    tau = 1
+    theta = 0
+    print(aoa.steering_vector(theta, tau))
+    assert False