diff --git a/openhrv/config.py b/openhrv/config.py index 023905e..dde8dd5 100644 --- a/openhrv/config.py +++ b/openhrv/config.py @@ -2,8 +2,10 @@ from math import ceil -HRV_MEAN_WINDOW: Final[float] = 15.0 # seconds IBI_MEDIAN_WINDOW: Final[int] = 11 # samples +EWMA_WEIGHT_CURRENT_SAMPLE: Final[float] = ( + 0.1 # set in range [0, 1], 1 being minimal smoothing, see https://en.wikipedia.org/wiki/Exponential_smoothing +) MIN_BREATHING_RATE: Final[float] = 4.0 # breaths per minute MAX_BREATHING_RATE: Final[float] = 7.0 # breaths per minute @@ -16,16 +18,12 @@ MIN_PLOT_IBI: Final[int] = 300 MAX_PLOT_IBI: Final[int] = 1500 - -HRV_BUFFER_SIZE: Final[int] = 10 # samples # IBI buffer must hold enough samples such that even if IBIs (on average) were # MIN_IBI long, there'd be enough samples to display for IBI_HISTORY_DURATION seconds. IBI_HISTORY_DURATION: Final[int] = 60 # seconds IBI_BUFFER_SIZE: Final[int] = ceil(IBI_HISTORY_DURATION / (MIN_IBI / 1000)) # samples -MEANHRV_HISTORY_DURATION: Final[int] = 120 # seconds -MEANHRV_BUFFER_SIZE: Final[int] = ceil( - MEANHRV_HISTORY_DURATION / (MIN_IBI / 1000) -) # samples +HRV_HISTORY_DURATION: Final[int] = 120 # seconds +HRV_BUFFER_SIZE: Final[int] = ceil(HRV_HISTORY_DURATION / (MIN_IBI / 1000)) # samples COMPATIBLE_SENSORS: Final[list[str]] = ["Polar", "Decathlon Dual HR"] diff --git a/openhrv/model.py b/openhrv/model.py index 6c7696a..48fc0af 100644 --- a/openhrv/model.py +++ b/openhrv/model.py @@ -7,7 +7,6 @@ from openhrv.utils import get_sensor_address, sign, NamedSignal from openhrv.config import ( tick_to_breathing_rate, - MEANHRV_BUFFER_SIZE, HRV_BUFFER_SIZE, IBI_BUFFER_SIZE, MAX_BREATHING_RATE, @@ -16,13 +15,13 @@ IBI_MEDIAN_WINDOW, MIN_HRV_TARGET, MAX_HRV_TARGET, - HRV_MEAN_WINDOW, + EWMA_WEIGHT_CURRENT_SAMPLE, ) class Model(QObject): ibis_buffer_update = Signal(NamedSignal) - mean_hrv_update = Signal(NamedSignal) + hrv_update = Signal(NamedSignal) addresses_update = Signal(NamedSignal) pacer_rate_update = Signal(NamedSignal) hrv_target_update = Signal(NamedSignal) @@ -35,14 +34,16 @@ def __init__(self): self.ibis_seconds: deque[float] = deque( map(float, range(-IBI_BUFFER_SIZE, 1)), IBI_BUFFER_SIZE ) - self.mean_hrv_buffer: deque[float] = deque( - [-1] * MEANHRV_BUFFER_SIZE, MEANHRV_BUFFER_SIZE + self.hrv_buffer: deque[float] = deque([-1] * HRV_BUFFER_SIZE, HRV_BUFFER_SIZE) + self.hrv_seconds: deque[float] = deque( + map(float, range(-HRV_BUFFER_SIZE, 1)), HRV_BUFFER_SIZE ) - self.mean_hrv_seconds: deque[float] = deque( - map(float, range(-MEANHRV_BUFFER_SIZE, 1)), MEANHRV_BUFFER_SIZE - ) - self._hrv_buffer: deque[int] = deque([-1] * HRV_BUFFER_SIZE, HRV_BUFFER_SIZE) + # Exponentially Weighted Moving Average: + # - https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average + # - https://en.wikipedia.org/wiki/Exponential_smoothing + # - http://nestedsoftware.com/2018/04/04/exponential-moving-average-on-streaming-data-4hhl.24876.html + self.ewma_hrv: float = 1.0 self.sensors: list[QBluetoothDeviceInfo] = [] self.breathing_rate: float = float(MAX_BREATHING_RATE) self.hrv_target: int = math.ceil((MIN_HRV_TARGET + MAX_HRV_TARGET) / 2) @@ -80,9 +81,9 @@ def update_sensors(self, sensors: list[QBluetoothDeviceInfo]): ) def validate_ibi(self, ibi: int) -> int: - validated_ibi = ibi + validated_ibi: int = ibi if ibi < MIN_IBI or ibi > MAX_IBI: - median_ibi = math.ceil( + median_ibi: int = math.ceil( statistics.median( islice( self.ibis_buffer, @@ -97,11 +98,20 @@ def validate_ibi(self, ibi: int) -> int: validated_ibi = MAX_IBI else: validated_ibi = median_ibi - print(f"Correcting invalid IBI: {ibi} to {validated_ibi}") + print(f"Correcting outlier IBI {ibi} to {validated_ibi}") return validated_ibi + def validate_hrv(self, hrv: int) -> int: + validated_hrv: int = hrv + if hrv > MAX_HRV_TARGET: + validated_hrv = min(math.ceil(self.ewma_hrv), MAX_HRV_TARGET) + print(f"Correcting outlier HRV {hrv} to {validated_hrv}") + + return validated_hrv + def compute_local_hrv(self): + """https://doi.org/10.1038/s41598-019-44201-7 (Figure 2)""" self._duration_current_phase += self.ibis_buffer[-1] # 1: IBI rises, -1: IBI falls, 0: IBI constant current_ibi_phase: int = sign(self.ibis_buffer[-1] - self.ibis_buffer[-2]) @@ -115,38 +125,21 @@ def compute_local_hrv(self): self.update_hrv_buffer(local_hrv) seconds_current_phase: float = self._duration_current_phase / 1000 - self.update_mean_hrv_seconds(seconds_current_phase) + self.update_hrv_seconds(seconds_current_phase) self._duration_current_phase = 0 self._last_ibi_extreme = current_ibi_extreme self._last_ibi_phase = current_ibi_phase - def update_hrv_buffer(self, local_hrv): - if not self._hrv_buffer[0]: - return # wait until buffer is full - threshold = max(map(abs, self._hrv_buffer)) * 4 - if local_hrv > threshold: - print(f"Correcting outlier HRV {local_hrv} to {threshold}") - local_hrv = threshold - self._hrv_buffer.append(local_hrv) - self.update_mean_hrv_buffer() - - def update_mean_hrv_buffer(self): - seconds = self.ibis_seconds - hrv_buffer_seconds = islice(seconds, len(seconds) - HRV_BUFFER_SIZE, None) - start_idx_mean_window = next( - ( - i - for i, second in enumerate(hrv_buffer_seconds) - if second >= -HRV_MEAN_WINDOW - ), - -1, + def update_hrv_buffer(self, local_hrv: int): + self.ewma_hrv = ( + EWMA_WEIGHT_CURRENT_SAMPLE * self.validate_hrv(local_hrv) + + (1 - EWMA_WEIGHT_CURRENT_SAMPLE) * self.ewma_hrv ) - self.mean_hrv_buffer.append( - statistics.mean(islice(self._hrv_buffer, start_idx_mean_window, None)) - ) - self.mean_hrv_update.emit( - NamedSignal("MeanHrv", (self.mean_hrv_seconds, self.mean_hrv_buffer)) + + self.hrv_buffer.append(self.ewma_hrv) + self.hrv_update.emit( + NamedSignal("HeartRateVariability", (self.hrv_seconds, self.hrv_buffer)) ) def update_ibis_seconds(self, seconds: float): @@ -155,8 +148,8 @@ def update_ibis_seconds(self, seconds: float): ) self.ibis_seconds.append(0.0) - def update_mean_hrv_seconds(self, seconds: float): - self.mean_hrv_seconds = deque( - [i - seconds for i in self.mean_hrv_seconds], MEANHRV_BUFFER_SIZE + def update_hrv_seconds(self, seconds: float): + self.hrv_seconds = deque( + [i - seconds for i in self.hrv_seconds], HRV_BUFFER_SIZE ) - self.mean_hrv_seconds.append(0.0) + self.hrv_seconds.append(0.0) diff --git a/openhrv/view.py b/openhrv/view.py index 6b3df4d..862bc24 100644 --- a/openhrv/view.py +++ b/openhrv/view.py @@ -28,7 +28,7 @@ from openhrv.model import Model from openhrv.config import ( breathing_rate_to_tick, - MEANHRV_HISTORY_DURATION, + HRV_HISTORY_DURATION, IBI_HISTORY_DURATION, MAX_BREATHING_RATE, MIN_BREATHING_RATE, @@ -162,7 +162,7 @@ def __init__(self, model: Model): self.model = model self.model.ibis_buffer_update.connect(self.plot_ibis) - self.model.mean_hrv_update.connect(self.plot_hrv) + self.model.hrv_update.connect(self.plot_hrv) self.model.addresses_update.connect(self.list_addresses) self.model.pacer_rate_update.connect(self.update_pacer_label) self.model.hrv_target_update.connect(self.update_hrv_target) @@ -194,7 +194,7 @@ def __init__(self, model: Model): self.model.addresses_update.connect(self.logger.write_to_file) self.model.pacer_rate_update.connect(self.logger.write_to_file) self.model.hrv_target_update.connect(self.logger.write_to_file) - self.model.mean_hrv_update.connect(self.logger.write_to_file) + self.model.hrv_update.connect(self.logger.write_to_file) self.signals.annotation.connect(self.logger.write_to_file) self.ibis_widget = XYSeriesWidget( @@ -211,13 +211,13 @@ def __init__(self, model: Model): self.ibis_widget.y_axis.setRange(MIN_PLOT_IBI, MAX_PLOT_IBI) self.hrv_widget = XYSeriesWidget( - self.model.mean_hrv_seconds, self.model.mean_hrv_buffer, WHITE + self.model.hrv_seconds, self.model.hrv_buffer, WHITE ) self.hrv_widget.x_axis.setTitleText("Seconds") # The time series displays only the samples within the last - # MEANHRV_HISTORY_DURATION seconds, - # even though there are more samples in self.model.mean_hrv_seconds. - self.hrv_widget.x_axis.setRange(-MEANHRV_HISTORY_DURATION, 0) + # HRV_HISTORY_DURATION seconds, + # even though there are more samples in self.model.hrv_seconds. + self.hrv_widget.x_axis.setRange(-HRV_HISTORY_DURATION, 0) self.hrv_widget.y_axis.setTitleText("HRV (msec)") self.hrv_widget.y_axis.setRange(0, self.model.hrv_target) colorgrad = QLinearGradient(0, 0, 0, 1) # horizontal gradient diff --git a/test/app.py b/test/app.py index 8afc80b..38cfc8e 100644 --- a/test/app.py +++ b/test/app.py @@ -58,7 +58,7 @@ def __init__(self): super().__init__() # Polar sensor emits a (package of) IBI(s) about every second. # Here we "emit" / simulate IBI(s) in quicker succession in order to push the rendering. - self.mean_ibi = 300 + self.mean_ibi = 900 self.timer = QTimer() self.timer.setInterval(self.mean_ibi) self.timer.timeout.connect(self.simulate_ibi) @@ -78,10 +78,16 @@ def simulate_ibi(self): # in a sinusoidal pattern around `mean_ibi`, # in a range of `range_ibi`. breathing_rate = 6 - range_ibi = 40 # HRV must settle at this value + range_ibi = 100 # without noise, HRV settles at this value ibi = self.mean_ibi + (range_ibi / 2) * math.sin( 2 * math.pi * breathing_rate / 60 * time.time() ) + # add noise spikes + if randint(1, 30) == 1: + if randint(1, 2) == 1: + ibi += 500 + else: + ibi -= 500 self.ibi_update.emit(ibi)