In the previous post, we discussed preprocessing data with OpenSim. An important part of preprocessing is low-pass filtering of various pieces of data. In this post, we will explore many of the unanswered questions and ambiguities of low-pass filtering in OpenSim. Following stories and personal experience of strange and erroneous results coming from OpenSim filtering, I began to doubt the implementation in OpenSim. To verify the results, I set out to validate the filtering against 2 existing implementations in Python and Matlab.

OpenSim Low-pass Delirium

To begin, there are a multitude of code paths which perform slightly different versions of the same underlying low-pass filter in OpenSim. All the low-pass filtering code eventually arrives at this function:

/**
 * 3rd ORDER LOWPASS IIR BUTTERWORTH DIGITAL FILTER
 *
 * It is assumed that enough memory is allocated at sigf.
 * Note also that the first and last three data points are not filtered.
 *
 *  @param T Sample interval in seconds.
 *  @param fc Cutoff frequency in Hz.
 *  @param N Number of data points in the signal.
 *  @param sig The sampled signal.
 *  @param sigf The filtered signal.
 *
 * @return 0 on success, and -1 on failure.
 */
int Signal::
LowpassIIR(double T,double fc,int N,const double *sig,double *sigf)

Looking at this code reveals that low-pass filtering in OpenSim will always be a 3rd order Butterworth filter.

Moving forward, this function can be called from many diverging paths. The first path is used in some tools, such as the InverseDynamicsTool.

Path 1 - The Storage Class

The first way to reach the underlying low-pass function comes from the Storage::lowpassIIR(double aCutoffFrequency) method which confusingly has the exact same name:

void Storage::
lowpassIIR(double aCutoffFrequency)
{
    int size = getSize();
    double dtmin = getMinTimeStep();
    double avgDt = (_storage[size-1].getTime() - _storage[0].getTime()) / (size-1);

    if(dtmin<SimTK::Eps) {
        log_error("Storage.lowpassIIR: storage cannot be resampled.");
        return;
    }

    // RESAMPLE if the sampling interval is not uniform
    if ((avgDt - dtmin) > SimTK::Eps) {
        dtmin = resample(dtmin, 5);
        size = getSize();
    }

    if(size<(4)) {
        log_error("Storage.lowpassIIR: too few data points to filter.");
        return;
    }

    // LOOP OVER COLUMNS
    int nc = getSmallestNumberOfStates();
    double *signal=NULL;
    Array<double> filt(0.0,size);
    for(int i=0;i<nc;i++) {
        getDataColumn(i,signal);
        Signal::LowpassIIR(dtmin,aCutoffFrequency,size,signal,&filt[0]);
        setDataColumn(i,filt);
    }

    // CLEANUP
    delete[] signal;
}

The Storage::lowpassIIR method is used in the InverseDynamicsTool and AnalyzeTool for low-pass filtering. Looking at the implementation we find:

            if(_lowpassCutoffFrequency>=0) {
                log_info("Low-pass filtering coordinates data with a cutoff "
                         "frequency of {}...", _lowpassCutoffFrequency);
                _coordinateValues->pad(_coordinateValues->getSize()/2);
                _coordinateValues->lowpassIIR(_lowpassCutoffFrequency);
            }

Here the data is being first padded and then resampled if the sampling rate is not sufficient. This is the exact opposite behavior of the second code path, where the data is first resampled and then padded if necessary.

This means that setting the lowpassCutoffFrequency value to a negative number will disable it. Setting it to a positive number will enable the 3rd order Butterworth low-pass filter at that frequency, but the resulting data will NOT be re-trimmed to original length. Additionally, there is currently no way to output the intermediary results of the filtered coordinates. This would be very useful for plotting and validation that the filtering is working properly.

Path 2 - The TableUtilities Utility Method

The second path comes from the TableUtilities::filterLowpass(TimeSeriesTable& table, double cutoffFreq, bool padData) method:

void TableUtilities::filterLowpass(
        TimeSeriesTable& table, double cutoffFreq, bool padData) {
    OPENSIM_THROW_IF(cutoffFreq < 0, Exception,
            "Cutoff frequency must be non-negative; got {}.", cutoffFreq);

    const auto& time = table.getIndependentColumn();
    const auto uniform = isUniform(time);
    const auto& uniformlySampled = uniform.first;
    const auto& dtMin = uniform.second;

    OPENSIM_THROW_IF(
            dtMin < SimTK::Eps, Exception, "Storage cannot be resampled.");

    // Resample if the sampling interval is not uniform.
    if (!uniformlySampled) {
        log_warn("Table not uniformly sampled! Resampling with interval: {}", dtMin);
        table = resampleWithInterval(table, dtMin);
    }

    // Apply padding after resampling, so we preserve the original initial time.
    if (padData) { pad(table, (int)table.getNumRows() / 2); }
    // Uniform sampling should be preserved after padding.
    OPENSIM_THROW_IF(!isUniform(table.getIndependentColumn()).first, Exception,
            "Internal error: non-uniform sampling after padding.");

    const int numRows = (int)table.getNumRows();
    OPENSIM_THROW_IF(numRows < 4, Exception,
            "Expected at least 4 rows to filter, but got {} rows.", numRows);

    SimTK::Vector filtered(numRows);
    for (int icol = 0; icol < (int)table.getNumColumns(); ++icol) {
        SimTK::VectorView column = table.getDependentColumnAtIndex(icol);
        Signal::LowpassIIR(dtMin, cutoffFreq, numRows,
                column.getContiguousScalarData(),
                filtered.updContiguousScalarData());
        table.updDependentColumnAtIndex(icol) = filtered;
    }
}

I am now wondering why this method exists in the first place, since it is nearly the same as the above Storage::lowpassIIR method. Perhaps it is because TimeSeriesTable and Storage do not have compatible APIs, which also does not make sense. Perhaps you have also noticed the mysterious 3rd boolean argument to this function which will come into play later.

Path 3 - The TabOpLowPassFilter Class

The third path comes from the TabOpLowPassFilter class:

/// Apply a low-pass filter to the trajectory.
class OSIMSIMULATION_API TabOpLowPassFilter : public TableOperator {
    OpenSim_DECLARE_CONCRETE_OBJECT(TabOpLowPassFilter, TableOperator);

public:
    OpenSim_DECLARE_PROPERTY(cutoff_frequency, double,
            "Low-pass cutoff frequency (Hz) (default is -1, which means no "
            "filtering).");
    OpenSim_DECLARE_PROPERTY(trim_to_original_time_range, bool,
            "Trim the rows of the output table to match the original table's "
            "time range after filtering (default: true).");
    TabOpLowPassFilter() { 
        constructProperty_cutoff_frequency(-1); 
        constructProperty_trim_to_original_time_range(true);
    }
    TabOpLowPassFilter(double cutoffFrequency) : TabOpLowPassFilter() {
        set_cutoff_frequency(cutoffFrequency);
    }
    TabOpLowPassFilter(double cutoffFrequency, bool trimToOriginalTimeRange) : 
            TabOpLowPassFilter() {
        set_cutoff_frequency(cutoffFrequency);
        set_trim_to_original_time_range(trimToOriginalTimeRange);
    }
    void operate(TimeSeriesTable& table, const Model* model = nullptr)
            const override {
        if (get_cutoff_frequency() != -1) {
            OPENSIM_THROW_IF(get_cutoff_frequency() <= 0, Exception,
                    "Expected cutoff frequency to be positive, but got {}.",
                    get_cutoff_frequency());

            const auto& times = table.getIndependentColumn();
            double initialTime = times.front();
            double finalTime = times.back();
            TableUtilities::filterLowpass(
                    table, get_cutoff_frequency(), true);
            if (get_trim_to_original_time_range()) {
                table.trim(initialTime, finalTime);
            }
        }
    }
};

You can observe, this uses the above method from Path 2, but adds padding by default and has an optional argument to trim the data back to its original length. Supposedly, the following is also an option, but I have gotten a segmentation fault when I’ve tried it:

TableProcessor proc = TableProcessor("file.sto") | TabOpLowPassFilter(6);

This method appears to be used mostly in Moco processing pipelines, but it is again attempting to achieve the same goal of low-pass filtering something in yet another slightly different way.

Filter Comparison Tests

The myriad of code paths left me befuddled. Principally, I wanted to verify that low-pass filtering in OpenSim was operating the same as in Python SciPy and Matlab.

To test this hypothesis, I created some simple tests to run Butterworth low-pass filtering at 3, 6, 12, and 20 Hz for the same data in all three platforms. Python was used to generate a perfect sine wave and then noise was added with a fixed random seed. The result was saved to a file and then loaded into each program to ensure the same signal was filtered. The signal was generated with the following code:

    frequency = 2  # Frequency of the sine wave in Hz
    amplitude = 1  # Amplitude of the sine wave
    sampling_rate = 100  # Samples per second
    duration = 1  # Duration in seconds

    t = np.linspace(0, duration, int(sampling_rate * duration))

    # Generate sine wave without noise
    sine_wave = amplitude * np.sin(2 * np.pi * frequency * t)

    # Fix the seed for reproducibility
    rng = np.random.default_rng(42)

    # Add Gaussian noise (mean, std)
    noise = rng.normal(0, 0.3, sine_wave.shape)
    sine_wave_with_noise = sine_wave + noise

To run the tests, I created a repository with the implementations for all 3 languages. The specific low-pass filter implementations are described below:

1. Python SciPy

The SciPy Python library performs Butterworth low-pass filtering as follows:

    normalized_cutoff = cutoff / (sampling_rate / 2)
    b, a = butter(filter_order, normalized_cutoff, btype="low")
    filtered_signal = filtfilt(b, a, noisy_signal)

2. Matlab

Matlab performs Butterworth low-pass filtering as follows:

    normalized_cutoff = cutoff / (sampling_rate / 2);
    [b, a] = butter(4, normalized_cutoff, 'low');
    filtered_signals(i, :) = filtfilt(b, a, noisy_signal);

3. OpenSim

For this experiment, the OpenSim low-pass filtering is tested as follows:

    OpenSim::TableUtilities::filterLowpass(upd_table, cutoff,true);

You can observe that the Matlab and Python APIs are identical whereas the OpenSim version is drastically different. Perhaps the OpenSim API could be implemented in a similar way in the future?

Results

Initially, I called the filterLowpass method in OpenSim without the mysterious 3rd boolean argument (which defaults to false). This proved to be a critical mistake.

Failure 1 - No Padding

Low-pass filter fail 3rd order

I learned that for the infinite impulse response (IIR) Butterworth filter to function properly, it needs an infinite amount of data. Since this is physically impossible, it is usually implemented by mirroring the signal in front of it and behind it to create a new “padded” signal. By default, the method does not do this and will produce significantly invalid data near the start and end of the signal.

Failure 2 - Wrong Filter Order

Low-pass filter comparison 4th order

Frequently lauded for its superior performance, the 4th order Butterworth low-pass filter with a frequency of 6 Hz is the filter of choice in many disciplines. Through visual acquiescence of the aforementioned figure, one can indubitably validate its superior qualities. While it stands victorious, unfortunately OpenSim currently only supports 3rd order filtering. While the differences between 3rd and 4th order are not strikingly grand, one can notice minor discrepancies between Python/Matlab and the OpenSim results.

Success - 3rd Order Filtering

Low-pass filter comparison 3rd order

After much experimentation, the signals finally aligned for the tested frequencies when all libraries were tested using 3rd order Butterworth low-pass filters. I am now satisfied in reporting that the 3rd order Butterworth low-pass implementations in Python, Matlab, and OpenSim produce the same results when the implicit assumptions of OpenSim are accepted. This highlights the importance of verifying that the programming implementations of methods and algorithms are properly functioning according to the underlying system assumptions.

Finally, supporting only 3rd order filtering is a major limitation in OpenSim’s implementation. The API would benefit from consolidating the disparate code paths into one, clear way to perform low-pass filtering. Additionally, following an API design similar or identical to existing implementations in Matlab and Python would improve clarity and understanding. Bespoke method development may seem promising during development, but it ultimately hinders understanding when it requires reverse engineering the implementation to verify its effectiveness.