Coverage for src/susi/reduc/demodulation/mod_states.py: 68%
145 statements
« prev ^ index » next coverage.py v7.5.0, created at 2025-08-11 10:03 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4In this module the mod state detector can be found
6@author: hoelken, iglesias (changed algorithm on 20241106)
7"""
8from __future__ import annotations
10import numpy as np
11import os
12import matplotlib.pyplot as plt
15from ...base import InsufficientDataException, Globals, Logging
16from ...plot import matplot_backend
17from ...utils import Collections
19logger = Logging.get_logger()
22class ModStateDetector:
23 """
24 Class will analyse the header information of the complete batch to detect the state of each frame and to
25 detect missing frames.
26 """
28 #: Constant defines the expected minimal delta that occurs on a roll-over
29 ROLLOVER_DELTA = -100
30 #: Constant that defines the allowed deviation from the mean TS delta
31 TSTMP_MEAN_DEVIATION = 1.31
32 #: Constant that defines when we consider two frames lost
33 TSTMP_DOUBLE_MEAN_DEVIATION = 2.2
34 #: Constant that defines when to terminate the program due to data corruption
35 TSTMP_CRITICAL_MEAN_DEVIATION = 3.33
37 def __init__(self, config, batch):
38 #: The config data item
39 self.config = config
40 #: A loaded (and sorted) batch of all the frames to demodulate (only header information is needed)
41 self.batch = batch
42 #: holds the angles per index
43 self.angles = np.empty(0)
44 #: holds the time stamps per index
45 self.tstmps = np.empty(0)
46 #: holds the time stamp deltas
47 self.ts_deltas = np.empty(0)
48 #: holds the median of the time stamp deltas
49 self.tstmp_median = np.empty(0)
50 #: holds the state per index
51 self.states = list()
52 #: holds the integer total number of modulation states
53 self.no_mod_cycles = np.empty(0)
55 def analyze(self) -> ModStateDetector:
56 """
57 Executes the analysis
58 :return: self
59 """
60 self.__get_orig_angles()
61 self.__get_orig_tstmps()
62 self.__correct_missing_frames()
63 self.__get_ini_modstate()
64 self.__build_state_map()
65 return self
67 def __get_orig_angles(self) -> None:
68 # Fetches the (uncorrected) index of angles from the fits batch
69 self.angles = Collections.as_float_array(self.batch.header_field(self.config.spol.pmu_ang_field))
70 self.no_mod_cycles = len(self.angles) // Globals.MOD_CYCLE_FRAMES
71 if self.no_mod_cycles < self.config.spol.min_mod_cycles:
72 req_time = self.config.spol.min_mod_cycles * Globals.MOD_CYCLE_FRAMES * Globals.SAMP_S
73 logger.warning(
74 'Not enough mod cycles, need at least %s (%s s) but got %s',
75 self.config.spol.min_mod_cycles,
76 req_time,
77 self.no_mod_cycles,
78 )
79 raise InsufficientDataException('Not enough mod cycles found to demodulate the data.')
81 def __get_orig_tstmps(self) -> None:
82 # Fetches the timestamp of all frames
83 self.tstmps = Collections.as_float_array(
84 self.batch.header_field(self.config.base.timestamp_field), dtype=np.float64
85 )
86 self.ts_deltas = np.diff(self.tstmps)
87 self.tstmp_median = np.median(self.ts_deltas)
89 def __correct_missing_frames(self) -> None:
90 # Check for timestamps deltas that are too large
91 tota_frames_lost = np.sum(self.ts_deltas > self.tstmp_median * ModStateDetector.TSTMP_MEAN_DEVIATION)
92 if tota_frames_lost > 0:
93 logger.warning('Found %s frames lost out a total of %s', tota_frames_lost, len(self.ts_deltas))
94 block_gaps = []
95 for i in range(0, self.no_mod_cycles):
96 block_start = i * Globals.MOD_CYCLE_FRAMES
97 block_end = block_start + Globals.MOD_CYCLE_FRAMES
98 gap_id, gap_frames = self.__check_timestamps(block_start, block_end)
99 if gap_id is not None:
100 block_gaps.append([i, gap_id, gap_frames])
101 if len(block_gaps) == 0:
102 return None
103 logger.warning('Found %s mod cycle(s) with missing frames.', len(block_gaps))
104 self._delete_mod_cycles(block_gaps)
106 def __check_timestamps(self, start_id, end_id):
107 ts_deltas = self.ts_deltas[start_id:end_id]
108 if any(delta > self.tstmp_median * ModStateDetector.TSTMP_CRITICAL_MEAN_DEVIATION for delta in ts_deltas):
109 logger.critical(
110 'Found too many frames lost from "%s" to "%s".',
111 self.batch.file_by(start_id),
112 self.batch.file_by(end_id),
113 )
114 self.__plot_on_error()
115 raise InsufficientDataException('Too many frames lost. See log for details.')
116 gap_id = np.sort(np.argwhere(ts_deltas > self.tstmp_median * ModStateDetector.TSTMP_MEAN_DEVIATION).flatten())
117 if len(gap_id) > 1:
118 raise InsufficientDataException('Found more than one timestamp gap in a single mod cycle. Not supported.')
119 if len(gap_id) == 0:
120 return None, None
121 gap_id = gap_id[0]
122 gap_frames = round(ts_deltas[gap_id] / self.tstmp_median) - 1
123 return gap_id, gap_frames
125 def _delete_mod_cycles(self, block_gaps):
126 # Eliminate the mod cycles with missing frames to ensure a continuous modulation state sequence
127 for i, gap in enumerate(block_gaps):
128 # remove frames to ensure the full block with missing frames is removed
129 # so the next block starts at the same modulation state as the first block
130 mod_cyc_id = gap[0]
131 if mod_cyc_id == 0:
132 logger.debug('First mod cycle has missing frames. Removing it.')
133 gap_id = gap[1]
134 gap_frames = gap[2]
135 start_id = mod_cyc_id * Globals.MOD_CYCLE_FRAMES
136 no_frames_to_remove = Globals.MOD_CYCLE_FRAMES - gap_frames
137 end_id = start_id + no_frames_to_remove
138 logger.warning('Processing mod cycle #%s that has %s missing frame(s)', mod_cyc_id, gap_frames)
139 self.__remove_block(start_id, end_id)
140 # corrects the mod_cyc_id of the following blocks with gaps
141 for j in range(i + 1, len(block_gaps)):
142 block_gaps[j][0] -= 1
144 def __remove_block(self, start_id, end_id):
145 logger.warning('Removing %s files from %s to %s', end_id - start_id, start_id, end_id)
146 self.angles = np.delete(self.angles, range(start_id, end_id))
147 self.no_mod_cycles -= 1
148 for i in range(start_id, end_id):
149 logger.debug('\t#%s: %s [IGNORED]', i, self.batch.file_by(start_id))
150 self.batch.batch.pop(start_id)
152 def __get_ini_modstate(self) -> None:
153 # Determines the location of modulation state 0.
154 # Defined as the state after the 360->0 roll over of the mean angle.
155 # No frame lost is assumed !!
156 # uses only a multiple of Globals.MOD_CYCLE_FRAMES
157 angles = self.angles[0 : self.no_mod_cycles * Globals.MOD_CYCLE_FRAMES]
158 angles_avg = angles.reshape(-1, Globals.MOD_CYCLE_FRAMES).mean(axis=0)
159 angles_avg = np.append(angles_avg[-1], angles_avg)
160 ang_diff = np.diff(angles_avg)
161 logger.debug('self.no_mod_cycles %s', self.no_mod_cycles)
162 self.ini_modstate_roll = np.argmax(np.abs(ang_diff))
163 if self.ini_modstate_roll == 0:
164 ini_modstate = 0
165 else:
166 ini_modstate = Globals.MOD_CYCLE_FRAMES - self.ini_modstate_roll
167 logger.info('Initial modulation state found is: %s', ini_modstate)
169 def __build_state_map(self) -> None:
170 # Create an array with length/mod_cycle_frames blocks of 0,..,11
171 self.states = np.tile(
172 np.arange(Globals.MOD_CYCLE_FRAMES), np.int64(np.ceil(len(self.angles) / Globals.MOD_CYCLE_FRAMES))
173 )
174 # shift the first state according to the initial modulation state location
175 self.states = np.roll(self.states, self.ini_modstate_roll)
176 self.states = self.states[0 : len(self.angles)]
178 def __plot_on_error(self):
179 # TODO Move to src/susi/analyse ?
180 # Plots relevant header quantities to help debugging lost frames or state ambiguity issues
181 # timestamp vs frame (as input, not necessarily sorted)
182 if self.config.data.log_dir is None:
183 return
185 with matplot_backend('Agg'):
186 fig = plt.figure(figsize=(15, 7))
187 fig.suptitle('')
188 plt.plot(self.tstmps, '*')
189 plt.xlabel('Frame')
190 plt.xlabel('PMU_ANG [deg]')
191 ofile = self.config.data.log_dir + '/susi_mod_states_warning_TIMESTSAMPS.png'
192 plt.savefig(ofile)
193 plt.close()
194 # tmestamp difference
195 fig = plt.figure(figsize=(15, 7))
196 fig.suptitle('')
197 plt.plot(np.diff(self.tstmps), '*')
198 plt.xlabel('Frame')
199 plt.xlabel('Timestamp difference [us]')
200 ofile = self.config.data.log_dir + '/susi_mod_states_warning_TIMESTSAMPS_DIFF.png'
201 plt.savefig(ofile)
202 plt.close()
203 # PMU_ANG vs timestamp
204 if len(self.states) > 0:
205 fig = plt.figure(figsize=(15, 7))
206 fig.suptitle('')
207 plt.plot(self.tstmps, self.angles, '*')
208 plt.xlabel(r'Timestamps [$\mu$s]')
209 plt.ylabel('PMU_ANG [deg]')
210 ofile = self.config.data.log_dir + '/susi_mod_states_warning_PMU_ANG.png'
211 plt.savefig(ofile)
212 plt.close()
213 # Identified mod state vs timestamp
214 if len(self.states) > 0 and len(self.states) == len(self.tstmps):
215 fig = plt.figure(figsize=(15, 7))
216 fig.suptitle('')
217 plt.plot(self.tstmps, self.states, '*')
218 plt.xlabel(r'Timestamps [$\mu$s]')
219 plt.ylabel('Detected modulation state')
220 ofile = self.config.data.log_dir + '/susi_mod_states_warning_MOD_STATES.png'
221 plt.savefig(ofile)
222 plt.close()
223 logger.warning('Saved some usefull plots to %s', os.path.dirname(ofile) + '/susi_mod_states_warning_...png')