14803779026

Committed 02 May 2025 09:35PM UTC coverage: 15.881% (-0.4%) from 16.234%

Build # 14803779026

Build Type

push

github

Committed by

paulmthompson

Commit Message

working hilbert transform with armadillo

Coverage Stats

1 of 50 new or added lines in 3 files covered. (2.0%)

2 existing lines in 1 file now uncovered.

320 of 2015 relevant lines covered (15.88%)

1.59 hits per line

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

0.0

/src/WhiskerToolbox/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 <complex>

std::shared_ptr<AnalogTimeSeries> hilbert_phase(
        AnalogTimeSeries const * analog_time_series,
        HilbertPhaseParams const & phaseParams) {
    auto phase_series = std::make_shared<AnalogTimeSeries>();

    // Get input data
    auto const & timestamps = analog_time_series->getTimeSeries(); // Copy because we will use these to construct output
    auto const & values = analog_time_series->getAnalogTimeSeries();
    auto output_timestamps = std::vector<size_t>(timestamps.back());
    std::iota(output_timestamps.begin(), output_timestamps.end(), 0);

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

    // Convert to arma::vec for processing
    auto signal = convertAnalogTimeSeriesToMlpackArray(analog_time_series, output_timestamps);
    signal.replace(arma::datum::nan, 0.0);

    // Calculate sampling rate from timestamps
    double dt = 0;
    if (timestamps.size() > 1) {
        dt = (timestamps[timestamps.size()-1] - timestamps[0]) / (timestamps.size() - 1);
    }
    double Fs = 1.0 / dt;

    // Bandpass filter design
    // Convert Hz to normalized frequency (0 to 1, where 1 is Nyquist frequency)
    double low_norm = phaseParams.lowFrequency / (Fs/2.0);
    double high_norm = phaseParams.highFrequency / (Fs/2.0);

    // Check if frequencies are valid
    if (low_norm >= 1.0 || high_norm >= 1.0 || low_norm <= 0 || high_norm <= 0) {
        std::cerr << "Invalid cutoff frequencies for bandpass filter" << std::endl;
        //return phase_series;
    }

    // Use Armadillo's built-in filter functions
    // For simplicity, using a basic butterworth filter design
    arma::vec b, a;

    // Design bandpass filter (4th order Butterworth)
    arma::vec filtered_signal = signal.as_col(); // Placeholder for filtered signal

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

    // Set negative frequencies to zero
    int const halfway = 1 + X.n_elem / 2;
    for (int i = halfway; i < X.n_elem; i++) {
        X(i) = std::complex<double>(0, 0);
    }

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

    // Extract phase
    arma::vec phase(analytic_signal.n_elem);
    for (int i = 0; i < analytic_signal.n_elem; i++) {
        phase(i) = std::arg(analytic_signal(i));
    }

    // Convert back to standard vector and set in the output time series
    auto phase_vec = arma::conv_to<std::vector<float>>::from(phase);
    phase_series->setData(phase_vec, output_timestamps);

    std::cout << "Total size of phase is " << phase_vec.size() << std::endl;
    std::cout << "Max value is " << phase_series->getMaxValue() << std::endl;
    std::cout << "Min value is " << phase_series->getMinValue() << std::endl;

    return phase_series;
}

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 true only if get_if succeeded AND the contained shared_ptr is not null.
    return ptr_ptr && *ptr_ptr;
}

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

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

    if (!ptr_ptr || !(*ptr_ptr)) {
        std::cerr << "HilbertPhaseOperation::execute called with incompatible variant type or null data." << std::endl;
        return {};// Return empty
    }

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

    HilbertPhaseParams currentParams;

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

        if (specificParams) {
            currentParams = *specificParams;
            std::cout << "Using parameters provided by UI." << std::endl;
        } else {
            std::cerr << "Warning: HilbertPhaseOperation received incompatible parameter type (dynamic_cast failed)! Using default parameters." << std::endl;
        }
    } else {
        // No parameter object was provided (nullptr). This might be expected if the
        // operation can run with defaults or was called programmatically.
        std::cout << "HilbertPhaseOperation received null parameters. Using default parameters." << std::endl;
    }

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

    if (!result) {
        std::cerr << "HilbertPhaseOperation::execute: 'calculate_intervals' failed to produce a result." << std::endl;
        return {};// Return empty
    }

    std::cout << "HilbertPhaseOperation executed successfully using variant input." << std::endl;
    return result;
}

1
2	#include "analog_hilbert_phase.hpp"
3
4	#include "AnalogTimeSeries/Analog_Time_Series.hpp"
5	#include "utils/armadillo_wrap/analog_armadillo.hpp"
6
7	#include <armadillo>
8
9	#include <complex>
10
11	std::shared_ptr<AnalogTimeSeries> hilbert_phase(	×
12	AnalogTimeSeries const * analog_time_series,
13	HilbertPhaseParams const & phaseParams) {
14	auto phase_series = std::make_shared<AnalogTimeSeries>();	×
15
16	// Get input data
NEW 17	auto const & timestamps = analog_time_series->getTimeSeries(); // Copy because we will use these to construct output	×
NEW 18	auto const & values = analog_time_series->getAnalogTimeSeries();	×
NEW 19	auto output_timestamps = std::vector<size_t>(timestamps.back());	×
NEW 20	std::iota(output_timestamps.begin(), output_timestamps.end(), 0);	×
21
NEW 22	if (timestamps.empty()) {	×
NEW 23	return phase_series;	×
24	}
25
26	// Convert to arma::vec for processing
NEW 27	auto signal = convertAnalogTimeSeriesToMlpackArray(analog_time_series, output_timestamps);	×
NEW 28	signal.replace(arma::datum::nan, 0.0);	×
29
30	// Calculate sampling rate from timestamps
NEW 31	double dt = 0;	×
NEW 32	if (timestamps.size() > 1) {	×
NEW 33	dt = (timestamps[timestamps.size()-1] - timestamps[0]) / (timestamps.size() - 1);	×
34	}
NEW 35	double Fs = 1.0 / dt;	×
36
37	// Bandpass filter design
38	// Convert Hz to normalized frequency (0 to 1, where 1 is Nyquist frequency)
NEW 39	double low_norm = phaseParams.lowFrequency / (Fs/2.0);	×
NEW 40	double high_norm = phaseParams.highFrequency / (Fs/2.0);	×
41
42	// Check if frequencies are valid
NEW 43	if (low_norm >= 1.0 \|\| high_norm >= 1.0 \|\| low_norm <= 0 \|\| high_norm <= 0) {	×
NEW 44	std::cerr << "Invalid cutoff frequencies for bandpass filter" << std::endl;	×
45	//return phase_series;
46	}
47
48	// Use Armadillo's built-in filter functions
49	// For simplicity, using a basic butterworth filter design
NEW 50	arma::vec b, a;	×
51
52	// Design bandpass filter (4th order Butterworth)
NEW 53	arma::vec filtered_signal = signal.as_col(); // Placeholder for filtered signal	×
54
55	// Perform FFT
NEW 56	auto X = arma::fft(signal).eval();	×
57
58	// Set negative frequencies to zero
NEW 59	int const halfway = 1 + X.n_elem / 2;	×
NEW 60	for (int i = halfway; i < X.n_elem; i++) {	×
NEW 61	X(i) = std::complex<double>(0, 0);	×
62	}
63
64	// Compute inverse FFT to get the analytic signal
NEW 65	auto analytic_signal = arma::ifft(X).eval();	×
66
67	// Extract phase
NEW 68	arma::vec phase(analytic_signal.n_elem);	×
NEW 69	for (int i = 0; i < analytic_signal.n_elem; i++) {	×
NEW 70	phase(i) = std::arg(analytic_signal(i));	×
71	}
72
73	// Convert back to standard vector and set in the output time series
NEW 74	auto phase_vec = arma::conv_to<std::vector<float>>::from(phase);	×
NEW 75	phase_series->setData(phase_vec, output_timestamps);	×
76
NEW 77	std::cout << "Total size of phase is " << phase_vec.size() << std::endl;	×
NEW 78	std::cout << "Max value is " << phase_series->getMaxValue() << std::endl;	×
NEW 79	std::cout << "Min value is " << phase_series->getMinValue() << std::endl;	×
80
UNCOV 81	return phase_series;	×
UNCOV 82	}	×
83
84	std::string HilbertPhaseOperation::getName() const {	×
85	return "Hilbert Phase";	×
86	}
87
88	std::type_index HilbertPhaseOperation::getTargetInputTypeIndex() const {	×
89	return typeid(std::shared_ptr<AnalogTimeSeries>);	×
90	}
91
92	bool HilbertPhaseOperation::canApply(DataTypeVariant const & dataVariant) const {	×
93	if (!std::holds_alternative<std::shared_ptr<AnalogTimeSeries>>(dataVariant)) {	×
94	return false;	×
95	}
96
97	auto const * ptr_ptr = std::get_if<std::shared_ptr<AnalogTimeSeries>>(&dataVariant);	×
98
99	// Return true only if get_if succeeded AND the contained shared_ptr is not null.
100	return ptr_ptr && *ptr_ptr;	×
101	}
102
103	DataTypeVariant HilbertPhaseOperation::execute(DataTypeVariant const & dataVariant, TransformParametersBase const * transformParameters) {	×
104
105	auto const * ptr_ptr = std::get_if<std::shared_ptr<AnalogTimeSeries>>(&dataVariant);	×
106
107	if (!ptr_ptr \|\| !(*ptr_ptr)) {	×
108	std::cerr << "HilbertPhaseOperation::execute called with incompatible variant type or null data." << std::endl;	×
109	return {};// Return empty	×
110	}
111
112	AnalogTimeSeries const * analog_raw_ptr = (*ptr_ptr).get();	×
113
114	HilbertPhaseParams currentParams;	×
115
116	if (transformParameters != nullptr) {	×
117	HilbertPhaseParams const * specificParams = dynamic_cast<HilbertPhaseParams const *>(transformParameters);	×
118
119	if (specificParams) {	×
120	currentParams = *specificParams;	×
121	std::cout << "Using parameters provided by UI." << std::endl;	×
122	} else {
NEW 123	std::cerr << "Warning: HilbertPhaseOperation received incompatible parameter type (dynamic_cast failed)! Using default parameters." << std::endl;	×
124	}
125	} else {
126	// No parameter object was provided (nullptr). This might be expected if the
127	// operation can run with defaults or was called programmatically.
128	std::cout << "HilbertPhaseOperation received null parameters. Using default parameters." << std::endl;	×
129	}
130
131	std::shared_ptr<AnalogTimeSeries> result = hilbert_phase(analog_raw_ptr,	×
132	currentParams);	×
133
134	if (!result) {	×
135	std::cerr << "HilbertPhaseOperation::execute: 'calculate_intervals' failed to produce a result." << std::endl;	×
136	return {};// Return empty	×
137	}
138
139	std::cout << "HilbertPhaseOperation executed successfully using variant input." << std::endl;	×
140	return result;	×
141	}	×

paulmthompson / WhiskerToolbox / 14803779026

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