MUSIC eigenvalue visualisation

One of the issues I've been facing was that the range of the eigenvalues
was quite large, and I couldn't easily see a natural threshold to set to
separate the signal and noise subspaces.

This visualisation provides a histogram for this scenario, to assist in
choosing the parameter (which is different in different configurations).

where_fi/cli/__init__.py

@@ -57,8 +57,7 @@ def heatmap() -> None:
         logger.info(f"Processed CSI data with shape {processed.shape}")
         processed_tensor = torch.tensor(processed, device=device)
         aoa.update(processed_tensor)
-        heatmap = aoa.heatmap()
-        visualise_data(heatmap, visualise.DataType.HEATMAP)
+        aoa.heatmap(visualiser=visualise_data)
 
     globals.csi_producer(csi_callback=callback)
     logger.info("Finished processing CSI data")

where_fi/processing/aoa.py

@@ -1,11 +1,13 @@
 import logging
 from datetime import datetime
+from typing import Any, Callable
 
 import numpy as np
 import numpy.typing as npt
 import torch
 
 from ..config import config
+from ..visualise import server as visualise
 
 logger = logging.getLogger(__name__)
 
@@ -111,16 +113,24 @@ class AoA:
         assert steering.shape == (N, self.N_rx // 2, self.N_subcarriers // 2)
         return steering.reshape(N, -1)
 
-    def evaluate(self, theta: torch.Tensor, tof: torch.Tensor) -> torch.Tensor:
+    def evaluate(
+        self,
+        theta: torch.Tensor,
+        tof: torch.Tensor,
+        visualiser: None
+        | Callable[[npt.NDArray[Any], visualise.DataType], None] = None,
+    ) -> torch.Tensor:
         R = torch.mean(self.historical_autocorr, dim=0)
 
         # The smallest eigenvectors span the noise subspace,
         # and the largest span the signal subspace.
         logger.debug(f"Calculating eigenvectors of R: {R.shape}")
         eigvals, eigvecs = torch.linalg.eig(R)
+        if visualiser:
+            visualiser(eigvals.numpy(), visualise.DataType.EIGENVALUES)
         assert isinstance(eigvals, torch.Tensor)
         assert isinstance(eigvecs, torch.Tensor)
-        logger.info(f"Eigenvalues: {eigvals}")
+        logger.debug(f"Eigenvalues: {eigvals}")
         E_n = eigvecs[:, torch.abs(eigvals) < config.music.eigval_threshold]
 
         logger.debug(f"Signal subspace: {E_n.shape}")
@@ -136,7 +146,10 @@ class AoA:
         c: torch.Tensor = 1 / (steering_h @ E_n @ E_n_H @ steering)
         return torch.abs(c)[:, 0, 0]
 
-    def heatmap(self) -> npt.NDArray[np.float32]:
+    def heatmap(
+        self,
+        visualiser: None | Callable[[npt.NDArray[Any], visualise.DataType], None],
+    ) -> npt.NDArray[np.float32]:
         thetas = np.linspace(
             0, np.pi, config.music.heatmap.theta_resolution, dtype=np.float32
         )
@@ -153,12 +166,16 @@ class AoA:
         evaluated = self.evaluate(
             torch.tensor(thetas_mesh.reshape(-1)),
             torch.tensor(tofs_mesh.reshape(-1)),
+            visualiser=visualiser,
         )
         logger.debug(f"Evaluated heatmap: {evaluated.shape}")
         heatmap: npt.NDArray[np.float32] = evaluated.reshape(
             config.music.heatmap.tof_resolution,
             config.music.heatmap.theta_resolution,
         ).numpy(force=True)
+
+        if visualiser:
+            visualiser(heatmap, visualise.DataType.HEATMAP)
         return heatmap
 
 

where_fi/visualise/frontend/src/components/Figure.vue

@@ -42,6 +42,24 @@ const plotData = computed(() => {
             type: 'heatmap' as const,
           },
         ]
+      case 'histogram':
+        const bin_start = data.value.figure.histogram.bins.slice(0, -1)
+        const bin_end = data.value.figure.histogram.bins.slice(1)
+        const bin_center = bin_start.map((start, i) => (start + bin_end[i]) / 2)
+        const bin_width = bin_start.map((start, i) => bin_end[i] - start)
+        console.log('Histogram data:', bin_center, bin_width)
+        console.log(
+          'Sizes:',
+          bin_center.length,
+          bin_width.length,
+          data.value.figure.histogram.data.length,
+        )
+        return data.value.figure.histogram.data.map((series) => ({
+          x: bin_center,
+          y: series.data,
+          width: bin_width,
+          type: 'bar' as const,
+        }))
     }
   }
   return []
@@ -49,8 +67,8 @@ const plotData = computed(() => {
 
 const layout = {
   title: { text: figure.title },
-  xaxis: { title: { text: figure.xLabel } },
-  yaxis: { title: { text: figure.yLabel } },
+  xaxis: { title: { text: figure.xLabel }, type: figure.logx ? 'log' : undefined },
+  yaxis: { title: { text: figure.yLabel }, type: figure.logy ? 'log' : undefined },
   height: 700,
 }
 

where_fi/visualise/protos/figure.proto

@@ -2,12 +2,16 @@ syntax = "proto3";
 
 import "figure_type/line.proto";
 import "figure_type/heatmap.proto";
+import "figure_type/histogram.proto";
 
 message Figure {
     string uuid = 1;
     string title = 2;
     string x_label = 3;
     string y_label = 4;
+
+    bool logx = 5;
+    bool logy = 6;
 }
 
 message FigureData {
@@ -15,6 +19,7 @@ message FigureData {
     oneof figure {
         LineChartData line = 2;
         HeatmapData heatmap = 3;
+        HistogramData histogram = 4;
     }
 }
 

where_fi/visualise/protos/figure_type/histogram.proto

@@ -0,0 +1,9 @@
+syntax = "proto3";
+
+message HistogramSeries { repeated float data = 1; }
+
+message HistogramData {
+    string uuid = 1;
+    repeated HistogramSeries data = 2;
+    repeated float bins = 3;
+}

where_fi/visualise/server/__init__.py

@@ -13,7 +13,7 @@ import numpy as np
 import numpy.typing as npt
 
 from .generated import figure_pb2, figure_pb2_grpc
-from .generated.figure_type import heatmap_pb2, line_pb2
+from .generated.figure_type import heatmap_pb2, histogram_pb2, line_pb2
 
 logger = logging.getLogger(__name__)
 
@@ -26,6 +26,7 @@ class DataType(Enum):
     UNWRAPPED_PHASE = 2
     PROCESSED_CSI = 3
     HEATMAP = 4
+    EIGENVALUES = 5
 
 
 @dataclass
@@ -47,6 +48,13 @@ figures = {
         x_label="Subcarrier",
         y_label="Unwrapped phase",
     ),
+    "aoa_eig": figure_pb2.Figure(
+        uuid=str(uuid.uuid4()),
+        title="AoA Eigenvalues",
+        x_label="Eigenvalue",
+        y_label="Frequency",
+        logx=True,
+    ),
     "processed_phase": figure_pb2.Figure(
         uuid=str(uuid.uuid4()),
         title="Preprocessed CSI Phase",
@@ -130,6 +138,20 @@ def add_data(dtype: DataType, new_data: npt.NDArray[np.complex128]) -> None:
             linechart = line_pb2.LineChartData(lines=lines)
             for q in clients.get(uuid, []):
                 q.put(figure_pb2.FigureData(uuid=uuid, line=linechart))
+        case DataType.EIGENVALUES:
+            uuid = figures["aoa_eig"].uuid
+            magn = np.abs(new_data)
+            bins = np.logspace(np.log10(magn.min()), np.log10(magn.max()), 10)
+            hist, _ = np.histogram(magn, bins=bins)
+            logger.info(f"bins: {bins}, hist: {hist}")
+            logger.info(
+                f"Updating histogram with {len(bins)} bins, and {len(hist)} bars"
+            )
+            logger.info(bins)
+            series = [histogram_pb2.HistogramSeries(data=hist)]
+            histogram = histogram_pb2.HistogramData(data=series, bins=bins)
+            for q in clients.get(uuid, []):
+                q.put(figure_pb2.FigureData(uuid=uuid, histogram=histogram))
 
 
 def listen_for_data(data_queue: "mp.Queue[VisualiserData]") -> None: