使用 Numpy 将高斯滤波器应用于一维数据 "by hands"

Question

我有一个非均匀采样数据，我正尝试对其应用高斯滤波器。我正在使用 python 的 numpy 库来解决这个问题。数据为XY类型，如下所示：

[[ -0.96       390.63523024]
[ -1.085      390.68523024]
[ -1.21       390.44023023]
...
[-76.695      390.86023024]
[-77.105      392.51023024]
[-77.155      392.10023024]]

而here是整个*.npz文件的link。

这是我的方法：

我从定义高斯函数开始
然后我开始使用 while 循环沿 X 轴扫描数据
在循环的每一步中：
- I select 两个截止长度内的一部分数据
- 移动 selected 数据部分的 X 轴，使其围绕 0
- 在每个点计算我的高斯函数，乘以相应的Y值，求和并除以元素数
移至下一点

代码如下所示：

import numpy as np
import matplotlib.pyplot as plt

xy = np.load('1D_data.npz')['arr_0']

def g_func(xx, w=1.0):
    a = 0.47 * w
    return (1 / a) * np.exp((xx / a) ** 2 * (-np.pi))

x, y, x_, y_ = xy[:, 0], xy[:, 1], [], []

counter, xi, ww = 0, x[0], 1.0
while xi > np.amin(x):

    curr_x = x[(x < xi) & (x >= xi - 2 * ww)]
    g, ysel = [], []
    for i, els in enumerate(curr_x):
        xil = els - curr_x[0] + abs(curr_x[0] - curr_x[-1]) / 2
        g.append(g_func(xil, ww))
        ysel.append(y[counter + i])

    y_.append(np.sum(np.multiply(g, ysel)) / len(g))
    x_.append(xi)

    counter += 1
    xi = x[counter]

plt.plot(x, y, '-k')
plt.plot(x_, y_, '-r')
plt.show()

虽然输出看起来不正确。（见下图）即使丢弃边缘，卷积也非常嘈杂，值似乎与数据不对应。我可能做错了什么？

Answer 1

您在代码中犯了一个错误：

g乘y_sel前，y_sel不居中

之所以y_sel要居中，是因为我们要将高斯加权的相对差异添加到中心的条目中。如果直接将 g 与 y_sel 相乘，不仅 window 内的相邻条目的值，而且中心条目的值将由高斯加权。这肯定会显着改变函数值。

下面是我使用 numpy

的解决方案

def g_func(xx, w=1.0):
    mean = np.mean(xx)
    a = 0.47 * w
    return (1 / a) * np.exp(((xx-mean) / a) ** 2 * (-np.pi))


def get_convolution(array,half_window_size):
    
    array = np.concatenate((np.repeat(array[0],half_window_size),
                            array,
                            np.repeat(array[-1],half_window_size)))
    window_inds = [list(range(ind-half_window_size,ind+half_window_size+1)) \
                   for ind in range(half_window_size,len(array)-half_window_size)]
    
    return np.take(array,window_inds)

xy = np.load('1D_data.npz')['arr_0']
x, y = xy[:, 0], xy[:, 1]

half_window_size = 4
x_conv = np.apply_along_axis(g_func,axis=1,arr=get_convolution(x,half_window_size=half_window_size))
y_conv = get_convolution(y,half_window_size=half_window_size)
y_mean = np.mean(y_conv,axis=1)
y_centered = y_conv - y_mean[:,None]
smoothed = np.sum(x_conv*y_centered,axis=1) / (half_window_size*2) + y_mean

fig,ax = plt.subplots(figsize=(10,6))
ax.plot(x, y, '-k')
ax.plot(x, smoothed, '-r')

运行代码，输出为

更新

为了统一w和half_window_size，这里有一种可能，思路是让高斯的标准差为2*half_window_size

def g_func(xx):
    
    std = len(xx)
    mean = np.mean(xx)
    return 1 / (std*np.sqrt(2*np.pi)) * np.exp(-1/2*((xx-mean)/std)**2)


def get_convolution(array,half_window_size):
    
    array = np.concatenate((np.repeat(array[0],half_window_size),
                            array,
                            np.repeat(array[-1],half_window_size)))
    window_inds = [list(range(ind-half_window_size,ind+half_window_size+1)) \
                   for ind in range(half_window_size,len(array)-half_window_size)]
    
    return np.take(array,window_inds)

xy = np.load('1D_data.npz')['arr_0']
x, y = xy[:, 0], xy[:, 1]

half_window_size = 4
x_conv = np.apply_along_axis(g_func,axis=1,arr=get_convolution(x,half_window_size=half_window_size))
y_conv = get_convolution(y,half_window_size=half_window_size)
y_mean = np.mean(y_conv,axis=1)
y_centered = y_conv - y_mean[:,None]
smoothed = np.sum(x_conv*y_centered,axis=1) / (half_window_size*2) + y_mean

fig,ax = plt.subplots(figsize=(10,6))
ax.plot(x, y, '-k')
ax.plot(x, smoothed, '-r')

使用 Numpy 将高斯滤波器应用于一维数据 "by hands"

Applying Gaussian filter to 1D data "by hands" using Numpy

python

numpy

convolution