3
$\begingroup$

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.

enter image description here

  1. 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.
  2. mwe_filtering.py contains a MWE that does the filtering and plots the raw and filtered outputs. This is also included below.

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

enter image description here

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)
$\endgroup$
  • $\begingroup$ Welcome to DSP.SE! Nice question! It might take a little to work through. $\endgroup$ – Peter K. Dec 5 '15 at 2:22
  • $\begingroup$ @MattL. The RPT file in the Google Drive link seems to be a simple text file. $\endgroup$ – Peter K. Dec 5 '15 at 11:50
2
$\begingroup$

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.

enter image description here

$\endgroup$
  • $\begingroup$ this seems to take me in the right direction, but it doesn't get me all the way. First, it clips the real peak in the loading curve; likewise there is a real peak on the unloading (around x = 0.02 in your plot), but I can live with this not being captured after filtering. Next, the filter does not help with the shape at the end of the loading $\endgroup$ – Sankara Subramanian Dec 5 '15 at 20:20
  • $\begingroup$ @SankaraSubramanian: There's a trade-off between denoising and capturing each and every detail of the data. What's wrong with the shape at the end of the loading curve? There is no dip as in your smoothed curve. $\endgroup$ – Matt L. Dec 5 '15 at 20:24
  • $\begingroup$ in your plot, actually there is not much wrong with the shape at the end of the loading. But, if I increase the window size to get more smoothing, the curve at the end of the loading begins to drop as it begins to get influenced by the unloading data. And I do want more smoothing than in your plot. $\endgroup$ – Sankara Subramanian Dec 5 '15 at 20:30
  • $\begingroup$ It occurs to me that one might treat these data as a combination of four different smaller ones: smooth loading, noisy loading, noisy unloading and smooth unloading. Guided by the physics, I am tempted to retain the first and fourth data pieces as they are, without any filtering of just minimal filtering and perhaps, filter only the middle two sections ensuring that the load does not drop at the end of the loading section. $\endgroup$ – Sankara Subramanian Dec 5 '15 at 20:33
  • $\begingroup$ @SankaraSubramanian: I think you have to smooth the loading and unloading data separately. $\endgroup$ – Matt L. Dec 5 '15 at 20:33
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ trying out some of these. newbie question; how do i 'cross-fade' the results? $\endgroup$ – Sankara Subramanian Dec 7 '15 at 2:52
  • $\begingroup$ is there a way to weight points differently in the filtering process if I have different levels of confidence in their accuracy? $\endgroup$ – Sankara Subramanian Dec 7 '15 at 2:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.