如何平滑这个嘈杂的信号?

信息处理 低通滤波器 平滑
2022-02-02 11:08:19

我有一个力学问题的数值解。我正在计算施加在物体上的力作为其变形的函数,并且在问题中存在不稳定性,由于这种不稳定性,随着位移的增加(负载),力会突然从一种行为跳到另一种行为。同样,当位移减小(卸载)时,力再次跳跃,但由于摩擦,加载时的力曲线与卸载时不同。

我用来计算力的数值技术会在输出中引入噪声。我有兴趣平滑输出,以便与实验进行有意义的比较。由于我在 DSP 方面几乎没有背景,所以我只是在破解这个问题。我在 StackOverflow 上发现了一个类似的问题,并使用了在那里找到的代码。由于我想避免过滤后的输出滞后,因此我使用filtfilt而不是lfilter该答案中使用的那样。

链接包含四个文件:

  1. pauchard_F_u.png显示了实验获得的力-位移曲线,该曲线显示了我在第一段中描述的特征;该图中的实心圆点对应于增加位移或加载部分,而空心圆对应于卸载部分。

在此处输入图像描述

  1. force_displacement.rpt包含作为时间函数的原始位移和力数据。在对应于 0.1 秒总持续时间的 1001 个等距时间点上可获得数据。
  2. mwe_filtering.py包含一个 MWE,它进行过滤并绘制原始和过滤后的输出。这也包括在下面。

  3. trial_force_disp_smoothing.png包含原始数据和平滑数据的图。在这些图中,黑色代表装载部分,蓝色代表卸载部分;符号代表原始数据,而曲线代表平滑输出。

在此处输入图像描述

该解决方案看起来不错,但是这里仍然有一些我不太喜欢的工件。

  1. 平滑力在加载部分的末端(即最大位移附近的黑色曲线)下降,而在原始数据(和实验)中,情况并非如此。第二次下降在物理上没有意义。

  2. 卸载曲线在零位移时明显偏离零(接近零位移的蓝色曲线)。同样,这在物理上没有意义,在零位移时,我们应该有零力。

  3. 在实验中确实观察到负载曲线中的第一次力下降,并在平滑数据中捕获,但过滤器平滑负载突然下降(黑线与黑色符号相比)有点太多了。如果我们能一直捕捉到负载下降点的原始曲线会更好。对于卸载(蓝色)曲线,情况并不那么明显。

问题:有人可以提出一种方法来避免平滑数据中的这些伪影并与实验数据取得更好的定性一致性吗?


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)
2个回答

我认为非线性滤波器会给你一个更接近你期望的结果。最简单的非线性滤波器之一是中值滤波器,它可以在保留跳跃的同时去除噪声。下图显示了您的数据和窗口长度为的中值滤波器的输出20样品。可以增加窗口长度以获得更平滑的效果,当然您可以为两条曲线使用不同的窗口长度。

在此处输入图像描述

filtfilt 的问题在于,此过程假设可能非常不正确的数据(零?)超过实际数据的两端。您可以通过多种不同的方式处理此问题。

在末端使用与数据向量中间不同的过滤器(可能是中位数?),然后交叉淡化结果。

或者,如果您喜欢当前的过滤器,您可以在每一端制作“可能不那么不正确”的数据,以便线性相位过滤器过程使用。

在每一端镜像数据(至少为线性相位滤波器的瞬态响应长度),然后在滤波后切断那些添加的端。

或者使用超定多项式回归在每一端拟合曲线(比如使用 8 个左右的点来拟合三次)。然后通过使用经过末端的拟合曲线外推(通过过滤器的长度)来弥补点。然后过滤包括添加的末端,然后丢弃这些末端。