Savitzky-Golay平滑滤波器的最小二乘拟合原理综述

DOI：10．3969／j．issn．1001—3824．201 1．01．017 Savitzky．Golay平滑滤波器的最小二乘拟合原理综述 蔡天净 ，唐 瀚 (1．重庆邮电大学 通信与信息 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 学院，重庆400065；2．中南大学 数学科学与计算技术学院，湖南 长沙 410075) 摘 要：介绍了Savitzky—Golay滤波器的推导方法—— 多项式的最小二乘拟合法及其推导过程，以及如何由Savitzky 和 Golay提 出的多项式卷积计算方法进行最小二乘拟合计算，对 Savitzky—Golay滤波器的二维算法进行了简单介 绍，对其一维和二维的MATLAB代码进行了分析处理，并将 Savitzky—Golay滤波器同其他低通滤波器进行了简单比 较，简要 说明 关于失联党员情况说明岗位说明总经理岗位说明书会计岗位说明书行政主管岗位说明书 了其优势以及一些应用方向。 关键词：Savitzky—Golay滤波器；数据平滑；多项式最小二乘拟合；卷积 0 引 言 1 Savitzky-Golay滤波器算法及推导 Savitzky—Golay滤波器是一种特殊的低通滤波 器，又称 Savitzy—Golay平滑器。低通滤波器的明显 用途是平滑噪声数据，噪声是用来描述所观察现象 提取信息中附加的不易区别的任意错误，而数据平 滑能消除所有带有较大误差障碍的数据点，或者从 图形中作出初步而又粗糙的简单参数估算。 Savitzky—Golay滤波器最初由 Savitzky A和 G0． 1ay M于 1964年提出，被广泛地运用于数据流平滑 除噪，是一种在时域内基于多项式，通过移动窗口 利用最小二乘法进行最佳拟合的方法。这是一种 直接处理来 自时间域内数据平滑问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 的方法，而不 是像通常的滤波器那样先在频域中定义特性后再 转换到时间域。通过这种方法，计算机的唯一功能 就是充当一个平滑噪声起伏的滤波器并尽量保证 原始数据的不失真。在这个过程中，计算机只需运 行相对小型的程序，减少了对电脑内存和数据处理 能力的要求，因此这种方法相对来说更加简单、快 速，而且相对于其他类似的平均方法而言，这种方 法更能保 留相对极大值、极小值和宽度等分布 特性 。 收稿日期：2010-03—19 1．1 最小二乘法拟合与Savitzky-Golay算法的导出 在平面坐标系中，用一条曲线来拟合一组数 据，不妨假设这条曲线为 Y=00+0l +I22X +03 +。 ，当每一个点的横坐标代入这个曲线方程后， 所得的值与该点的纵坐标之差的平方之和最小时， 这条曲线的拟合度最高，从而可以确定所有的系数 n (i=0，1，2，3，4)。 如图 1所示，有 9个数据点被在左边的括号包 住，如果这些点均在一条曲线附近，则能近似地被 所示的能描述这条曲线的方程所表示，然后就可以 根据数值计算方法用具体 的过程把所有 的系数 (i)确定出来。将中间点的横坐标代回到方程中， 得到图中圆圈所表示的点，用这个点代替原来 的 点。这个过程中，获得的该点值的大小是在最小二 乘法和这样一组观察数据点基础上的最佳拟合。这 一 过程可以对有9个数据点的每一组数重复进行，每 进行一次就去掉最左边的一个点然后加上最右边的 一 个点，一直到最右边的区域都确定出来。通常来 说，每一组9个点确定的曲线系数都不同 ]。 现讨论由Savtzky和 Golay导出的方法。设一组 数据为 (i)，i的取值为2m+1个连续的整值，即i =一m，⋯，0，⋯，m。现 构 造 一 个 n阶 多 项 式 (n≤2m+1)来拟合这一组数据 一 63 — 图1 9点表示的滑动平均多项式 = ∑b．ki =b加+boli+bn2i +⋯6 (1) = 0 设点的误差平方和是 E=∑If,一 ( )] = ∑[∑b．ki 一 ( )] (2) 为使 E最小，可令 E对各系数的导数为0，即 _0， ，l，2，⋯，n (3) 得 蓑=彘{ m【 n 2 ．Z[ 6 一 ( = m =O 一 0，r=0，1，⋯，n (4) 即 ∑b ∑i =∑ (i)i (5) 令 S =∑i¨ 1 F，=∑b,e,S (6) l 0 J 给定需要拟合的单边点数m，多项式的阶次n及待拟 合的数据 (一m)，⋯，x(o)，⋯， (m)，则可求出F，， 将s 代入式F，=∑b．S 中，系数b加，b ⋯，b 就可求出，因此多项式 可以确定。 由于r+k是奇数时有 5 =∑i¨ =0 (7) — — 64 —— DIGITAL COMMUNICATION／201 1．1 所以，只有当r+k为偶数时，才能存在 =∑b S ，r=0，1，2⋯凡的表达形式 引。 总的来说，当n和S均为偶数或者 n和s均为奇 数时有 b =b( ) (8) 在实际应用中，往往并不需要把系数 b ，b ⋯b 全部求出，分析可得 I 0=b棚=a 孔。 。 以 l =n!b = dt I ：0 即 S!b =a 。 Savitzky和 Golay使用了一种简单的多项式卷 积方法，并制定了卷积系数表 j，通过使用卷积系 数表中的卷积系数来计算出系数 b棚，b ⋯b 的大 小，方便而且迅速。 例如，用 1个四次多项式f／=b加+b li+bn2i + 6 i +bn4i 来拟合9个点，先使用表求 b加。如表 1 所示，找到表 1中9个点对应的一列系数为15，一55， 30，135，179，135，30，一55，15，对应的分母为429。 表 1 Savitzky．Golay卷积系数表 110 — 198 — 16O l1O 5 — 3O 75 13l 75 — 3O 5 设多项式拟合的9个数据值的大小分别为 Y Y一3，Y一2，Y一1，Yo，Yl，Y2，Y3，Y4 则有 5 5 4 5 3 5 2 5 4 3 2 l ∞ ∞ 枷∞啪 ∞圳 0 8 1 ∞ " ∞ ∞ m m 3 6 6 6 3 l 一 一 l ^ l 一0。2 3 4 5 6 — l5y_4—55y 一3+30y一2+135y—l+179yo+135y1+30y2—5 ! “40— 429 640= a4o (10) 参考以上方法利用卷积系数表 可分别求得 Ⅱ ， 042，043，口44的大／J、。 又 昔， 署 ⋯ 6 =著 =并 则可求得 b ，b以，b ，b 的大小。 使用卷积系数表时，先根据所求系数 b 找到相 应的表，每一个表中都有一个A 形式的标识，其中i 表示拟合多项式权数， 表示导数的阶数，找到相应 的表后，则根据拟合需拟合的点的个数查表确定其 卷积系数以及分母大小 。 1．2 一维算法实现 Savitzky-Golay方法存在一个主要的缺点，即这 种计算方法截断了左右两端各 m个点，这对于在两 端为零值的大量光谱来说影响不大，但是对于在两 端有极限值的数据组，该算法是不适用的。 事实上如果进行一个规范的最小二乘运算，产 生出相应的系数，用来算前 m个点中任何一个点的 拟合都 比计算 中间的那个点简单。这种方法被 Khan采用并使用矩阵方程进行全最小二乘拟合，但 是这样一来就失去了 Savitzky—Golay方法的最大优 势——即计算简便、迅速 ]。 仍然构造一个 乃阶多项式 )=∑ap (12) 令第 个数的横坐标为 (0≤i≤s一1)，同时令 。=0以简化参考，并且令它的值的大小为 g ，在 对 。 的求解过程中从简化了最小二乘的程度上优 化了 )对g 的拟合 s = )一gi = ∑ = 0，0≤ ≤n (13) 用如下的方式进行伪转置 -厂( )=g B = 1 1 ⋯ ⋯ ； ； ； ； 1 i ⋯ ： ⋯ j i 1 ⋯ ： ⋯ a= (0o，⋯。 ) Ba = g j 曰 Ba = B g = (B B)一 B g a cf 0 =塞 i ⋯ 一台(一c)! ： s(n BTg S=( “c ， 一 ⋯，⋯ n! ⋯ 、 = J 令 =S(B B) B ，这样， 便是一个矢量， 这个矢量与样本矢量g点乘的结果就是 处的 c阶 导数 。 接下来应用 Savitzky—Golay滤波器处理一个数 据长度为 31的信号，程序代码及结果如图2所示。 ． f 、 一 。1_】 O 200 400 6o0 8oo 1 oo0 1 2o0 l 4001 6ool 80o 2 000 (a)原始数据 一 ! ： 一 1 ⋯ ： ： l l I 。 0 200 400 600 8【x】1 o00 l 2oo l 4o0l 6()0l 800 2 000 (b)平滑后数据 图2 利用 Savitzky-Golay滤波器处理信号对比 1．3 二维算法实现 Savitzky-Golay滤波器对二维图像的处理同样 适用，此时可以用一个二维多项式来拟合图像的二 维区域，对图像进行平滑处理或者数字分化 。 构造一个 n阶多项式 n n—D ， )=∑∑0 Y (14) 则关于 的 m阶导数为 一 65 — 关于Y的 m阶导数为 o y mf = ny nx i 假定要平滑数据或分化一个7×7的图像块，图 像块数据 d(i)和它在坐标字中的位置如表2所示。 表 2 图像块数据 d(i)和它在坐标系中的位置 表2中：d( )是像素值，列向量 d代表表 2所 有的数据，即 d=(d(0)，d(1)，⋯d(48)) 假定要用 1个二维三次多项式来拟合这些数 据，多项式可设为 d(i) ，Y )= aoo + aloXi 十 口0l 十 n2o + Ⅱl1 Y +ao2y + 口30 + 口2l y +Ⅱ12 Y + ao3y + 口40 + 口31 y + a22 2 y2 + 口13 y + no4y4 (15) 式(15)中：x'f的系数是a ；( ，)， )是 d( )的像素 坐标。 建立矩阵方程 xa ： d f 16) 0 Yo XoYo 1 Yl X1Yl 札8 扎8 x48Y48 4 O 4 I ： ● 4 x48 a是多项式系数矢量，a=(a00，alo，aol，a2o，all，ao2， 030，⋯，。03，040，⋯ ，oo4)T 则有 a ： fxTx)～ d 接下来应用 Savitzky-Golay滤波器处理一个二维图 像信号 #001 function h： Sg _ 2d(x，Y，FIX，ny，d，flag_coupling) — — 66 —— DIGITAL COMMUNICATION／201 1．1 #002 if nargin<6 #0 03 flag — coupling=0； #004 end #005 if nargin<5 #006 d=0； #007 end #008 if nargin<4 #009 ny=1； }}010 end #0 11 if nargin<3 #012 nx=1； #0 13 end #014[X，Y]=mesh d(X，y)； #015[1y，lx]=size(X)； #016 if nx>Ix一1 Il ny>ly一1 #0 17 error (polynomial order too large!’； 加 18 end #o19 if d>nx 革}020 error (tlifferentiation order too large!，； 加 2l end #022 x=x(：)； #023 Y=Y(：)； #024 A=l j； #025 if flag _ coupling 峁026 for kx=nx：一1：0 #027 for ky=ny：一1：0 #028 A=[A X． kx． Y． ky]； !o29 end }1030 end 加31 else 加32 for k=nx：一1：0 #033 A=[A X． k]； }{034 end #0 35 for k=ny：一1：1 }}o36 A=[A Y． k]； }}037 end 嵌)38 end #039 h=inv(A }A) A ； #040 if flag _ coupling {f041 h=h((BX+1一d) (ny十1)，：)； #042 else #043 h=h(nx+1一d，：)； 抻344 end #045 h=reshape(h，ly，lx)； #046 h=prod(1：d)$h； 输入 ∑ ～ ， ， ； ， ， ， ‰ -芸 儿 i ． ． 1 中 ， 、 其 ： Sg 2d(一3：3，一4：4，2，3，2，1) 即可得结果 anS = 一 0．010 8 0．000 0 0．006 5 0．008 7 0．006 5 n 000 0 —0．010 8 007 2 0．000 0 —0．004 3 一n 0058 — 0．004 3 0．000 0 0．007 2 0．020 1 0．000 0 —0．012 1 — 0．016 1 —0．012 1 0．000 0 0．020 1 0．027 8 0．000 0 — 0．016 7 —0．022 3 —0．016 7 0．000 0 0．027 8 0．030 4 0．000 0 —0．018 2 一Q 024 3 —0．018 2 n 000 0 0．030 4 0．027 8 0．000 0 —0．016 7 —0．022 3 一n 016 7 0．000 0 0．027 8 0．020 1 0．000 0 一n 012 1 一n 016 1 一n012 1 0．000 0 0．020 1 0．007 2 0．000 0 —0．004 3 —0．005 8 — 0．004 3 0．000 0 0．007 2 —0．010 8 0．000 0 0．OO6 5 0．008 7 0．006 5 0．000 0 —0．010 8。 2 分析讨论 2．1 Savitzky．Golay滤波器与其它低通滤波器滤波 特性的对比 2．1．1 Savitzky—Golay滤波器与 Butterworth滤波器 Buttterworth滤波器被广泛使用并以其最大平坦的 特点而闻名，它是大小能根据方程计算的低通类型。 ． 1 lH(jW)l ——■— 1+f ＼ 、 ／ = 2 (17) 级数 凡越高，Butterworth滤波器特性越接近于 理想低通滤波器。 如图3所示。一至五阶(图中编号 1—5)巴特沃 斯低通滤波器频率响应曲线展示了一阶至五阶巴 特沃斯低通滤波器的频率响应曲线。可见阶数 n越 高，其幅频特性越好，低频检测信号保真度越高 J。 Savitzky-Golay滤波器与 Butterworth滤波器 的 比较如图4所示。 0 - 20 璺 一40 甜 一 一 60 - 80 一 l00 ＼ l＼＼＼ f ＼ ＼ll＼ 0．0l 0 1 1 10 1o0 频"~／rad sl 图3 一至五阶巴特沃斯低通滤波器频率响应曲线 2．1．2 Savitzky-Golay滤波器与移动平均滤波器 移动平均滤波器是通过输入信号的一些点来 产生输出信号中数据点的方法。这是一个需要重 复的过程，因此选 M个点为移动窗口时可采用如下 算法得到数据集的平均值 。 『 ∑ > 且 <Ⅳ一 Yi { L0 其他 (18) 式(18)中：s 是输出信号；Y 是输出信号；M是移 动窗口中的点数；Ⅳ是数据点总数。 作为一种替代，输入信号点可以选择围绕输出 点对称，如果这样的话， 就必须是奇数，用方程表 示为 > 且 <Ⅳ 一 其他 (19) Savitzky—Golay与移动平均滤波器的比较如图5 所示。 (h)巴特沃斯滤波结果 (c)Savi~ky—Golay滤波结果 图4 Savitzky-Golay与 Butterworth滤波器的比较 2．2 Savitzky-Golay滤波器的优点 Savitzky—Golay滤波器有如下几大优点： 1)利用最小二乘的多项式拟合方法非常清晰 易懂，并且在计算上来说，多项式卷积的操作 比最 小二乘的计算可操作性更强； 2)滤波系数只需要在对应的卷积系数表中进 行查找，很容易获得； 一 67 一 ∑ 0 ，●●●●●J‘l，●，L = y 3)Savitzky—Golay滤波器可以有任意的长度，因 此有利于采样频率通常很低的生物学或者生物力 学的数据处理 。 (b)巴特沃斯滤波结果 (c)Savitzky—Golay滤波结果 图5 Savitzky-Golay与移动平均滤波器的比较 2．3 实际应用状况 2．3．1 修复损坏的传感器 可以使用由合适维数滤波器产生的平滑量代 替传感器中损坏的位置，以此来做一个最好的估 测，同时可以忽略损坏的元素来计算空问导数。 2．3．2 图像的金字塔模型构建和多栅极问题 用于采样图像为多尺度分析建立分辨率金字塔 已经使用广泛，但如果同时使用 Savitzky-Golay方法， 这个过程可以有效地扩展成包括限制极数的多栅极 高维问题。网格几何形态不仅仅只是矩形，它也可能 是基于三角形的，只要它们能在这个区域上形成一个 重复的图案，从而能够被平滑、分化或采样。 2．3．3 时间序列问题 可以以平面直角坐标系和时间建立一系列的 图像模型，图形中一块正方形的区域是最普通的， 除非有强制性的原因去采用另外的外貌或形状，而 且可能会用了3个时间点和一个空间的立体模型来 故居边缘探测，以及使用一个线性的或者时域模型 来呈现两幅画之间极小的变化 ’m 。 3 结 论 Savitzky—Golay滤波器是一种时域内的低通滤波 器，它的基本思路是基于多项式，通过移动窗口利用 最小二乘法对数据进行最佳拟合，并且它的最初提出 者 Savitzky和 Golay提供了一种利用多项式卷积的计 — — 68 —— DIGITAL COMMUNICATION／201 1．1 算方法来计算滤波系数，这种方法简单易懂且计算简 便、迅速，可操作性非常强。Savitzky．Golay滤波器的 最初计算方法有一定的局限性，会导致部分数据丢 失，但后人对该方法进行了多次改进和完善口 ，并 且拓展到二维甚至多维的应用 -l6_，在高光谱分 析，优化生物数据图边缘等领域的应用也有一定的 研究成果，是一种实用性很强，除噪效果非常好的 滤波器 参考文献： [1] WIKIPEDIA．Savitzky Golay smoothing filter[EB／OL]． (2010-01—15)[2010-02—15]．http：／／en．wikipedia． org／wiki／Savitzky—Golay — smoothing_ filter． [2] SAVITZKY A，GOLAY M J．Smoothing and differentia— tion of data by simplified leastsquares procedures[J]． Analytical Chemistry，1964(7)：1627—1639． [3] WILLIAM H P，SAUL A T，WILLAM H V，et a1．C++ 数值算法[M]．胡建伟，译．北京：电子工业 出版社， 2005：117—130． [4] GORRY P A．General least—squares smoothing and differ- entiation by the convolution(Savitzky—Golay)[J]．Ana— lytical Chemistry，1990(3)：570—573． [5] THORNLEY D J．Anisotropie multidimensional Savitzky Golay kernels for smoothing，differentiation and reeon— struetion[J]．Departmental Technical Reports，2006 (7)：1123—1127． [6] KRUMM J．Savitzky—Golay Filters for 2D Images．[EB／OL]． (21301-08-01)[2010-02—15]．http：／／research．microsoft．corn／ en—us／urn／people／'jeknmma／SavGol／SavGo1．htm． 『7] NACHAIYAPHUM K，SUJITJORN S，RUGMAI S．A— daptive wiener filter based numerical filter with an appli— cation to beam position monitoring[J]．Wseas transae— tions on electronics，2008(2)：40—52． [8] GUINON J L，ORTEGA E，ANTON J G，et a1．Moving Average and Savitzky—Golay Smoothing[EB／OL]．(2006— 1O一19)[2010—10—20]．http：∥ ．omicronrseh．corn／ Art5 一 SG％ 20Smoothing／Using％ 20SG％ 20Smooth— ing．htm． [9] LUO J W，YING K，HE P，et a1．Properties of Savitzky一 lay digital differentiators．Digital Signal Processing[EB／ OL]．(2OO6—10—20)[2010—10-20]．http：∥WWW．es．wright． edu／— phe／Research／DigitalSignalprocessing-05．pdf． [10]BAKKALI S．Using savitzky—golay fliting method to opti— mize surface phosphate deposit disturbances[J]．Ingenie— rias，Abril—Junio，2007(7)：62-67． (下转第82页) 大的提高，算法效果理想。 4．2．2 手写测试 把算法程序封装成一个单独的函数供驱动程 序调用，结合 LCD画点函数在Pc机上交叉编译后， 动态加载到 ARM9开发板上。当触摸笔在屏上滑 过时，笔迹在液晶屏上显示出来。 图3是触摸笔在屏上书写汉字时，液晶屏上显 示效果的截图。图 3(a)是未做数据处理时字体效 果，飞点现象明显，图3(b)是经过消抖滤波算法处 理后的效果，线条更趋向于平滑，飞点基本消除。 ■ ■ (a)未做数据处理 (b)消抖滤波算法处理 图3 触摸屏数据处理效果对比图 5 结束语 结合实际的硬件平台，详细地介绍了基于 C／ 0s．II操作系统下触摸屏驱动程序的开发流程 。 同时，针对当前触摸屏数据采样中存在的飞点和抖 动问题，从软硬件 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 角度提出了切实可行的解决 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 ，增强了触摸屏的抗干扰性能，提高了触摸屏 定位识别精度。 参考文献： [1] [2] [3] [4] [5] [6] 蔡海蛟，危峻．便携式红外相机 中触摸屏原理与应用 [J]．计算机工程与设计 ，2008，29(7)：1808—1813． 史蕊，蔡浩，王振 ，等．基于S3C44BOX IxC／OS—II的触 摸屏设计[J]．电测与仪表，2007，44(494)：50-53， 张毅，王海涛．基于 $3C2410A的WinCE5．0下触摸屏 驱动的实现[J]．重庆邮电大学 学报 ：自然科学版， 2008，20(6)：742-745． 张颖超，杨乐，叶小玲，等．基于 Linux的触摸屏[J]．仪 表技术与传感器，2008(2)：71-73． 李海青，李树广．Windows CE．NET下触摸屏操作系统 优化与研究[J]．测控技术，2009，28(4)：91．94． 张毅，王海涛．基于 $3C2410A的 WinCE5．0下触摸屏 驱动的实现 [J]．重庆邮电大学学报 ：自然科学版， 2oo8，20(16)：742-745． 作者简介： 陈 勇(1963一)，男，博士，教授，主要研究方向为智能 仪器仪表 和 自动化装置 与系统，E—mail：ehenyong@equpt． edu．en． (上接第68页) [11 J ARVOR D，JONATHAN M，MEIRELLES M，et a1．Com— pafison of muhitemporal MODIS·-EVI smoothing algo·· rithms and its contribution to crop monitoring[EB／OL]． (2008-01-20)[2010—10—20]．http：／／ieeexplore．ieee． org／stamp／stamp．jsp?arnumber=04779155． [1 2]HARGITYAI S．Savitzky—Golay least—squares polynomial filters in ECG signal processing[EB／OL]．(2005-01 J01) [2010—10-20]．http：∥ieeexplore．ieee．or#xpls／abs— al1．isp?amumber=1588216． [13]ZHAO H，YANG Z W，DI L P，et a1．Crop phenology date estimation based on NDVI derived from the reeon— structed MODIS daily sul~aee reflectance data[EB／OL]． (2009-03．15)[2010—10—20]．http：∥ieeexplore．ieee． or#stamp／stamp．jsp?arnumber=05293522． [14]WAYT H J，KHAN T R．Integrated Savitzky—Golay filter from inverse taylor series approach[EB／OL J．(2007—10一 一 82 一 DIGITAL COMMUNICATION／201 1．1 15)[2010—10—20]．http：／／ieeexplore．ieee．org／stamp／ stamp．jsp?arnumber=04288597． [15]RUFFIN C，KING R L．The analysist of hyperspeetral data using savitzky—golay flitering[EB／OL]．(1999-01— 01)[2010—10-20]，http：∥www．ece．msstate．edu／～rk— ing／downloads／savitzky％ 20golay％202．pdf． [1 6]CHINRUNGRUENG C，TOONKUM P．Real—time speckle reduction and coherence enhancement of ultrasound ima— ges based on mixture of anisotropic Savitzky—Golay filters [EB／OL]．(2004—10—20)[2010—10-201．http：∥ieeex— plore．ieee．0r stamp／stamp．jsp? tp= &amumber= 1466820&userType=inst． 作者简介： 蔡天净(1990一)，女，安徽安庆人，本科，主要研究方向 为通信 工程 ，E—mail：anderenya@yahoo．corn．ca。