diff --git a/where_fi/cli/globals.py b/where_fi/cli/globals.py index d487bba..b298a3b 100644 --- a/where_fi/cli/globals.py +++ b/where_fi/cli/globals.py @@ -1,15 +1,20 @@ from pathlib import Path from .. import collection -from ..collection import file, ingest +from ..collection import file, ingest, raytracing csi_producer: collection.CSIProducer = collection.NoopCSIProducer() is_live = True -def main(from_file: Path | None = None) -> None: +def main(from_file: Path | None = None, from_environment: Path | None = None) -> None: global csi_producer, is_live + if from_file and from_environment: + raise ValueError("Cannot specify both a data file and an environment file") if from_file: csi_producer = file.FileCSIPRoducer(path=from_file) + elif from_environment: + environment = raytracing.Environment.from_config(from_environment) + csi_producer = raytracing.SimulatedCSIProducer(environment) else: csi_producer = ingest.RealtimeCSIProducer() diff --git a/where_fi/collection/__init__.py b/where_fi/collection/__init__.py index 8b54e48..43f6ec0 100644 --- a/where_fi/collection/__init__.py +++ b/where_fi/collection/__init__.py @@ -3,7 +3,7 @@ from typing import Iterator import numpy as np import numpy.typing as npt -from . import file, ingest +from . import file, ingest, raytracing from .protocols import CSIProducer, MergedCSI @@ -21,4 +21,4 @@ class NoopCSIProducer: CSIMatrix = npt.NDArray[np.complex64] -__all__ = ["file", "ingest", "NoopCSIProducer", "CSIProducer"] +__all__ = ["file", "ingest", "raytracing", "NoopCSIProducer", "CSIProducer"] diff --git a/where_fi/collection/raytracing/__init__.py b/where_fi/collection/raytracing/__init__.py new file mode 100644 index 0000000..c992209 --- /dev/null +++ b/where_fi/collection/raytracing/__init__.py @@ -0,0 +1,237 @@ +import logging +import time +from dataclasses import dataclass +from datetime import datetime +from pathlib import Path +from typing import Iterator + +import numpy as np +import numpy.typing as npt + +from where_fi.collection.protocols import CSIProducer, MergedCSI +from where_fi.config import config + +C = 299_792_458 +MAX_DEPTH = 2 +LOSS_EXPONENT = 0.8 + +logger = logging.getLogger(__name__) + + +@dataclass +class PathComponent: + delay: float + phase: float + attenuation: float + frequency: float + + +class ChannelImpulseResponse: + def __init__(self, path_components: list[PathComponent]) -> None: + self.path_components = path_components + + @staticmethod + def delayed( + other: "ChannelImpulseResponse", delay: float, reflect: bool = False + ) -> "ChannelImpulseResponse": + new_path_components: list[PathComponent] = [] + distance = delay * C + + for component in other.path_components: + wavelength = C / component.frequency + new_path_components.append( + PathComponent( + delay=component.delay + delay, + phase=component.phase + + 2 * np.pi * component.frequency * delay + + (np.pi if reflect else 0), + attenuation=component.attenuation + * (np.exp(-distance * LOSS_EXPONENT)), + frequency=component.frequency, + ) + ) + return ChannelImpulseResponse(new_path_components) + + def __add__(self, other: "ChannelImpulseResponse") -> "ChannelImpulseResponse": + return ChannelImpulseResponse(self.path_components + other.path_components) + + +class Object: + def __init__(self, x: float, y: float) -> None: + self.x = x + self.y = y + self.cir = ChannelImpulseResponse([]) + + def distance(self, other: "Object") -> float: + return ((self.x - other.x) ** 2 + (self.y - other.y) ** 2) ** 0.5 + + def reflect( + self, + incoming_cir: ChannelImpulseResponse, + target: list["Object"], + depth: int = 0, + ) -> None: + self.cir = self.cir = self.cir + incoming_cir + + if depth > MAX_DEPTH: + return + + for obj in target: + if id(obj) == id(self): + continue + distance = self.distance(obj) + assert distance > 0 + delay = distance / C + obj.reflect( + ChannelImpulseResponse.delayed(incoming_cir, delay, reflect=True), + target, + depth + 1, + ) + + +class Transmitter(Object): + def __init__(self, x: float, y: float) -> None: + super().__init__(x, y) + + +DELTA_T = 10 +GAMMA = np.pi / 4 + + +class Receiver(Object): + def __init__(self, x: float, y: float, ideal: bool = True) -> None: + super().__init__(x, y) + self.delta_t = 0 if ideal else DELTA_T + self.gamma = 0 if ideal else GAMMA + + def get_cfr(self) -> npt.NDArray[np.complex64]: + """ + Calculate the CSI matrix for this rx-tx pair. This is computed by the Fourier + Transform of the Channel Impulse Response (CIR). The CIR is calculated as in + [1], [2]. + + To get the FT of the CIR, we use the sifting property of the Dirac delta + + [1] - https://tns.thss.tsinghua.edu.cn/wst/docs/pre + [2] - https://dl.acm.org/doi/10.1145/2543581.2543592, Equation 4 + """ + cfr = np.zeros(len(config.subcarrier_frequencies), dtype=np.complex64) + for i_sub, subcarrier in enumerate(config.subcarrier_frequencies): + cfr[i_sub] = sum( + [ + component.attenuation + * np.exp(1j * component.phase) + * np.exp(-1j * subcarrier * component.delay) + for component in self.cir.path_components + ] + ) * np.exp( + 1j + * ( + 2 + * np.pi + * (i_sub / len(config.subcarrier_frequencies)) + * self.delta_t + + self.gamma + ) + ) + return cfr + + +class Environment: + def __init__(self, objects: list[Object]) -> None: + self.transmitters = [obj for obj in objects if isinstance(obj, Transmitter)] + self.receivers = [obj for obj in objects if isinstance(obj, Receiver)] + self.objects = [ + obj for obj in objects if obj not in self.transmitters + self.receivers + ] + + @staticmethod + def from_config(filename: Path) -> "Environment": + """ + Read a config file that includes a scene description and create an envionment + based on that + + The config file should be a text file where each line corresponds to an object. + The first word of each line should be the type of object (TX/RX/OBJ), followed + by the x and y coordinates of the object. + """ + objects: list[Object] = [] + with open(filename, "r") as f: + for line in f.readlines(): + if line.startswith("#"): + continue + parts = line.split(" ") + x, y = map(float, parts[1:]) + if parts[0] == "TX": + objects.append(Transmitter(x, y)) + elif parts[0] == "RX": + objects.append(Receiver(x, y)) + else: + objects.append(Object(x, y)) + + return Environment(objects) + + def get_csi(self) -> npt.NDArray[np.complex64]: + """ + Calculate the Channel State Information (CSI) matrix for the simulated + environment. + + The returned matrix is of shape (num_subcarriers, num_receivers, + num_transmitters) + + This is calculated by finding all paths leading to each receiver, and + calculating the CFR evaluated at each subcarrier. + """ + csi = np.zeros( + ( + len(config.subcarrier_frequencies), + len(self.receivers), + len(self.transmitters), + ), + dtype=np.complex64, + ) + for i_tx, transmitter in enumerate(self.transmitters): + for obj in self.objects + self.receivers + self.transmitters: + obj.cir = ChannelImpulseResponse([]) + transmitter.reflect( + ChannelImpulseResponse( + [ + PathComponent( + delay=0, + phase=0, + # attenuation=100000000, + attenuation=100, + frequency=subcarrier, + ) + for subcarrier in config.subcarrier_frequencies + ] + ), + self.objects + self.receivers, + ) + for i_rx, receiver in enumerate(self.receivers): + logger.info(f"RX {i_rx} paths: {len(receiver.cir.path_components)}") + csi[:, i_rx, i_tx] = receiver.get_cfr() + return csi + + +class SimulatedCSIProducer(CSIProducer): + pass + + def __init__(self, environment: Environment) -> None: + self.environment = environment + + def __call__(self) -> Iterator[MergedCSI]: + while True: + start = datetime.now() + csi = self.environment.get_csi() + yield MergedCSI( + frames={}, + matrix=csi, + ) + time.sleep( + max( + 0, + 1 / config.collection_sample_rate + - (datetime.now() - start).total_seconds(), + ) + ) diff --git a/where_fi/config/models.py b/where_fi/config/models.py index df86b24..179e06b 100644 --- a/where_fi/config/models.py +++ b/where_fi/config/models.py @@ -77,10 +77,22 @@ class Config(BaseModel): return int((self.channel_width * 1e6) // self.delta_f) - 8 @property - def delta_f(self) -> int: + def _delta_f_no_skipping(self) -> int: if self.frame_format == "HESU": - return 78_125 * self.preprocessing.subcarrier_step - return 312_500 * self.preprocessing.subcarrier_step + return 78_125 + return 312_500 + + @property + def delta_f(self) -> int: + return self._delta_f_no_skipping * self.preprocessing.subcarrier_step + + @property + def subcarrier_frequencies(self) -> list[int]: + num_subcarriers = self.channel_width * 1_000_000 // self._delta_f_no_skipping + return [ + self.central_freq_hz + i * self._delta_f_no_skipping + for i in range(-num_subcarriers // 2, num_subcarriers // 2) + ] @model_validator(mode="after") def channels(self) -> Self: