探索【一维插值法】高还原度保留细节

15195775117

1、"还原度"不高的scipy.interpolate.interp1d

工作中常有将折线平滑的需求，如果数据量比较大，插值法看起来像那么回事，
但当数据量较小的时候，就会发现某些插值法与原始数据有差距，会影响我们后续的识别计算，

刚开始我用的一维插值法是scipy.interpolate.interp1d，
效果如下：

观察之，发现这种插值法并不合适，
因为，一般我们希望：
原始折线上的峰值在插值曲线上也是峰值，
原始折线上的谷值在插值曲线上也是谷值，
但上图中，
原始的点并不在插值曲线上峰谷的位置，
插值曲线上的峰谷位置也没有原始点。
所以，插值曲线对原始序列的还原度不高。

上图源代码：
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
# 注意，interp1d中间的是数字1，不是字母L！

# 画图数据：
x=np.linspace(0,1,10) ; y=np.random.random(10)

# 根据a,b的关系生成个函数，类似一个拟合函数,kind表示样条插值的幂，一般用3就挺好
f=interp1d(x,y,kind=5)

x2=np.linspace(0,1,100) # 将x加密，生成x2
y2=f(x2) # 根据x2和f计算出y2

# 画图看效果：
plt.figure(figsize=(6,3),dpi=200)
plt.scatter(x,y,marker="o",s=200)
plt.plot(x2,y2,"r")
plt.show()

15195775117

2、高还原度的算法

思路：先把原始序列加密成阶梯型序列，例如：
[1,2,3]-->[1,1,1,2,2,2,3,3,3]
再将阶梯序列用高斯平滑（scipy.ndimage.gaussian_filter1d）处理一下。

下图中，红色圆点是原始数据，蓝色阶梯折线是加密的阶梯数据，红色曲线是高斯平滑后的数据。

可见，原始数据都在插值曲线的峰谷上，比较符合一般认知。



源代码：

'''
该函数实现“高还原度”的一维平滑插值
'''
import numpy as np
from scipy.ndimage import gaussian_filter1d

def detailSmooth1D(x,y,accuracy):
    # 不断对x和y加密，直到到达精度要求为止
    for repeatTimes in np.arange(2,200,2):

        # 序列加密：

        # “散点数据”y-->“阶梯数据”y2：
        # 注意：首尾的台阶只有一半
        y2=[]
        y2.extend(y[0].repeat(repeatTimes/2))
        y2.extend(y[1:-1].repeat(repeatTimes))
        y2.extend(y[-1].repeat(repeatTimes/2))

        # 将x序列加密：
        x2= np.linspace(np.min(x),np.max(x),len(y2))

        # 平滑y2：
        y3= gaussian_filter1d(y2,3)

        tag=True
        # 一旦发生某个插的值与原始值差距超过accuracy，
        # tag就赋予False，此次循环停止，切入下个循环值，
        # 如果一直没有失败，tag依旧是True，这时获得的结果就可以返回了，
        # repeatTimes的循环也就该停止了

        # 对于每一段“阶梯”，求插的值与原始值的差距，
        # 如果差距都在误差限内，则通过考验
        for i in range(0,len(y)-1):
            # 该阶梯范围的y2的值应该是一样的：
            # 阶梯的起始位置：
            pos1=int((i+0.5)*repeatTimes)
            pos2=int((i+1.5)*repeatTimes)
            stage_y2=y2[pos1:pos2][0]
            # 该阶梯范围的y3的值：
            stage_y3=y3[pos1:pos2]
            # y3的值与y2的值的差的最小值：
            df=np.abs(stage_y3-stage_y2)
            df=np.min(df)
            # 一旦出现某一个插出来的值与原始值差距超出误差限，
            # 就没必要再算了，宣判失败：
            if np.abs(df)>accuracy:
                tag=False
                break

        # 如果通过考验，就没必要再加密了，上层循环停止
        if tag:
            print('最小倍数=',repeatTimes)
            break
    return {'x':x2,'y':y3}

if __name__=="__main__":
    import matplotlib.pyplot as plt
    from datetime import datetime
    t1=datetime.now()

    # 准备参数x,y,accuracy
    y=np.arange(10)
    y=np.array(y,dtype=float)
    np.random.shuffle(y) # 洗牌打乱
    x=np.arange(len(y))+100

    # 插值后的值与实际值的差小于accuracy
    accuracy=0.1

    result=detailSmooth1D(x,y,accuracy)
    x2=result['x']
    y2=result['y']

    plt.plot(x,y,'bo')
    plt.plot(x2,y2,'r')
    plt.show()

    t2=datetime.now()
    dt=t2-t1
    print('运行时间=%0.1f'%(dt.microseconds/1000),'毫秒')

由于算法是通过阶梯曲线来平滑，所以accuracy越小，曲线的峰谷越平：