I have a numerical solution to a problem in mechanics. I am computing the force applied to an object as a function of its deformation and in the problem, there is an instability, due to which the force abruptly jumps from one behaviour to another as the displacement is increased (loading). Likewise, when the displacement is decreased (unloading), the force again jumps, but due to friction, the force curve during loading is not the same as unloading.
The numerical technique I use to compute the force introduces noise into the output. I am interested in smoothing the output so that I can have a meaningful comparison with experiment. Since I have little to no background in DSP, I am just hacking my way through this. I found a similar question on StackOverflow and I used the code I found there. Since I wanted to avoid lag in the filtered output, I used filtfilt
instead of lfilter
as used in that answer.
This link contains four files:
pauchard_F_u.png
shows an experimentally obtained force-displacement curve which shows the features I describe in the first paragraph; the solid dots in this figure correspond to increasing displacement or the loading portion, while the hollow circles correspond to the unloading portion.
force_displacement.rpt
contains the raw displacement and force data as functions of time. Data is available at 1001 equally spaced time points corresponding to a total time duration of 0.1 seconds.mwe_filtering.py
contains a MWE that does the filtering and plots the raw and filtered outputs. This is also included below.trial_force_disp_smoothing.png
contains the plot of the raw and smoothed data. In these plots, black represents the loading portion and blue the unloading portion; symbols represent the raw data while curves represent smoothed output.
The solution looks okay, but, there are a still a couple of artifacts here that I don't really like.
The smoothed force dips at the end of the loading section (i.e. black curve near maximum displacement) whereas in the raw data (and the experiment), this is not so. This second dip is not physically meaningful.
The unloading curve deviates significantly from zero at zero displacement (blue curve near zero displacement). Again, this is not physically meaningful, at zero displacement, we should have zero force.
The first dip in force in the loading curve is indeed observed in the experiment and is captured in the smoothed data, but the filter smooths the abrupt load drop (black line compared to black symbols) a bit too much. It would be better if we can capture the raw curve all the way up to the point at which the load drops. The situation is not so pronounced for the unloading (blue) curve.
Question: Can someone suggest a way to that avoids these artifacts in the smoothed data and achieves a better qualitative agreement with the experimental data?
mwe_filtering.py
""" Created: 11/6/2015
Butterworth filtering and filtfilt modeled after these two sources
https://stackoverflow.com/questions/25191620/creating-lowpass-filter-in-scipy-understanding-methods-and-units
http://dsp.stackexchange.com/questions/19084/applying-filter-in-scipy-signal-use-lfilter-or-filtfilt
Note that use of filtfilt removes the lag that lfilter introduces
"""
def butter_lowpass(cutoff, fs, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = butter(order, normal_cutoff, btype='low', analog=False)
return b, a
def butter_lowpass_filtfilt(data, cutoff, fs, order=5):
b, a = butter_lowpass(cutoff, fs, order=order)
y = filtfilt(b, a, data)
return y
if __name__ == "__main__":
import copy
import glob
from matplotlib import pyplot as plt
import pandas
from scipy.signal import argrelmax, butter, filtfilt, freqz, lfilter
import sys
# MAGIC NUMBERS
RADIUS = 25.0
THICKNESS = 0.790569415042095
ORDER = 6
FS = 10000
CUTOFF = 70
radius = RADIUS
thickness = THICKNESS
thickness_normalized = np.array(thickness) / radius
legend_string = r"$t/R = " + str(np.round(thickness_normalized)) + "$"
# For the Butterworth filter
order = ORDER
fs = FS
cutoff = CUTOFF
b, a = butter_lowpass(cutoff, fs, order)
rpt_file = 'force_disp.rpt'
my_df = pandas.read_csv(rpt_file, delim_whitespace = True)
disp = my_df.iloc[:,1]
forc = my_df.iloc[:,2]
forc_filtered = butter_lowpass_filtfilt(forc, cutoff, fs, order)
plt.figure()
plt.hold(True)
plt.plot(disp[:500:3], forc[:500:3],
marker = "o",
ms = 5,
ls = ':',
color = 'black',
label = legend_strings[count],)
plt.plot(disp[500::3], forc[500::3],
marker = "o",
ms = 5,
ls = ':',
mec = "blue",
mfc = "white",
label = legend_strings[count],)
plt.plot(disp[:500:3], forc_filtered[:500:3],
ls = '-',
color = 'black',
label = legend_strings[count],)
plt.plot(disp[500::3], forc_filtered[500::3],
ls = '-',
color = "blue",
label = legend_strings[count],)
# plt.legend(loc = 2, frameon = False, fontsize = 'xx-large')
plt.xlabel(r'$u/h$', fontsize = 20)
plt.ylabel(r'$F R/ E h^3$', fontsize = 20)