[Bug] Low Pass Filtering of Time

[alexbeattie42]

When using the Low Pass filter TabOpLowPassFilter the time column is also filtered/shifted.

A simple reproduction is available here in the LowPassFilterTime folder

This code is used in the example

auto tableProcessor =
            OpenSim::TableProcessor(coordinates_file_name) |
            OpenSim::TabOpLowPassFilter(6) |
            OpenSim::TabOpUseAbsoluteStateNames();
    auto coordinatesRadians = tableProcessor.processAndConvertToRadians(model);

When starting with the input data file that contains rows of this form:

time    pelvis_tilt pelvis_list pelvis_rotation pelvis_tx   pelvis_ty   pelvis_tz   hip_flexion_r   hip_adduction_r hip_rotation_r  knee_angle_r    ankle_angle_r   subtalar_angle_r    mtp_angle_r hip_flexion_l   hip_adduction_l hip_rotation_l  knee_angle_l    ankle_angle_l   subtalar_angle_l    mtp_angle_l lumbar_extension    lumbar_bending  lumbar_rotation
0   0.01068669390326228 0.02227099665722097 -23.8604585223402   0   0.95    0   -0.00247198488240275    -0.03048687110098423    0.002135644701990358    -0.007002076054592501   0.02117478415987187 0.008410610568575982    0   0.002375731643235115    0.02905620892918622 0.002304254953310731    0.005173967979601628    -0.02002546339977476    -0.0088756853696829 0   0   0   0
0.025   0.1420605864218891  0.1168697208777443  -23.89205323097089  0   0.95    0   0.07214083262762791 -0.2419351688763855 0.1610492605076858  -0.2643792819274658 0.1078648902039214  0.3187753232449059  0   -0.1647218090057443 0.1422293467617664  -0.6820281123974768 0.01164676605859697 0.1173102673695643  0.303094359792636   0   0   0   0
0.05    0.1420605864218891  0.1168697208777443  -23.89205323097089  0   0.95    0   0.07214083262762791 -0.2419351688763855 0.1610492605076858  -0.2643792819274658 0.1078648902039214  0.3187753232449059  0   -0.1647218090057443 0.1422293467617664  -0.6820281123974768 0.01164676605859697 0.1173102673695643  0.303094359792636   0   0   0   0
0.07500000000000001 0.1322986490885017  0.1201732180342828  -23.88929863760237  0   0.95    0   0.06633588355747383 -0.2445139197759301 0.1626747009690541  -0.268446881042728  0.1096859306392073  0.3159168223854317  0   -0.167505976939765  0.1432673537173976  -0.6874535109512012 0.00851513840238282 0.1174678435804044  0.3018689695747808  0   0   0   0

After running the above table processor with low-pass filtering it becomes:

time    /jointset/ground_pelvis/pelvis_tilt/value   /jointset/ground_pelvis/pelvis_list/value   /jointset/ground_pelvis/pelvis_rotation/value   /jointset/ground_pelvis/pelvis_tx/value /jointset/ground_pelvis/pelvis_ty/value /jointset/ground_pelvis/pelvis_tz/value /jointset/hip_r/hip_flexion_r/value /jointset/hip_r/hip_adduction_r/value   /jointset/hip_r/hip_rotation_r/value    /jointset/knee_r/knee_angle_r/value /jointset/ankle_r/ankle_angle_r/value   /jointset/subtalar_r/subtalar_angle_r/value /jointset/mtp_r/mtp_angle_r/value   /jointset/hip_l/hip_flexion_l/value /jointset/hip_l/hip_adduction_l/value   /jointset/hip_l/hip_rotation_l/value    /jointset/knee_l/knee_angle_l/value /jointset/ankle_l/ankle_angle_l/value   /jointset/subtalar_l/subtalar_angle_l/value /jointset/mtp_l/mtp_angle_l/value   /jointset/back/lumbar_extension/value   /jointset/back/lumbar_bending/value /jointset/back/lumbar_rotation/value
-12.15000000000011  -7.271659835996715  0.2526807067493628  -15.15721143271405  0   0.95    0   -2.760892184433378  -0.6271791115638784 0.8841443912025075  53.7439857124721    6.961750540897951   -21.93690907795083  0   -14.61827063823635  -6.776824211826309  -1.579267426163963  16.03636930005162   3.271257781401269   0.1649098762727437  0   0   0   0
-12.12500000000012  -7.325743401683875  0.3089371250227219  -15.07269135460541  0   0.9500000000000002  0   -2.204062995250402  -0.9561946108658287 0.6900905517703898  53.35390635425459   7.31932982891661    -22.61952115545578  0   -15.16053580077654  -6.778735867460504  -1.483030436742771  16.40478087320183   3.509599495807006   0.2380661807937111  0   0   0   0
-12.10000000000013  -7.363149948187758  0.4015403719472346  -14.99965330359036  0   0.9500000000000002  0   -1.730682902884436  -1.367459766304516  0.3998644921687521  52.92655186460718   7.537301309052862   -23.48053002248338  0   -15.67671652273414  -6.75083378741887   -1.35412892872716   16.7482220997471    3.766258726018362   0.2991510694091759  0   0   0   0
-12.07500000000014  -7.402813791766203  0.5851813714239973  -14.90088472785716  0   0.95    0   -0.9004002480721484 -1.975026070077736  -0.1943940979493001 52.06568810954337   7.926730690341419   -25.10962542159387  0   -16.44998216418404  -6.673104977819559  -1.160426551931182  17.1794484812042    4.24661367221694    0.406039046434963   0   0   0   0

If you run the same code and comment out the low-pass filter step, the time remains unchanged as expected. The low-pass filter is filtering or shifting the time column which shouldn't be happening.

Disabling padding by setting this line to false seems to fix the problem. That leads me to believe the issue is caused here where the time column is padded since here (where the padding is done) prepends data (so I suppose it goes backwards in time?)

[alexbeattie42]

This seems to be related to this issue

[nickbianco]

Yes, your issue is likely the same one in #3287. Been meaning to fix this but it kept slipping my mind.

Just submitted a fix at https://github.com/opensim-org/opensim-core/pull/3969.

[nickbianco]

@alexbeattie42 the fix at #3969 has been merged. Let me know if that resolves this issue.

[alexbeattie42]

@nickbianco I have just tested the fix and it seems to have solved it. The only problem I am noticing now is that the first data point is removed from the output. So I think the trimming is chopping one too many elements off the result after filtering.

[nickbianco]

@alexbeattie42 thanks for reporting. My guess is that this issue related to how TimeSeriesTable handles trimming. TimeSeriesTable::trim() computes indices based on the initial and final time provided:

https://github.com/opensim-org/opensim-core/blob/db5dcd45aa2160e8e143046a06e3a05de6b54151/OpenSim/Common/TimeSeriesTable.h#L467

But getRowIndexAfterTime or getRowIndexBeforeTime might grab and index just after (or before) the index for the original time, e.g.,

https://github.com/opensim-org/opensim-core/blob/db5dcd45aa2160e8e143046a06e3a05de6b54151/OpenSim/Common/TimeSeriesTable.h#L346

It would be helpful to take a look at what's happening with your data. Willing to take a look? If not, I could try to give it a shot in the next day or two.

[alexbeattie42]

Thanks for the guidance! It looks like it is a floating point precision issue. The data that I am using starts with uniformly sampled decimal time points every 0.025 seconds.

After processing with the TabOpLowPassFilter there are floating point precision errors in the time column. To me the expected behavior would be that the time column is unaltered since I do not want to resample the data.

I have determined the problem is that the filterLowpass is deciding to resample the data here.

The SimTK docs say "If Real is double (the normal case) then Eps ~= 1e-16". The result of dtAvg - dtMin is 8.20524e-15 which is larger than SimTK::Eps (which is 2.22045e-16 on my machine). So even though the file is already uniformly sampled, the code silently tries to resample it since the determined sampling rate error is 8.20524e-15 > 2.22045e-16.

Manually searching for where 0 should be yields that the 0 time value is now -4.2543746303636e-12. Resulting data after filtering:

-4.2543746303636e-12    0.0001865179945371941   0.0003887022192010853   -0.4164435622503228 0   0.9499999999999995  0   -4.314427538659543e-05  -0.0005320962789695338  3.727403142175347e-05   -0.0001222092811723213  0.0003695697018353934   0.0001467928459787541   0   4.146433961397945e-05   0.0005071265137617322   4.021683694910552e-05   9.030277662670975e-05   -0.0003495102708057906  -0.000154909933621432   0   0   0   0
0.0249999999957371  0.001444431434548449    0.001321552278548245    -0.4167450479826717 0   0.9499999999999995  0   0.0006707204198965602   -0.002604717315420156   0.001594071484128392    -0.002647996736714108   0.001221114680527505    0.003172171719565231    0   -0.001598795396258941   0.00161591734612298 -0.006660622999522429   0.0001451134537261808   0.0009925676347088941   0.00288996565013442 0   0   0   0
0.04999999999572857 0.002297312405166353    0.001972911593669966    -0.4169484981715055 0   0.9499999999999994  0   0.00115427729936877 -0.00404274258526328    0.002673666564447302    -0.004403017216715688   0.001813219179243035    0.00526057732241384 0   -0.00273877801207574    0.002385180873038865    -0.01130663935391225    0.0001771377685193153   0.001921152994475303    0.004993938565327437    0   0   0   0

The things that I believe would be good to include are:

1. A warning or log message to the user that the data is being resampled at the code location above
2. The ability to disable automatic resampling of the data
3. A less aggressive default tolerance for the sampling interval (or the ability to change it as a parameter). For example SimTK::SqrtEps should be ~1e-8 if Real is double which is suitable for this case
4. Documentation updates about the default behavior of the lowpass filter. I was unable to find any documentation that specifies that by default the lowpass filter mirror pads the signal, resamples it, and then chops off the padding at the end. This would be very useful information to know when deciding to use the filter because it improves the credibility and trustworthiness of the filter implementation for the end user.

[nickbianco]

Nice job finding the root issue @alexbeattie42! I largely agree with proposed changes.

We will try to prioritize this in our development plans since this likely affects a lot of users. If you're interested in submitting a PR yourself, I'd be happy to review it.