16080442806

Committed 04 Jul 2025 08:10PM UTC coverage: 75.191% (+0.6%) from 74.575%

Build # 16080442806

Build Type

push

github

Committed by

paulmthompson

Commit Message

remove concrete filters header

Run Details

14051 of 18687 relevant lines covered (75.19%)

1040.01 hits per line

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

89.0

/src/WhiskerToolbox/DataManager/utils/filter/filter.cpp

#include "filter.hpp"
#include "AnalogTimeSeries/Analog_Time_Series.hpp"
#include "new_filter_bridge.hpp"

#include <algorithm>
#include <cmath>
#include <iostream>
#include <map>
#include <numeric>
#include <stdexcept>

// ========== FilterOptions Validation ==========

bool FilterOptions::isValid() const {
    return getValidationError().empty();
}

std::string FilterOptions::getValidationError() const {
    if (order < 1 || order > max_filter_order) {
        return "Filter order must be between 1 and " + std::to_string(max_filter_order);
    }

    if (sampling_rate_hz <= 0.0) {
        return "Sampling rate must be positive";
    }

    if (cutoff_frequency_hz <= 0.0) {
        return "Cutoff frequency must be positive";
    }

    double nyquist = sampling_rate_hz / 2.0;
    if (cutoff_frequency_hz >= nyquist) {
        return "Cutoff frequency must be less than Nyquist frequency (" +
               std::to_string(nyquist) + " Hz)";
    }

    // Additional validation for band filters
    if (response == FilterResponse::BandPass || response == FilterResponse::BandStop) {
        // For RBJ filters, we use cutoff_frequency_hz as center frequency and q_factor
        if (type == FilterType::RBJ) {
            if (cutoff_frequency_hz >= nyquist) {
                return "Center frequency must be less than Nyquist frequency";
            }
            if (q_factor <= 0.0) {
                return "Q factor must be positive for RBJ band filters";
            }
        } else {
            // For IIR filters, we need both cutoff frequencies
            if (high_cutoff_hz <= cutoff_frequency_hz) {
                return "High cutoff frequency must be greater than low cutoff frequency";
            }
            if (high_cutoff_hz >= nyquist) {
                return "High cutoff frequency must be less than Nyquist frequency";
            }
        }
    }

    // Validation for filter-specific parameters
    if (type == FilterType::ChebyshevI && passband_ripple_db <= 0.0) {
        return "Chebyshev I passband ripple must be positive";
    }

    if (type == FilterType::ChebyshevII && stopband_ripple_db <= 0.0) {
        return "Chebyshev II stopband ripple must be positive";
    }

    if (type == FilterType::RBJ && q_factor <= 0.0) {
        return "RBJ Q factor must be positive";
    }

    return "";// Valid
}

// ========== Sampling Rate Estimation ==========

double estimateSamplingRate(
        AnalogTimeSeries const * analog_time_series,
        std::optional<TimeFrameIndex> start_time,
        std::optional<TimeFrameIndex> end_time) {
    if (!analog_time_series) {
        return 0.0;
    }

    size_t num_samples = analog_time_series->getNumSamples();
    if (num_samples < 2) {
        return 0.0;
    }

    // Get time series for analysis
    auto time_indices = analog_time_series->getTimeSeries();

    // Determine analysis range
    size_t start_idx = 0;
    size_t end_idx = time_indices.size();

    if (start_time.has_value()) {
        auto start_data_idx = analog_time_series->findDataArrayIndexGreaterOrEqual(start_time.value());
        if (start_data_idx.has_value()) {
            start_idx = start_data_idx.value().getValue();
        }
    }

    if (end_time.has_value()) {
        auto end_data_idx = analog_time_series->findDataArrayIndexLessOrEqual(end_time.value());
        if (end_data_idx.has_value()) {
            end_idx = std::min(end_data_idx.value().getValue() + 1, time_indices.size());
        }
    }

    if (end_idx <= start_idx + 1) {
        return 0.0;
    }

    // Calculate time differences
    std::vector<double> time_diffs;
    time_diffs.reserve(end_idx - start_idx - 1);

    for (size_t i = start_idx + 1; i < end_idx; ++i) {
        double dt = static_cast<double>(time_indices[i].getValue() - time_indices[i - 1].getValue());
        if (dt > 0) {
            time_diffs.push_back(dt);
        }
    }

    if (time_diffs.empty()) {
        return 0.0;
    }

    // Use median time difference for robustness
    std::sort(time_diffs.begin(), time_diffs.end());
    double median_dt = time_diffs[time_diffs.size() / 2];

    // Assume time indices are in units that give sampling rate of 1/dt
    // This is a heuristic - users should specify sampling rate explicitly
    return 1.0 / median_dt;
}

// ========== Main Filtering Functions ==========

FilterResult filterAnalogTimeSeries(
        AnalogTimeSeries const * analog_time_series,
        TimeFrameIndex start_time,
        TimeFrameIndex end_time,
        FilterOptions const & options) {
    return filterAnalogTimeSeriesNew(analog_time_series, start_time, end_time, options);
}

FilterResult filterAnalogTimeSeries(
        AnalogTimeSeries const * analog_time_series,
        FilterOptions const & options) {
    return filterAnalogTimeSeriesNew(analog_time_series, options);
}

// ========== Filter Defaults ==========

namespace FilterDefaults {

FilterOptions lowpass(double cutoff_hz, double sampling_rate_hz, int order) {
    FilterOptions options;
    options.type = FilterType::Butterworth;
    options.response = FilterResponse::LowPass;
    options.cutoff_frequency_hz = cutoff_hz;
    options.sampling_rate_hz = sampling_rate_hz;
    options.order = order;
    return options;
}

FilterOptions highpass(double cutoff_hz, double sampling_rate_hz, int order) {
    FilterOptions options;
    options.type = FilterType::Butterworth;
    options.response = FilterResponse::HighPass;
    options.cutoff_frequency_hz = cutoff_hz;
    options.sampling_rate_hz = sampling_rate_hz;
    options.order = order;
    return options;
}

FilterOptions bandpass(double low_cutoff_hz, double high_cutoff_hz, double sampling_rate_hz, int order) {
    FilterOptions options;
    options.type = FilterType::Butterworth;
    options.response = FilterResponse::BandPass;
    options.cutoff_frequency_hz = low_cutoff_hz;
    options.high_cutoff_hz = high_cutoff_hz;
    options.sampling_rate_hz = sampling_rate_hz;
    options.order = order;
    return options;
}

FilterOptions notch(double center_freq_hz, double sampling_rate_hz, double q_factor) {
    FilterOptions options;
    options.type = FilterType::RBJ;
    options.response = FilterResponse::BandStop;
    options.cutoff_frequency_hz = center_freq_hz;
    options.high_cutoff_hz = center_freq_hz;  // Set equal to center frequency to force Q-factor path
    options.sampling_rate_hz = sampling_rate_hz;
    options.q_factor = q_factor;
    options.order = 2;  // RBJ filters are always 2nd order
    return options;
}

}// namespace FilterDefaults

1	#include "filter.hpp"
2	#include "AnalogTimeSeries/Analog_Time_Series.hpp"
3	#include "new_filter_bridge.hpp"
4
5	#include <algorithm>
6	#include <cmath>
7	#include <iostream>
8	#include <map>
9	#include <numeric>
10	#include <stdexcept>
11
12	// ========== FilterOptions Validation ==========
13
14	bool FilterOptions::isValid() const {	93✔
15	return getValidationError().empty();	93✔
16	}
17
18	std::string FilterOptions::getValidationError() const {	106✔
19	if (order < 1 \|\| order > max_filter_order) {	106✔
20	return "Filter order must be between 1 and " + std::to_string(max_filter_order);	12✔
21	}
22
23	if (sampling_rate_hz <= 0.0) {	100✔
24	return "Sampling rate must be positive";	18✔
25	}
26
27	if (cutoff_frequency_hz <= 0.0) {	94✔
28	return "Cutoff frequency must be positive";	6✔
29	}
30
31	double nyquist = sampling_rate_hz / 2.0;	92✔
32	if (cutoff_frequency_hz >= nyquist) {	92✔
33	return "Cutoff frequency must be less than Nyquist frequency (" +	8✔
34	std::to_string(nyquist) + " Hz)";	8✔
35	}
36
37	// Additional validation for band filters
38	if (response == FilterResponse::BandPass \|\| response == FilterResponse::BandStop) {	90✔
39	// For RBJ filters, we use cutoff_frequency_hz as center frequency and q_factor
40	if (type == FilterType::RBJ) {	12✔
41	if (cutoff_frequency_hz >= nyquist) {	3✔
42	return "Center frequency must be less than Nyquist frequency";	×
43	}
44	if (q_factor <= 0.0) {	3✔
45	return "Q factor must be positive for RBJ band filters";	×
46	}
47	} else {
48	// For IIR filters, we need both cutoff frequencies
49	if (high_cutoff_hz <= cutoff_frequency_hz) {	9✔
50	return "High cutoff frequency must be greater than low cutoff frequency";	6✔
51	}
52	if (high_cutoff_hz >= nyquist) {	7✔
53	return "High cutoff frequency must be less than Nyquist frequency";	×
54	}
55	}
56	}
57
58	// Validation for filter-specific parameters
59	if (type == FilterType::ChebyshevI && passband_ripple_db <= 0.0) {	88✔
60	return "Chebyshev I passband ripple must be positive";	6✔
61	}
62
63	if (type == FilterType::ChebyshevII && stopband_ripple_db <= 0.0) {	86✔
64	return "Chebyshev II stopband ripple must be positive";	6✔
65	}
66
67	if (type == FilterType::RBJ && q_factor <= 0.0) {	84✔
68	return "RBJ Q factor must be positive";	6✔
69	}
70
71	return "";// Valid	246✔
72	}
73
74	// ========== Sampling Rate Estimation ==========
75
76	double estimateSamplingRate(	4✔
77	AnalogTimeSeries const * analog_time_series,
78	std::optional<TimeFrameIndex> start_time,
79	std::optional<TimeFrameIndex> end_time) {
80	if (!analog_time_series) {	4✔
81	return 0.0;	1✔
82	}
83
84	size_t num_samples = analog_time_series->getNumSamples();	3✔
85	if (num_samples < 2) {	3✔
86	return 0.0;	2✔
87	}
88
89	// Get time series for analysis
90	auto time_indices = analog_time_series->getTimeSeries();	1✔
91
92	// Determine analysis range
93	size_t start_idx = 0;	1✔
94	size_t end_idx = time_indices.size();	1✔
95
96	if (start_time.has_value()) {	1✔
97	auto start_data_idx = analog_time_series->findDataArrayIndexGreaterOrEqual(start_time.value());	×
98	if (start_data_idx.has_value()) {	×
99	start_idx = start_data_idx.value().getValue();	×
100	}
101	}
102
103	if (end_time.has_value()) {	1✔
104	auto end_data_idx = analog_time_series->findDataArrayIndexLessOrEqual(end_time.value());	×
105	if (end_data_idx.has_value()) {	×
106	end_idx = std::min(end_data_idx.value().getValue() + 1, time_indices.size());	×
107	}
108	}
109
110	if (end_idx <= start_idx + 1) {	1✔
111	return 0.0;	×
112	}
113
114	// Calculate time differences
115	std::vector<double> time_diffs;	1✔
116	time_diffs.reserve(end_idx - start_idx - 1);	1✔
117
118	for (size_t i = start_idx + 1; i < end_idx; ++i) {	100✔
119	double dt = static_cast<double>(time_indices[i].getValue() - time_indices[i - 1].getValue());	99✔
120	if (dt > 0) {	99✔
121	time_diffs.push_back(dt);	99✔
122	}
123	}
124
125	if (time_diffs.empty()) {	1✔
126	return 0.0;	×
127	}
128
129	// Use median time difference for robustness
130	std::sort(time_diffs.begin(), time_diffs.end());	1✔
131	double median_dt = time_diffs[time_diffs.size() / 2];	1✔
132
133	// Assume time indices are in units that give sampling rate of 1/dt
134	// This is a heuristic - users should specify sampling rate explicitly
135	return 1.0 / median_dt;	1✔
136	}	1✔
137
138	// ========== Main Filtering Functions ==========
139
140	FilterResult filterAnalogTimeSeries(	1✔
141	AnalogTimeSeries const * analog_time_series,
142	TimeFrameIndex start_time,
143	TimeFrameIndex end_time,
144	FilterOptions const & options) {
145	return filterAnalogTimeSeriesNew(analog_time_series, start_time, end_time, options);	1✔
146	}
147
148	FilterResult filterAnalogTimeSeries(	24✔
149	AnalogTimeSeries const * analog_time_series,
150	FilterOptions const & options) {
151	return filterAnalogTimeSeriesNew(analog_time_series, options);	24✔
152	}
153
154	// ========== Filter Defaults ==========
155
156	namespace FilterDefaults {
157
158	FilterOptions lowpass(double cutoff_hz, double sampling_rate_hz, int order) {	12✔
159	FilterOptions options;	12✔
160	options.type = FilterType::Butterworth;	12✔
161	options.response = FilterResponse::LowPass;	12✔
162	options.cutoff_frequency_hz = cutoff_hz;	12✔
163	options.sampling_rate_hz = sampling_rate_hz;	12✔
164	options.order = order;	12✔
165	return options;	12✔
166	}
167
168	FilterOptions highpass(double cutoff_hz, double sampling_rate_hz, int order) {	2✔
169	FilterOptions options;	2✔
170	options.type = FilterType::Butterworth;	2✔
171	options.response = FilterResponse::HighPass;	2✔
172	options.cutoff_frequency_hz = cutoff_hz;	2✔
173	options.sampling_rate_hz = sampling_rate_hz;	2✔
174	options.order = order;	2✔
175	return options;	2✔
176	}
177
178	FilterOptions bandpass(double low_cutoff_hz, double high_cutoff_hz, double sampling_rate_hz, int order) {	2✔
179	FilterOptions options;	2✔
180	options.type = FilterType::Butterworth;	2✔
181	options.response = FilterResponse::BandPass;	2✔
182	options.cutoff_frequency_hz = low_cutoff_hz;	2✔
183	options.high_cutoff_hz = high_cutoff_hz;	2✔
184	options.sampling_rate_hz = sampling_rate_hz;	2✔
185	options.order = order;	2✔
186	return options;	2✔
187	}
188
189	FilterOptions notch(double center_freq_hz, double sampling_rate_hz, double q_factor) {	2✔
190	FilterOptions options;	2✔
191	options.type = FilterType::RBJ;	2✔
192	options.response = FilterResponse::BandStop;	2✔
193	options.cutoff_frequency_hz = center_freq_hz;	2✔
194	options.high_cutoff_hz = center_freq_hz; // Set equal to center frequency to force Q-factor path	2✔
195	options.sampling_rate_hz = sampling_rate_hz;	2✔
196	options.q_factor = q_factor;	2✔
197	options.order = 2; // RBJ filters are always 2nd order	2✔
198	return options;	2✔
199	}
200
201	}// namespace FilterDefaults

paulmthompson / WhiskerToolbox / 16080442806

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