17335530769

Committed 29 Aug 2025 09:57PM UTC coverage: 66.478% (+0.3%) from 66.194%

Build # 17335530769

Build Type

push

github

Committed by

paulmthompson

Commit Message

update testing for analog interval threshold

Run Details

113 of 116 new or added lines in 3 files covered. (97.41%)

103 existing lines in 6 files now uncovered.

27064 of 40711 relevant lines covered (66.48%)

1114.57 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

93.51

/src/DataManager/transforms/AnalogTimeSeries/AnalogHilbertPhase/analog_hilbert_phase.cpp

#include "analog_hilbert_phase.hpp"

#include "AnalogTimeSeries/Analog_Time_Series.hpp"
#include "utils/armadillo_wrap/analog_armadillo.hpp"

#include <armadillo>

#include <algorithm>
#include <complex>
#include <iostream>
#include <numbers>
#include <numeric>//std::iota
#include <vector>

/**
     * @brief Represents a continuous chunk of data in the time series
     */
struct DataChunk {
    DataArrayIndex start_idx;         // Start index in original timestamps
    DataArrayIndex end_idx;           // End index in original timestamps (exclusive)
    TimeFrameIndex output_start;      // Start position in output array
    TimeFrameIndex output_end;        // End position in output array (exclusive)
    std::vector<float> values;        // Values for this chunk
    std::vector<TimeFrameIndex> times;// Timestamps for this chunk
};

/**
     * @brief Detects discontinuities in the time series and splits into continuous chunks
     * @param analog_time_series Input time series
     * @param threshold Maximum gap size before considering it a discontinuity
     * @return Vector of continuous data chunks
     */
std::vector<DataChunk> detectChunks(AnalogTimeSeries const * analog_time_series, size_t threshold) {
    std::vector<DataChunk> chunks;

    auto const & timestamps = analog_time_series->getTimeSeries();
    auto const & values = analog_time_series->getAnalogTimeSeries();

    if (timestamps.empty()) {
        return chunks;
    }

    DataArrayIndex chunk_start = DataArrayIndex(0);
    TimeFrameIndex last_time = timestamps[0];

    for (DataArrayIndex i = DataArrayIndex(1); i < DataArrayIndex(timestamps.size()); ++i) {
        TimeFrameIndex current_time = timestamps[i.getValue()];
        TimeFrameIndex gap = current_time - last_time;

        // If gap is larger than threshold, end current chunk and start new one
        if (gap.getValue() > static_cast<int64_t>(threshold)) {

            std::vector<float> chunk_values;
            std::vector<TimeFrameIndex> chunk_times;
            chunk_values.reserve(i - chunk_start);
            chunk_times.reserve(i - chunk_start);
            for (DataArrayIndex j = chunk_start; j < i; ++j) {
                chunk_values.push_back(values[j.getValue()]);
                chunk_times.push_back(timestamps[j.getValue()]);
            }

            // End chunk at last valid time + 1 (exclusive)
            DataChunk chunk{.start_idx = chunk_start,
                            .end_idx = i,
                            .output_start = timestamps[chunk_start.getValue()],
                            .output_end = last_time + TimeFrameIndex(1),
                            .values = std::move(chunk_values),
                            .times = std::move(chunk_times)};

            chunks.push_back(std::move(chunk));
            chunk_start = i;
        }
        last_time = current_time;
    }

    // Add final chunk
    std::vector<float> final_chunk_values;
    std::vector<TimeFrameIndex> final_chunk_times;
    final_chunk_values.reserve(timestamps.size() - chunk_start.getValue());
    final_chunk_times.reserve(timestamps.size() - chunk_start.getValue());
    for (DataArrayIndex j = chunk_start; j < DataArrayIndex(timestamps.size()); ++j) {
        final_chunk_values.push_back(values[j.getValue()]);
        final_chunk_times.push_back(timestamps[j.getValue()]);
    }

    DataChunk final_chunk{.start_idx = chunk_start,
                          .end_idx = DataArrayIndex(timestamps.size()),
                          .output_start = timestamps[chunk_start.getValue()],
                          .output_end = timestamps.back() + TimeFrameIndex(1),
                          .values = std::move(final_chunk_values),
                          .times = std::move(final_chunk_times)};

    chunks.push_back(std::move(final_chunk));

    return chunks;
}

/**
     * @brief Processes a single continuous chunk using Hilbert transform
     * @param chunk The data chunk to process
     * @param phaseParams Parameters for the calculation
     * @return Vector of phase values for this chunk
     */
std::vector<float> processChunk(DataChunk const & chunk, HilbertPhaseParams const & phaseParams) {
    if (chunk.values.empty()) {
        return {};
    }

    std::cout << "Processing chunk with " << chunk.values.size() << " values" << std::endl;

    // First check for NaN values and remove them
    std::vector<float> clean_values;
    std::vector<TimeFrameIndex> clean_times;
    clean_values.reserve(chunk.values.size());
    clean_times.reserve(chunk.times.size());

    for (size_t i = 0; i < chunk.values.size(); ++i) {
        if (!std::isnan(chunk.values[i])) {
            clean_values.push_back(chunk.values[i]);
            clean_times.push_back(chunk.times[i]);
        }
    }

    // If all values were NaN, return empty vector
    if (clean_values.empty()) {
        return std::vector<float>(static_cast<size_t>(chunk.output_end.getValue() - chunk.output_start.getValue()), 0.0f);
    }

    // Convert to arma::vec for processing
    arma::vec signal(clean_values.size());
    std::copy(clean_values.begin(), clean_values.end(), signal.begin());

    // Calculate sampling rate from timestamps
    double dt = 1.0 / 1000.0;// Default to 1kHz if we can't calculate
    if (clean_times.size() > 1) {
        // Find minimum time difference between consecutive samples
        double min_dt = std::numeric_limits<double>::max();
        for (size_t i = 1; i < clean_times.size(); ++i) {
            double diff = static_cast<double>(clean_times[i].getValue() - clean_times[i - 1].getValue());
            if (diff > 0) {
                min_dt = std::min(min_dt, diff);
            }
        }
        if (min_dt != std::numeric_limits<double>::max()) {
            dt = min_dt / 1000.0;// Convert from TimeFrameIndex units to seconds
        }
    }
    double const Fs = 1.0 / dt;

    // Validate frequency parameters (but continue processing regardless)
    double const nyquist = Fs / 2.0;
    if (phaseParams.lowFrequency <= 0 || phaseParams.highFrequency <= 0 ||
        phaseParams.lowFrequency >= nyquist || phaseParams.highFrequency >= nyquist ||
        phaseParams.lowFrequency >= phaseParams.highFrequency) {
        // Log warning but continue processing
        std::cerr << "hilbert_phase: Invalid frequency parameters for chunk. "
                  << "Low: " << phaseParams.lowFrequency
                  << ", High: " << phaseParams.highFrequency
                  << ", Nyquist: " << nyquist << std::endl;
    }

    // Perform FFT
    arma::cx_vec X = arma::fft(signal).eval();

    // Create analytic signal by zeroing negative frequencies
    arma::uword const n = X.n_elem;
    arma::uword const halfway = (n + 1) / 2;

    // Zero out negative frequencies
    for (arma::uword i = halfway; i < n; ++i) {
        X(i) = std::complex<double>(0.0, 0.0);
    }

    // Double the positive frequencies (except DC and Nyquist if present)
    for (arma::uword i = 1; i < halfway; ++i) {
        X(i) *= 2.0;
    }

    // Compute inverse FFT to get the analytic signal
    arma::cx_vec const analytic_signal = arma::ifft(X).eval();

    // Extract phase using vectorized operation
    arma::vec const phase = arma::arg(analytic_signal);

    // Convert back to standard vector
    std::vector<float> phase_values = arma::conv_to<std::vector<float>>::from(phase);

    // Create output vector for this chunk only
    auto chunk_size = static_cast<size_t>(chunk.output_end.getValue() - chunk.output_start.getValue());
    std::vector<float> output_phase(chunk_size, 0.0f);

    // Fill in the original points (excluding NaN values)
    for (size_t i = 0; i < clean_times.size(); ++i) {
        auto output_idx = static_cast<size_t>(clean_times[i].getValue() - chunk.output_start.getValue());
        if (output_idx < output_phase.size()) {
            output_phase[output_idx] = phase_values[i];
        }
    }

    // Interpolate small gaps within this chunk only
    for (size_t i = 1; i < clean_times.size(); ++i) {
        int64_t gap = clean_times[i].getValue() - clean_times[i - 1].getValue();
        if (gap > 1 && static_cast<size_t>(gap) <= phaseParams.discontinuityThreshold) {
            // Linear interpolation for points in between
            float phase_start = phase_values[i - 1];
            float phase_end = phase_values[i];

            // Handle phase wrapping
            if (phase_end - phase_start > static_cast<float>(std::numbers::pi)) {
                phase_start += 2.0f * static_cast<float>(std::numbers::pi);
            } else if (phase_start - phase_end > static_cast<float>(std::numbers::pi)) {
                phase_end += 2.0f * static_cast<float>(std::numbers::pi);
            }

            for (int64_t j = 1; j < gap; ++j) {
                float t = static_cast<float>(j) / static_cast<float>(gap);
                float interpolated_phase = phase_start + t * (phase_end - phase_start);

                // Wrap back to [-π, π]
                while (interpolated_phase > static_cast<float>(std::numbers::pi)) interpolated_phase -= 2.0f * static_cast<float>(std::numbers::pi);
                while (interpolated_phase <= -static_cast<float>(std::numbers::pi)) interpolated_phase += 2.0f * static_cast<float>(std::numbers::pi);

                auto output_idx = static_cast<size_t>((clean_times[i - 1].getValue() + j) - chunk.output_start.getValue());
                if (output_idx < output_phase.size()) {
                    output_phase[output_idx] = interpolated_phase;
                }
            }
        }
    }

    return output_phase;
}

///////////////////////////////////////////////////////////////////////////////

std::shared_ptr<AnalogTimeSeries> hilbert_phase(
        AnalogTimeSeries const * analog_time_series,
        HilbertPhaseParams const & phaseParams) {
    return hilbert_phase(analog_time_series, phaseParams, [](int) {});
}

std::shared_ptr<AnalogTimeSeries> hilbert_phase(
        AnalogTimeSeries const * analog_time_series,
        HilbertPhaseParams const & phaseParams,
        ProgressCallback progressCallback) {

    // Input validation
    if (!analog_time_series) {
        std::cerr << "hilbert_phase: Input AnalogTimeSeries is null" << std::endl;
        return std::make_shared<AnalogTimeSeries>();
    }

    auto const & timestamps = analog_time_series->getTimeSeries();
    if (timestamps.empty()) {
        std::cerr << "hilbert_phase: Input time series is empty" << std::endl;
        return std::make_shared<AnalogTimeSeries>();
    }

    if (progressCallback) {
        progressCallback(5);
    }

    // Detect discontinuous chunks
    auto chunks = detectChunks(analog_time_series, phaseParams.discontinuityThreshold);
    if (chunks.empty()) {
        std::cerr << "hilbert_phase: No valid chunks detected" << std::endl;
        return std::make_shared<AnalogTimeSeries>();
    }

    // Determine total output size based on last chunk's end
    auto const & last_chunk = chunks.back();
    auto total_size = static_cast<size_t>(last_chunk.output_end.getValue());

    // Create output vectors with proper size
    std::vector<float> output_data(total_size, 0.0f);
    std::vector<TimeFrameIndex> output_times;
    output_times.reserve(total_size);
    for (size_t i = 0; i < total_size; ++i) {
        output_times.push_back(TimeFrameIndex(static_cast<int64_t>(i)));
    }

    // Process each chunk
    size_t total_chunks = chunks.size();
    for (size_t i = 0; i < chunks.size(); ++i) {
        auto const & chunk = chunks[i];

        // Process chunk
        auto chunk_phase = processChunk(chunk, phaseParams);

        // Copy chunk results to output
        if (!chunk_phase.empty()) {
            auto start_idx = static_cast<size_t>(chunk.output_start.getValue());
            size_t end_idx = std::min(start_idx + chunk_phase.size(), output_data.size());
            std::copy(chunk_phase.begin(),
                      chunk_phase.begin() + static_cast<long int>(end_idx - start_idx),
                      output_data.begin() + static_cast<long int>(start_idx));
        }

        // Update progress
        if (progressCallback) {
            int progress = 5 + static_cast<int>((90.0f * static_cast<float>(i + 1)) / static_cast<float>(total_chunks));
            progressCallback(progress);
        }
    }

    // Create result
    auto result = std::make_shared<AnalogTimeSeries>(std::move(output_data), std::move(output_times));

    if (progressCallback) {
        progressCallback(100);
    }

    return result;
}

///////////////////////////////////////////////////////////////////////////////

std::string HilbertPhaseOperation::getName() const {
    return "Hilbert Phase";
}

std::type_index HilbertPhaseOperation::getTargetInputTypeIndex() const {
    return typeid(std::shared_ptr<AnalogTimeSeries>);
}

bool HilbertPhaseOperation::canApply(DataTypeVariant const & dataVariant) const {
    if (!std::holds_alternative<std::shared_ptr<AnalogTimeSeries>>(dataVariant)) {
        return false;
    }

    auto const * ptr_ptr = std::get_if<std::shared_ptr<AnalogTimeSeries>>(&dataVariant);
    return ptr_ptr && *ptr_ptr;
}

std::unique_ptr<TransformParametersBase> HilbertPhaseOperation::getDefaultParameters() const {
    return std::make_unique<HilbertPhaseParams>();
}

DataTypeVariant HilbertPhaseOperation::execute(
        DataTypeVariant const & dataVariant,
        TransformParametersBase const * transformParameters) {
    return execute(dataVariant, transformParameters, nullptr);
}

DataTypeVariant HilbertPhaseOperation::execute(
        DataTypeVariant const & dataVariant,
        TransformParametersBase const * transformParameters,
        ProgressCallback progressCallback) {

    auto const * ptr_ptr = std::get_if<std::shared_ptr<AnalogTimeSeries>>(&dataVariant);

    if (!ptr_ptr || !(*ptr_ptr)) {
        std::cerr << "HilbertPhaseOperation::execute: Invalid input data variant" << std::endl;
        return {};
    }

    AnalogTimeSeries const * analog_raw_ptr = (*ptr_ptr).get();

    HilbertPhaseParams currentParams;

    if (transformParameters != nullptr) {
        auto const * specificParams =
                dynamic_cast<HilbertPhaseParams const *>(transformParameters);

        if (specificParams) {
            currentParams = *specificParams;
        } else {
            std::cerr << "HilbertPhaseOperation::execute: Incompatible parameter type, using defaults" << std::endl;
        }
    }

    std::shared_ptr<AnalogTimeSeries> result = hilbert_phase(
            analog_raw_ptr, currentParams, progressCallback);

    if (!result) {
        std::cerr << "HilbertPhaseOperation::execute: Phase calculation failed" << std::endl;
        return {};
    }

    return result;
}

1	#include "analog_hilbert_phase.hpp"
2
3	#include "AnalogTimeSeries/Analog_Time_Series.hpp"
4	#include "utils/armadillo_wrap/analog_armadillo.hpp"
5
6	#include <armadillo>
7
8	#include <algorithm>
9	#include <complex>
10	#include <iostream>
11	#include <numbers>
12	#include <numeric>//std::iota
13	#include <vector>
14
15	/**
16	* @brief Represents a continuous chunk of data in the time series
17	*/
18	struct DataChunk {
19	DataArrayIndex start_idx; // Start index in original timestamps
20	DataArrayIndex end_idx; // End index in original timestamps (exclusive)
21	TimeFrameIndex output_start; // Start position in output array
22	TimeFrameIndex output_end; // End position in output array (exclusive)
23	std::vector<float> values; // Values for this chunk
24	std::vector<TimeFrameIndex> times;// Timestamps for this chunk
25	};
26
27	/**
28	* @brief Detects discontinuities in the time series and splits into continuous chunks
29	* @param analog_time_series Input time series
30	* @param threshold Maximum gap size before considering it a discontinuity
31	* @return Vector of continuous data chunks
32	*/
33	std::vector<DataChunk> detectChunks(AnalogTimeSeries const * analog_time_series, size_t threshold) {	24✔
34	std::vector<DataChunk> chunks;	24✔
35
36	auto const & timestamps = analog_time_series->getTimeSeries();	24✔
37	auto const & values = analog_time_series->getAnalogTimeSeries();	24✔
38
39	if (timestamps.empty()) {	24✔
40	return chunks;	×
41	}
42
43	DataArrayIndex chunk_start = DataArrayIndex(0);	24✔
44	TimeFrameIndex last_time = timestamps[0];	24✔
45
46	for (DataArrayIndex i = DataArrayIndex(1); i < DataArrayIndex(timestamps.size()); ++i) {	1,707✔
47	TimeFrameIndex current_time = timestamps[i.getValue()];	1,683✔
48	TimeFrameIndex gap = current_time - last_time;	1,683✔
49
50	// If gap is larger than threshold, end current chunk and start new one
51	if (gap.getValue() > static_cast<int64_t>(threshold)) {	1,683✔
52
53	std::vector<float> chunk_values;	10✔
54	std::vector<TimeFrameIndex> chunk_times;	10✔
55	chunk_values.reserve(i - chunk_start);	10✔
56	chunk_times.reserve(i - chunk_start);	10✔
57	for (DataArrayIndex j = chunk_start; j < i; ++j) {	50✔
58	chunk_values.push_back(values[j.getValue()]);	40✔
59	chunk_times.push_back(timestamps[j.getValue()]);	40✔
60	}
61
62	// End chunk at last valid time + 1 (exclusive)
63	DataChunk chunk{.start_idx = chunk_start,	10✔
64	.end_idx = i,
65	.output_start = timestamps[chunk_start.getValue()],	10✔
66	.output_end = last_time + TimeFrameIndex(1),	10✔
67	.values = std::move(chunk_values),	10✔
68	.times = std::move(chunk_times)};	30✔
69
70	chunks.push_back(std::move(chunk));	10✔
71	chunk_start = i;	10✔
72	}	10✔
73	last_time = current_time;	1,683✔
74	}
75
76	// Add final chunk
77	std::vector<float> final_chunk_values;	24✔
78	std::vector<TimeFrameIndex> final_chunk_times;	24✔
79	final_chunk_values.reserve(timestamps.size() - chunk_start.getValue());	24✔
80	final_chunk_times.reserve(timestamps.size() - chunk_start.getValue());	24✔
81	for (DataArrayIndex j = chunk_start; j < DataArrayIndex(timestamps.size()); ++j) {	1,691✔
82	final_chunk_values.push_back(values[j.getValue()]);	1,667✔
83	final_chunk_times.push_back(timestamps[j.getValue()]);	1,667✔
84	}
85
86	DataChunk final_chunk{.start_idx = chunk_start,	24✔
87	.end_idx = DataArrayIndex(timestamps.size()),
88	.output_start = timestamps[chunk_start.getValue()],	24✔
89	.output_end = timestamps.back() + TimeFrameIndex(1),	24✔
90	.values = std::move(final_chunk_values),	24✔
91	.times = std::move(final_chunk_times)};	72✔
92
93	chunks.push_back(std::move(final_chunk));	24✔
94
95	return chunks;	24✔
96	}	24✔
97
98	/**
99	* @brief Processes a single continuous chunk using Hilbert transform
100	* @param chunk The data chunk to process
101	* @param phaseParams Parameters for the calculation
102	* @return Vector of phase values for this chunk
103	*/
104	std::vector<float> processChunk(DataChunk const & chunk, HilbertPhaseParams const & phaseParams) {	34✔
105	if (chunk.values.empty()) {	34✔
106	return {};	×
107	}
108
109	std::cout << "Processing chunk with " << chunk.values.size() << " values" << std::endl;	34✔
110
111	// First check for NaN values and remove them
112	std::vector<float> clean_values;	34✔
113	std::vector<TimeFrameIndex> clean_times;	34✔
114	clean_values.reserve(chunk.values.size());	34✔
115	clean_times.reserve(chunk.times.size());	34✔
116
117	for (size_t i = 0; i < chunk.values.size(); ++i) {	1,741✔
118	if (!std::isnan(chunk.values[i])) {	1,707✔
119	clean_values.push_back(chunk.values[i]);	1,706✔
120	clean_times.push_back(chunk.times[i]);	1,706✔
121	}
122	}
123
124	// If all values were NaN, return empty vector
125	if (clean_values.empty()) {	34✔
126	return std::vector<float>(static_cast<size_t>(chunk.output_end.getValue() - chunk.output_start.getValue()), 0.0f);	×
127	}
128
129	// Convert to arma::vec for processing
130	arma::vec signal(clean_values.size());	34✔
131	std::copy(clean_values.begin(), clean_values.end(), signal.begin());	34✔
132
133	// Calculate sampling rate from timestamps
134	double dt = 1.0 / 1000.0;// Default to 1kHz if we can't calculate	34✔
135	if (clean_times.size() > 1) {	34✔
136	// Find minimum time difference between consecutive samples
137	double min_dt = std::numeric_limits<double>::max();	28✔
138	for (size_t i = 1; i < clean_times.size(); ++i) {	1,700✔
139	double diff = static_cast<double>(clean_times[i].getValue() - clean_times[i - 1].getValue());	1,672✔
140	if (diff > 0) {	1,672✔
141	min_dt = std::min(min_dt, diff);	1,672✔
142	}
143	}
144	if (min_dt != std::numeric_limits<double>::max()) {	28✔
145	dt = min_dt / 1000.0;// Convert from TimeFrameIndex units to seconds	28✔
146	}
147	}
148	double const Fs = 1.0 / dt;	34✔
149
150	// Validate frequency parameters (but continue processing regardless)
151	double const nyquist = Fs / 2.0;	34✔
152	if (phaseParams.lowFrequency <= 0 \|\| phaseParams.highFrequency <= 0 \|\|	34✔
153	phaseParams.lowFrequency >= nyquist \|\| phaseParams.highFrequency >= nyquist \|\|	33✔
154	phaseParams.lowFrequency >= phaseParams.highFrequency) {	31✔
155	// Log warning but continue processing
156	std::cerr << "hilbert_phase: Invalid frequency parameters for chunk. "
157	<< "Low: " << phaseParams.lowFrequency	4✔
158	<< ", High: " << phaseParams.highFrequency	4✔
159	<< ", Nyquist: " << nyquist << std::endl;	4✔
160	}
161
162	// Perform FFT
163	arma::cx_vec X = arma::fft(signal).eval();	34✔
164
165	// Create analytic signal by zeroing negative frequencies
166	arma::uword const n = X.n_elem;	34✔
167	arma::uword const halfway = (n + 1) / 2;	34✔
168
169	// Zero out negative frequencies
170	for (arma::uword i = halfway; i < n; ++i) {	880✔
171	X(i) = std::complex<double>(0.0, 0.0);	1,692✔
172	}
173
174	// Double the positive frequencies (except DC and Nyquist if present)
175	for (arma::uword i = 1; i < halfway; ++i) {	860✔
176	X(i) *= 2.0;	826✔
177	}
178
179	// Compute inverse FFT to get the analytic signal
180	arma::cx_vec const analytic_signal = arma::ifft(X).eval();	34✔
181
182	// Extract phase using vectorized operation
183	arma::vec const phase = arma::arg(analytic_signal);	34✔
184
185	// Convert back to standard vector
186	std::vector<float> phase_values = arma::conv_to<std::vector<float>>::from(phase);	34✔
187
188	// Create output vector for this chunk only
189	auto chunk_size = static_cast<size_t>(chunk.output_end.getValue() - chunk.output_start.getValue());	34✔
190	std::vector<float> output_phase(chunk_size, 0.0f);	102✔
191
192	// Fill in the original points (excluding NaN values)
193	for (size_t i = 0; i < clean_times.size(); ++i) {	1,740✔
194	auto output_idx = static_cast<size_t>(clean_times[i].getValue() - chunk.output_start.getValue());	1,706✔
195	if (output_idx < output_phase.size()) {	1,706✔
196	output_phase[output_idx] = phase_values[i];	1,706✔
197	}
198	}
199
200	// Interpolate small gaps within this chunk only
201	for (size_t i = 1; i < clean_times.size(); ++i) {	1,706✔
202	int64_t gap = clean_times[i].getValue() - clean_times[i - 1].getValue();	1,672✔
203	if (gap > 1 && static_cast<size_t>(gap) <= phaseParams.discontinuityThreshold) {	1,672✔
204	// Linear interpolation for points in between
205	float phase_start = phase_values[i - 1];	42✔
206	float phase_end = phase_values[i];	42✔
207
208	// Handle phase wrapping
209	if (phase_end - phase_start > static_cast<float>(std::numbers::pi)) {	42✔
210	phase_start += 2.0f * static_cast<float>(std::numbers::pi);	×
211	} else if (phase_start - phase_end > static_cast<float>(std::numbers::pi)) {	42✔
212	phase_end += 2.0f * static_cast<float>(std::numbers::pi);	7✔
213	}
214
215	for (int64_t j = 1; j < gap; ++j) {	874✔
216	float t = static_cast<float>(j) / static_cast<float>(gap);	832✔
217	float interpolated_phase = phase_start + t * (phase_end - phase_start);	832✔
218
219	// Wrap back to [-π, π]
220	while (interpolated_phase > static_cast<float>(std::numbers::pi)) interpolated_phase -= 2.0f * static_cast<float>(std::numbers::pi);	1,047✔
221	while (interpolated_phase <= -static_cast<float>(std::numbers::pi)) interpolated_phase += 2.0f * static_cast<float>(std::numbers::pi);	832✔
222
223	auto output_idx = static_cast<size_t>((clean_times[i - 1].getValue() + j) - chunk.output_start.getValue());	832✔
224	if (output_idx < output_phase.size()) {	832✔
225	output_phase[output_idx] = interpolated_phase;	832✔
226	}
227	}
228	}
229	}
230
231	return output_phase;	34✔
232	}	34✔
233
234	///////////////////////////////////////////////////////////////////////////////
235
236	std::shared_ptr<AnalogTimeSeries> hilbert_phase(	19✔
237	AnalogTimeSeries const * analog_time_series,
238	HilbertPhaseParams const & phaseParams) {
239	return hilbert_phase(analog_time_series, phaseParams, [](int) {});	19✔
240	}
241
242	std::shared_ptr<AnalogTimeSeries> hilbert_phase(	27✔
243	AnalogTimeSeries const * analog_time_series,
244	HilbertPhaseParams const & phaseParams,
245	ProgressCallback progressCallback) {
246
247	// Input validation
248	if (!analog_time_series) {	27✔
249	std::cerr << "hilbert_phase: Input AnalogTimeSeries is null" << std::endl;	2✔
250	return std::make_shared<AnalogTimeSeries>();	2✔
251	}
252
253	auto const & timestamps = analog_time_series->getTimeSeries();	25✔
254	if (timestamps.empty()) {	25✔
255	std::cerr << "hilbert_phase: Input time series is empty" << std::endl;	1✔
256	return std::make_shared<AnalogTimeSeries>();	1✔
257	}
258
259	if (progressCallback) {	24✔
260	progressCallback(5);	23✔
261	}
262
263	// Detect discontinuous chunks
264	auto chunks = detectChunks(analog_time_series, phaseParams.discontinuityThreshold);	24✔
265	if (chunks.empty()) {	24✔
266	std::cerr << "hilbert_phase: No valid chunks detected" << std::endl;	×
267	return std::make_shared<AnalogTimeSeries>();	×
268	}
269
270	// Determine total output size based on last chunk's end
271	auto const & last_chunk = chunks.back();	24✔
272	auto total_size = static_cast<size_t>(last_chunk.output_end.getValue());	24✔
273
274	// Create output vectors with proper size
275	std::vector<float> output_data(total_size, 0.0f);	72✔
276	std::vector<TimeFrameIndex> output_times;	24✔
277	output_times.reserve(total_size);	24✔
278	for (size_t i = 0; i < total_size; ++i) {	8,861✔
279	output_times.push_back(TimeFrameIndex(static_cast<int64_t>(i)));	8,837✔
280	}
281
282	// Process each chunk
283	size_t total_chunks = chunks.size();	24✔
284	for (size_t i = 0; i < chunks.size(); ++i) {	58✔
285	auto const & chunk = chunks[i];	34✔
286
287	// Process chunk
288	auto chunk_phase = processChunk(chunk, phaseParams);	34✔
289
290	// Copy chunk results to output
291	if (!chunk_phase.empty()) {	34✔
292	auto start_idx = static_cast<size_t>(chunk.output_start.getValue());	34✔
293	size_t end_idx = std::min(start_idx + chunk_phase.size(), output_data.size());	34✔
294	std::copy(chunk_phase.begin(),	102✔
295	chunk_phase.begin() + static_cast<long int>(end_idx - start_idx),	68✔
296	output_data.begin() + static_cast<long int>(start_idx));	68✔
297	}
298
299	// Update progress
300	if (progressCallback) {	34✔
301	int progress = 5 + static_cast<int>((90.0f * static_cast<float>(i + 1)) / static_cast<float>(total_chunks));	33✔
302	progressCallback(progress);	33✔
303	}
304	}	34✔
305
306	// Create result
307	auto result = std::make_shared<AnalogTimeSeries>(std::move(output_data), std::move(output_times));	24✔
308
309	if (progressCallback) {	24✔
310	progressCallback(100);	23✔
311	}
312
313	return result;	24✔
314	}	25✔
315
316	///////////////////////////////////////////////////////////////////////////////
317
318	std::string HilbertPhaseOperation::getName() const {	24✔
319	return "Hilbert Phase";	72✔
320	}
321
322	std::type_index HilbertPhaseOperation::getTargetInputTypeIndex() const {	24✔
323	return typeid(std::shared_ptr<AnalogTimeSeries>);	24✔
324	}
325
326	bool HilbertPhaseOperation::canApply(DataTypeVariant const & dataVariant) const {	4✔
327	if (!std::holds_alternative<std::shared_ptr<AnalogTimeSeries>>(dataVariant)) {	4✔
328	return false;	×
329	}
330
331	auto const * ptr_ptr = std::get_if<std::shared_ptr<AnalogTimeSeries>>(&dataVariant);	4✔
332	return ptr_ptr && *ptr_ptr;	4✔
333	}
334
335	std::unique_ptr<TransformParametersBase> HilbertPhaseOperation::getDefaultParameters() const {	4✔
336	return std::make_unique<HilbertPhaseParams>();	4✔
337	}
338
339	DataTypeVariant HilbertPhaseOperation::execute(	1✔
340	DataTypeVariant const & dataVariant,
341	TransformParametersBase const * transformParameters) {
342	return execute(dataVariant, transformParameters, nullptr);	1✔
343	}
344
345	DataTypeVariant HilbertPhaseOperation::execute(	4✔
346	DataTypeVariant const & dataVariant,
347	TransformParametersBase const * transformParameters,
348	ProgressCallback progressCallback) {
349
350	auto const * ptr_ptr = std::get_if<std::shared_ptr<AnalogTimeSeries>>(&dataVariant);	4✔
351
352	if (!ptr_ptr \|\| !(*ptr_ptr)) {	4✔
UNCOV 353	std::cerr << "HilbertPhaseOperation::execute: Invalid input data variant" << std::endl;	×
UNCOV 354	return {};	×
355	}
356
357	AnalogTimeSeries const * analog_raw_ptr = (*ptr_ptr).get();	4✔
358
359	HilbertPhaseParams currentParams;	4✔
360
361	if (transformParameters != nullptr) {	4✔
362	auto const * specificParams =
363	dynamic_cast<HilbertPhaseParams const *>(transformParameters);	4✔
364
365	if (specificParams) {	4✔
366	currentParams = *specificParams;	4✔
367	} else {
UNCOV 368	std::cerr << "HilbertPhaseOperation::execute: Incompatible parameter type, using defaults" << std::endl;	×
369	}
370	}
371
372	std::shared_ptr<AnalogTimeSeries> result = hilbert_phase(	4✔
373	analog_raw_ptr, currentParams, progressCallback);	4✔
374
375	if (!result) {	4✔
UNCOV 376	std::cerr << "HilbertPhaseOperation::execute: Phase calculation failed" << std::endl;	×
UNCOV 377	return {};	×
378	}
379
380	return result;	4✔
381	}	4✔

paulmthompson / WhiskerToolbox / 17335530769

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous