17270491352

Committed 27 Aug 2025 02:57PM UTC coverage: 65.333%. Remained the same

Build # 17270491352

Build Type

push

github

Committed by

paulmthompson

Commit Message

Merge branch 'main' of https://github.com/paulmthompson/WhiskerToolbox

Run Details

352 of 628 new or added lines in 92 files covered. (56.05%)

357 existing lines in 24 files now uncovered.

26429 of 40453 relevant lines covered (65.33%)

1119.34 hits per line

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

80.33

/src/DataManager/transforms/AnalogTimeSeries/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;
}

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

paulmthompson / WhiskerToolbox / 17270491352

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