Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Bug] Low Pass Filtering of Time #3968

Open
alexbeattie42 opened this issue Nov 15, 2024 · 7 comments · May be fixed by #3978
Open

[Bug] Low Pass Filtering of Time #3968

alexbeattie42 opened this issue Nov 15, 2024 · 7 comments · May be fixed by #3978
Assignees

Comments

@alexbeattie42
Copy link
Contributor

alexbeattie42 commented Nov 15, 2024

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
Copy link
Contributor Author

This seems to be related to this issue

@nickbianco
Copy link
Member

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 #3969.

@nickbianco nickbianco self-assigned this Nov 19, 2024
@nickbianco
Copy link
Member

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

@alexbeattie42
Copy link
Contributor Author

alexbeattie42 commented Nov 20, 2024

@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
Copy link
Member

@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:

start_index = this->getRowIndexAfterTime(newStartTime);

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

if (timeCol[candidate] < time-eps)

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
Copy link
Contributor Author

alexbeattie42 commented Nov 21, 2024

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
Copy link
Member

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants