How do I smooth this noisy signal?

Question

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:

1. 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.

2. 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.

3. mwe_filtering.py contains a MWE that does the filtering and plots the raw and filtered outputs. This is also included below.

4. 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.

1. 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.

2. 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.

3. 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)

Answer 1

I think that a non-linear filter will give you a result that comes closer to what you expect. One of the simplest non-linear filters is a median filter, which can remove noise while preserving jumps. The figure below shows your data and the outputs of a median filter with a window length of $20$ samples. The window length can be increased to get more smoothing, and of course you can use different window lengths for both curves.

Answer 2

The problem with filtfilt is that this process is assuming potentially very incorrect data (zeros?) past the two ends of your actual data. You can deal with this problem a number of different ways.

Use a different filter at the ends than in the middle of your data vector (perhaps median?), then cross-fade the results.

Or, if you like your current filter, you can make up "perhaps less incorrect" data past each end for the linear phase filter process to work with.

Mirror the data at each end (for at least the length of the transient response of your linear phase filter), and then cut off those added ends after filtering.

Or use over-determined polynomial regression to fit a curve (say use 8 or so points to fit a cubic) at each end. Then make up points by extrapolating using that fitted curve past the ends (by the length of the filter). Then filter including the added ends and afterwards throw away these ends.