python 数据平滑_数据平滑方法的原理和应用

一、简介

在实际的工程应用中,经常会遇到初始结果噪声太多的问题,比如信号强度抖动的太厉害,比如视频流中的bbox抖动的太厉害,比如光谱信号抖动的太厉害等等,这时候就需要一些简单的滑动平均算法。滑动平均其实是一个很朴素的方法,但是要与实际结合,构造出合适的平滑方式,是需要一些思考的。下面我将分别介绍滑动平均法(Moving Average)、指数滑动平均法(Exponential Mean Average)、SG滤波法(Savitzky Golay Filter)。

二、滑动平均法

简单来说,滑动平均法把前后时刻的一共2n+1个观测值做平均,得到当前时刻的滤波结果。这是一个比较符合直觉的平滑方法,在生活中、工作中很经常会用到,但是很少去思考这么做的依据是什么,下面我就来仔细分析一下其中的原理。

对于一个观测序列,我们有这样的假设:每一次的观测值是带有噪声的,而我们期望噪声的均值为0,方差为

,观测值和真实值之间的关系如下:

(1)

其中,

为观测值,

为真实值,

为噪声。为了降低噪声的影响,我们把相邻时刻的观测值相加后平均,公式如下:

(2)

表示

时刻的滤波结果,

表示

时刻的观测值,

代表滑动窗口半径。将公式(1)代入公式(2),可以得到

(3)

前面说到了,我们假设噪声的均值为0,所以

为0,那么我们得到的结果就是:

(4)

当观测数据的真实值变化较小时,或者变化为线性时,可以近似认为:

(5)

从上面的分析过程我们可以看到,当滑动窗口内的真实数据变化不大的时候,我们可以抑制掉很大一部分噪声,滤波结果近似真实值;当滑动窗口内的真实值变化较大时,这种滤波方式就会损失一部分精确度,滤波结果接近真实值的平均期望。所以,窗口的大小会对滤波结果有很大影响。窗口越大,滤波结果越平滑,但会一定程度上偏离真实值;窗口越小,滤波结果越接近观测值,但噪声偏大。

滑动平均法还有一个升级版本,也就是加权滑动平均法。实际场景中,每个观测值的重要程度不同,忽略每个观测值的置信度直接平均不能得到精确的结果,所以就需要给观测值加权。加权滑动平均法的公式如下:

(6)

时刻的权重。(6)式表示的是把每个观测值乘以权重后再平均。这种方法适用于观测值本身带有置信度的情况。注意,这里有一个小问题,如果置信度的取值范围是0到1之间,那么加权之后计算得到的观测值往往小于真实值,我来解释一下为什么。首先,我们假设观测值和真实值的均值是相等的,也就是

(7)

当我们把观测值乘以权重了之后,观测值和真实值的均值就不相等了,因为真实值的权重均值为1,而观测值的权重均值为

,是小于等于1的,最终的预测值也是小于等于真实值的,而且大概率是小于。所以我们需要对(6)式增加一个修正:

(8)

这样,得到的预测值就会更加合理了。

小结:滑动平均法使用的前提是,噪声的均值为0,真实值变化不大或线性变化的场景。如果真实值有较高频率的非线性突变的话,滑动平均法的效果就不够好了。同时,滑动平均法的窗口选取很重要,需要根据具体数据来选择。如果需要使用在线版本的滑动平均,那么就要把窗口前移,也就是把当前时刻的前n个观测值进行平均,但这样得到的结果会明显滞后于当前观测值,窗口越大,滞后的现象越严重。

class MovAvg(object):

def __init__(self, window_size=7):

self.window_size = window_size

self.data_queue = []

def update(self, data):

if len(self.data_queue) == self.window_size:

del self.data_queue[0]

self.data_queue.append(data)

return sum(self.data_queue)/len(self.data_queue)鼠标轨迹的滑动平均效果一维数据的滑动平均效果

三、指数滑动平均法

指数滑动平均法相当于加权滑动平均法的变体,主要区别在于,指数滑动平均法的权重是固定的、随时间推移呈指数衰减。指数滑动平均法的公式如下:

(9)

表示预测值,

表示衰减权重,通常我们设为固定值0.9,

表示观测值,这是一个递推公式。前面说了,指数滑动平均法的权重是随时间推移呈指数衰减的,那么上面的这个递推公式的指数体现在哪里呢?我们把(9)式进行延伸:

(10)

将(9)和(10)两式子联立,可得

(11)

发现没有,在(11)式中

的关系是

倍,而在(9)式中

的关系是

倍,呈指数衰减关系。同时,在初始时刻有如下关系:

(12)

根据这一关系和上述的递推公式,我们就能够得到整个算法的公式了。

由于这种指数衰减的特性,指数滑动平均法会比滑动平均法的实时性更强,更加接近当前时刻的观测值。在实际场景下,如果目标的波动较大时,指数滑动平均法会比滑动平均更加接近当前的真实值。那么是不是就说明,指数滑动平均法在任意场景下都比滑动平均法更好呢?不一定。我们来分析一下指数衰减法的误差项,这里为了简便表示,设定

,同时,将(1)式和(12)式代入公式(11),可得到误差项:

(13)

所以误差项也是呈指数衰减的,越接近当前时刻的误差项权重越大。假如在当前的工程场景中,误差是固定的分布,不受目标的观测值大小影响的话,那么指数滑动平均法会更接近真实值;假如误差会受目标观测值影响,比如我们观测的是一个连续运动的目标,中间突然出现了一个偏离很远的观测点,那么这个点为误检的概率相当大,也就是该观测值的误差比之前其他点的误差要大得多,此时指数加权平均法的结果就会波动较大,结果就不如滑动平均了。

小结:当误差不受观测值大小影响的话,指数滑动平均比滑动平均好;当误差随观测值大小变化时,滑动平均比指数滑动平均更好。

class ExpMovAvg(object):

def __init__(self, decay=0.9):

self.shadow = 0

self.decay = decay

self.first_time = True

def update(self, data):

if self.first_time:

self.shadow = data

self.first_time = False

return data

else:

self.shadow = self.decay*self.shadow+(1-self.decay)*data

return self.shadow鼠标轨迹的指数滑动平均效果一维数据的指数滑动平均效果

四、SG滤波法

SG滤波法(Savitzky Golay Filter)的核心思想也是对窗口内的数据进行加权滤波,但是它的加权权重是对给定的高阶多项式进行最小二乘拟合得到。它的优点在于,在滤波平滑的同时,能够更有效地保留信号的变化信息,下面我来介绍一下其原理。

我们同样对当前时刻的前后一共2n+1个观测值进行滤波,用k-1阶多项式对其进行拟合。对于当前时刻的观测值,我们用下面的公式进行拟合:

(14)

同样,对于前后时刻(如t-1、t+1、t-2、t+2等时刻)的预测值,我们同样可以用(14)式来计算,这样一共得到2n+1个式子,构成一个矩阵(似乎发不了矩阵,我放个图片吧):

要使得整个矩阵有解,必须满足 2n+1>k,这样我们才能够通过最小二乘法确定参数

...

。我们把上面的矩阵简化表示为下面公式:

(15)

各个参数下标表示它们各自的维度,如

表示有k行1列的参数。通过最小二乘法,我们可以求得

的解为:

(16)

上标trans表示转置。那么,模型的滤波值为:

(17)

最终可以得到滤波值和观测值之间的关系矩阵:

(18)

算出了B矩阵,我们就能够快速的将观测值转换为滤波值了。

小结:SG滤波法对于数据的观测信息保持的更好,在一些注重数据变化的场合会比较适用。

class SavGol(object):

def __init__(self, window_size=11, rank=2):

assert window_size % 2 == 1

self.window_size = window_size

self.rank = rank

self.size = int((self.window_size - 1) / 2)

self.mm = self.create_matrix(self.size)

self.data_seq = []

def create_matrix(self, size):

line_seq = np.linspace(-size, size, 2*size+1)

rank_seqs = [line_seq**j for j in range(self.rank)]

rank_seqs = np.mat(rank_seqs)

kernel = (rank_seqs.T * (rank_seqs * rank_seqs.T).I) * rank_seqs

mm = kernel[self.size].T

return mm

def update(self, data):

self.data_seq.append(data)

if len(self.data_seq) > self.window_size:

del self.data_seq[0]

padded_data = self.data_seq.copy()

if len(padded_data) < self.window_size:

left = int((self.window_size-len(padded_data))/2)

right = self.window_size-len(padded_data)-left

for i in range(left):

padded_data.insert(0, padded_data[0])

for i in range(right):

padded_data.insert(

len(padded_data), padded_data[len(padded_data)-1])

return (np.mat(padded_data)*self.mm).item()鼠标轨迹的SG滤波效果一维数据的SG滤波效果

附录

本文图片制作的相关代码。

import cv2

import numpy as np

import matplotlib.pyplot as plt

import imageio

# 一维数据滤波

ma, ema, sg = MovAvg(), ExpMovAvg(), SavGol()

data_list, data_ma, data_ema, data_sg = [], [], [], []

for i in range(200):

data = i+np.random.randint(-50, 50)

data_list.append(data)

data_ma.append(ma.update(data))

data_ema.append(ema.update(data))

data_sg.append(sg.update(data))

plt.plot(data_list, label='raw')

plt.plot(data_ma, label='ma')

plt.plot(data_ema, label='ema')

plt.plot(data_sg, label='sg')

plt.legend()

plt.show()

# 鼠标轨迹滤波

ma_x, ma_y = MovAvg(), MovAvg()

ema_x, ema_y = ExpMovAvg(), ExpMovAvg()

sg_x, sg_y = SavGol(), SavGol()

def draw_circle(event, x, y, flags, param):

if event == cv2.EVENT_MOUSEMOVE:

sx = np.random.randint(-50, 51)

sy = np.random.randint(-50, 51)

cv2.circle(show, (x+sx, y+sy), 5, (255, 255, 255), -1)

x, y = ma_x.update(x+sx), ma_y.update(y+sy)

cv2.circle(show, (int(x), int(y)), 5, (0, 0, 255), -1)

x, y = ema_x.update(x+sx), ema_y.update(y+sy)

cv2.circle(show, (int(x), int(y)), 5, (0, 255, 0), -1)

x, y = sg_x.update(x+sx), sg_y.update(y+sy)

cv2.circle(show, (int(x), int(y)), 5, (255, 0, 0), -1)

show = np.zeros((1024, 1024, 3), np.uint8)

cv2.namedWindow('image')

buff = []

while True:

cv2.setMouseCallback('image', draw_circle)

cv2.imshow('image', show)

save = cv2.resize(show, (512, 512))

save = cv2.cvtColor(save, cv2.COLOR_BGR2RGB)

buff.append(save)

if cv2.waitKey(100) == ord('q'):

break

cv2.destroyAllWindows()

imageio.mimwrite('test.gif', buff, 'GIF', duration=0.1)

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/472043.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

LeetCode 1920. 基于排列构建数组

文章目录1. 题目2. 解题1. 题目 给你一个 从 0 开始的排列 nums&#xff08;下标也从 0 开始&#xff09;。 请你构建一个 同样长度 的数组 ans &#xff0c;其中&#xff0c;对于每个 i&#xff08;0 < i < nums.length&#xff09;&#xff0c;都满足 ans[i] nums[nu…

程序员职业生涯的11个阶段程序人生

程序员的职业生涯是一段充满起伏的有趣经历。考虑到其陡峭的学习曲线&#xff0c;完全可以预见你将经历挫折、启蒙、骄傲自大这几个时期&#xff0c;以及穿插其间的各种心路历程。在这篇文章中让我们轻松一下&#xff0c;在作者诙谐的语言中来回顾这11个阶段&#xff1a; 阶段1…

unity python_Unity引擎内嵌python

Unity脚本using System.Collections;using System;using System.Collections.Generic;using UnityEngine;using System.Diagnostics; //需要添加这个名词空间&#xff0c;调用DataReceivedEventArgpublic class LoadPython : MonoBehaviour{string sArguments "UnityLoad…

LeetCode 1922. 统计好数字的数目(快速幂)

文章目录1. 题目2. 解题1. 题目 我们称一个数字字符串是 好数字 当它满足&#xff08;下标从 0 开始&#xff09;偶数 下标处的数字为 偶数 且 奇数 下标处的数字为 质数 &#xff08;2&#xff0c;3&#xff0c;5 或 7&#xff09;。 比方说&#xff0c;“2582” 是好数字&a…

《人性的弱点》

卡耐基-《人性的弱点》&#xff0c;讲做人要平和、真诚&#xff0c;沟通的成功在于尽量避免争辩&#xff0c;最常见的情况是在争辩中取得了胜利却失去了成功的机会。这本书适合长时间的品味&#xff0c;以至自觉地养成良好的习惯以及卓越的品格。没有什么励志书是可以让你一下子…

crc java_java实现CRC16 MODBUS校验算法

/*** 查表法计算CRC16校验**paramdata 需要计算的字节数组*/public static String getCRC3(byte[] data) {byte[] crc16_h {(byte) 0x00, (byte) 0xC1, (byte) 0x81, (byte) 0x40, (byte) 0x01, (byte) 0xC0, (byte) 0x80, (byte) 0x41, (byte) 0x01, (byte) 0xC0, (byte) 0x80…

python图片隐写_CTF 图像隐写Python脚本处理

CTF中经常会遇到很多图片的隐写题目需要使用脚本来解题&#xff0c;最常用到的就是使用python中的PIL库&#xff0c;所以如果要更好的解出图片隐写相关处理的题目&#xff0c;掌握好这个库的使用是必要的。本期就来给大家介绍下这个库的基本使用和几道图片题目的解题思路。0x00…

LeetCode 1933. 判断字符串是否可分解为值均等的子串

文章目录1. 题目2. 解题1. 题目 一个字符串的所有字符都是一样的&#xff0c;被称作等值字符串。 举例&#xff0c;"1111" 和 "33" 就是等值字符串。 相比之下&#xff0c;"123"就不是等值字符串。 规则&#xff1a;给出一个数字字符串s&…

java 资源锁定_如何在Java中创建时正确锁定资源

也许ConcurrentHashMap可以帮到你.顾名思义,它支持并发修改.要只创建一个新元素,您可以执行以下操作&#xff1a;private Map map new ConcurrentHashMap<>();private final Object lock new Object();public Thing getById(String id) {Thing t map.get(id);if (t n…

pythonb超分辨成像_深度原理与框架-图像超分辨重构-tensorlayer

图像超分辨重构的原理&#xff0c;输入一张像素点少&#xff0c;像素较低的图像&#xff0c; 输出一张像素点多&#xff0c;像素较高的图像而在作者的文章中&#xff0c;作者使用downsample_up, 使用imresize(img, []) 将图像的像素从原理的384&#xff0c;384降低到96&#xf…

LeetCode 1921. 消灭怪物的最大数量(排序)

文章目录1. 题目2. 解题1. 题目 你正在玩一款电子游戏&#xff0c;在游戏中你需要保护城市免受怪物侵袭。 给你一个 下标从 0 开始 且长度为 n 的整数数组 dist &#xff0c;其中 dist[i] 是第 i 个怪物与城市的 初始距离&#xff08;单位&#xff1a;米&#xff09;。 怪物以…

软件测试课程学习总结

一、知识结构 介绍&#xff1a; 1.Definition of Software testing: Software testing is any activity aimed at evaluating an attribute or capability of a program or system and determining that it meets its required results Software Testing is an empirical&#…

java listutils_Java的list自定义工具类ListUtils

/*** 将list中map的key为ID的值作为KEY在套一层*/public static Map> keyToID(List> datalist) {Map> res new HashMap>();for (Map map : datalist) {String id map.get("ID");res.put(id, map);}return res;}/*** 移除List中所有Map的某个元素** par…

vue 高德地图 不同区域显示不同颜色_高德百度哪家强?苹果Carplay第三方分屏功能评测...

几天前&#xff0c;苹果公司正式更新了iOS13.4版本。一个小版本系统更新&#xff0c;却让车主群热闹了起来。在这个版本中&#xff0c;苹果正式开放了Carplay分屏显示模式下对第三方地图的支持。车主们的热情&#xff0c;化为高德地图和百度地图微博下网友的催更。不过&#xf…

python 全局变量、局部变量

from 《流畅的python》 def f1(a):print(a)print(b) f1(3)# NameError: name b is not defineddef f1(a):print(a)print(b)b 5 # 全局变量 f1(3) # 输出正常python编译时&#xff0c;判断 b 是局部变量&#xff0c;因为在函数中给他赋值了 当打印 b 时&#xff0c;发现 b 没…

BZOJ 1452 [JSOI2009] Count

这道题好像有点简单的样子... absi找题目好厉害啊...确实是一道比较裸的2dBIT啊. 水掉吧. 附:2dBIT怎么做: 2dBIT就是BIT套BIT啦. 所以修改loop(xlowbit(x)){loop(ylowbit(y)){}} 查询loop(x-lowbit(x)){loop(y-lowbit(y)){}} 然后查询区间当然是用容斥... 假设查询(x11,y11)(x…

同花顺如何切换k线_30分钟线可分析出庄家的意图:教你如何用30分钟K线选股做超短线...

30分钟做超短线的好处1、兼具超短线和短线的优点&#xff0c;是联系超短周期和短周期的有利武器。2、30分钟线把一天分成8个部分&#xff0c;正好是一个神奇数字&#xff0c;自然界中很多神奇的规律不可不信&#xff01;3、对于庄家控盘的股票&#xff0c;30分钟线可分析出庄家…

LeetCode 1925. 统计平方和三元组的数目

文章目录1. 题目2. 解题1. 题目 一个 平方和三元组 (a,b,c) 指的是满足 a2b2c2a^2 b^2 c^2a2b2c2 的 整数 三元组 a&#xff0c;b 和 c 。 给你一个整数 n &#xff0c;请你返回满足 1 < a, b, c < n 的 平方和三元组 的数目。 示例 1&#xff1a; 输入&#xff1a;…

JAVA软件工程师应该具备的技能有哪些?

前言&#xff1a;有朋友问我&#xff1a;学历和能力哪个重要&#xff1f;我个人觉得能力大于学历&#xff0c;没有能力哪来的学历&#xff0c;学历只是证明能力的一方面。为此在能力方面畅谈java软件工程师必备的能力。作为一名合格的java工程师&#xff0c;不仅需要学历&#…

石头剪刀布程序流程图_机器学习终章:剪刀石头布猜拳机器人

8.1. 简介TensorFlow对于我们普通人来说一直是高冷的存在&#xff0c;都是大神们的狂欢。喵家最近将TensorFlow移植到喵家编程软件Kittenblock中&#xff0c;希望能让人工智能&#xff0c;机器学习尽快落地。能让普通用户也能用上这个好用的机器学习框架&#xff0c;解决生活中…