我有一个力学问题的数值解。我正在计算施加在物体上的力作为其变形的函数,并且在问题中存在不稳定性,由于这种不稳定性,随着位移的增加(负载),力会突然从一种行为跳到另一种行为。同样,当位移减小(卸载)时,力再次跳跃,但由于摩擦,加载时的力曲线与卸载时不同。
我用来计算力的数值技术会在输出中引入噪声。我有兴趣平滑输出,以便与实验进行有意义的比较。由于我在 DSP 方面几乎没有背景,所以我只是在破解这个问题。我在 StackOverflow 上发现了一个类似的问题,并使用了在那里找到的代码。由于我想避免过滤后的输出滞后,因此我使用filtfilt
而不是lfilter
该答案中使用的那样。
此链接包含四个文件:
pauchard_F_u.png
显示了实验获得的力-位移曲线,该曲线显示了我在第一段中描述的特征;该图中的实心圆点对应于增加位移或加载部分,而空心圆对应于卸载部分。
force_displacement.rpt
包含作为时间函数的原始位移和力数据。在对应于 0.1 秒总持续时间的 1001 个等距时间点上可获得数据。mwe_filtering.py
包含一个 MWE,它进行过滤并绘制原始和过滤后的输出。这也包括在下面。trial_force_disp_smoothing.png
包含原始数据和平滑数据的图。在这些图中,黑色代表装载部分,蓝色代表卸载部分;符号代表原始数据,而曲线代表平滑输出。
该解决方案看起来不错,但是这里仍然有一些我不太喜欢的工件。
平滑力在加载部分的末端(即最大位移附近的黑色曲线)下降,而在原始数据(和实验)中,情况并非如此。第二次下降在物理上没有意义。
卸载曲线在零位移时明显偏离零(接近零位移的蓝色曲线)。同样,这在物理上没有意义,在零位移时,我们应该有零力。
在实验中确实观察到负载曲线中的第一次力下降,并在平滑数据中捕获,但过滤器平滑负载突然下降(黑线与黑色符号相比)有点太多了。如果我们能一直捕捉到负载下降点的原始曲线会更好。对于卸载(蓝色)曲线,情况并不那么明显。
问题:有人可以提出一种方法来避免平滑数据中的这些伪影并与实验数据取得更好的定性一致性吗?
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)