Savitzky-Golay 滤波器详解及C/matlab语言程序设计
1.Savitzky-Golay 滤波器
Savitzky-Golay滤波器(通常简称为S-G滤波器)最初由Savitzky和Golay于1964年提出,发表于Analytical Chemistry 杂志。之后被广泛地运用于数据流平滑除噪,是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。这种滤波器最大的特点在于在滤除噪声的同时可以确保信号的形状、宽度不变。使用平滑滤波器对信号滤波时,实际上是拟合了信号中的低频成分,而将高频成分平滑出去了。如果噪声在高频端,那么滤波的结果就是去除了噪声,反之,若噪声在低频段,那么滤波的结果就是留下了噪声。最成功的应用在于心电图的基线漂移现象:基线漂移现象是有很低的频率信号引起的(体位的移动、呼吸作用),这样的低频信号有时我们也把它称为信号中的趋势项。可以利用平滑滤波器拟合该低频信号,再将其从原始的心电信号中减去,就可以抑制趋势项。总之,平滑滤波是光谱分析中常用的预处理方法之一。用Savitzky.Golay方法进行平滑滤波,可以提高光谱的平滑性,并降低噪音的干扰。S-G平滑滤波的效果,随着选取窗宽不同而不同,可以满足多种不同场合的需求。
2.信号的最小二乘平滑
信号的最小二乘平滑的基本思想可以通过下图来说明:一列数据x[n] 在图中用实心的圆点表示。考虑一组以n=0为中心的2M+1个数据,可以用如下的多项式来拟合:所以只需要获得拟合多项式的常数项。
Savitzky和Golay 发现计算a0相当于对原始数据进行一次FIR 滤波。也就是说可以利用卷积运算来实现。
也可以说是对输入数据进行了加权平均。图1中x 表示的就是加权系数。下面就来说说如何获得a0。
由基本的微积分知识可知,若要ε最小,ε对各个参数的偏导数都应为0,也就是:
因此,Savitzky-Golay卷积平滑算法是移动平滑算法的改进:
3.matlab仿真测试
noise = wgn(1,101,2); for n=1:1:101s(1,n) = 3*sin(0.2*n)+3*sin(0.4*n); end x = s + noise; %滤波器 % y = sgolayfilt(x,3,7); h = [-2, 3, 6, 7, 6, 3, -2 ]/21; y = zeros(1,101); y(1,1:3) = s(1,1:3); y(1,99:101) = s(1,99:101); for k=4:1:98y(1,k) = x(1,k-3)*h(1)+x(1,k-2)*h(2)+x(1,k-1)*h(3)+x(1,k)*h(4)+...x(1,k+1)*h(5)+x(1,k+2)*h(6)+x(1,k+3)*h(7); end plot(x,'b'); hold on plot(y,'r'); hold on plot(s,'g'); hold off
MATLAB自带的sgolayfilt滤波器的处理结果:依据算法思想进行卷积滤波处理结果:
4.算法的C语言实现
要用C 语言从头写起这个代码会非常的长,光一个 ployfit 函数实现起来就很麻烦。程序中直接调用了 GSL 的相关函数。先说说 polyfit 函数。多项式函数拟合可以转换成超定线性代数方程组的最小二乘解的问题。标准解法是 SVD 分解。不过GSL库中提供了gsl_multifit_linear函数进行线性模型的拟合,可以直接使用。下面给出代码。很简单,就不多做解释了。这个例子也可以作为 gsl_multifit_linear 函数用法的一个小例子。/** * 多项式拟合函数 * 根据 x y 拟合出一个N次多项式。返回多项式的系数。 */ gsl_vector* PolyFit(gsl_vector *x, gsl_vector *y, int N) { int i, j; int M = GSL_MIN(x->size, y->size); double chisq, xi; gsl_matrix *XX = gsl_matrix_alloc(M, N + 1); gsl_vector *c = gsl_vector_alloc(N + 1); gsl_matrix *cov = gsl_matrix_alloc(N + 1, N + 1); for(i = 0; i < M; i++) { gsl_matrix_set(XX, i, 0, 1.0); xi = gsl_vector_get(x, i); for(j = 1; j <= N; j++) { gsl_matrix_set(XX, i, j, pow(xi, j)); } } gsl_multifit_linear_workspace *workspace = gsl_multifit_linear_alloc(M, N + 1); gsl_multifit_linear(XX, y, c, cov, &chisq, workspace); gsl_matrix_free(XX); gsl_matrix_free(cov); gsl_multifit_linear_free(workspace); return c; }
有了 PolyFit 函数就可以计算 SG 滤波器的系数了。Matlab 中的 polyval 函数 可以用 gsl_poly_eval 来实现。
/** * 计算 Savitzky-Golay 滤波器系数 * SG 滤波器的阶数为 2M+1,多项式的最高次数为 N */ gsl_vector* SG_FilterCreate(int M, int N /* Poly Order */) { int i; gsl_vector *x = gsl_vector_alloc(2 * M + 1); gsl_vector *y = gsl_vector_alloc(2 * M + 1); gsl_vector_set_zero(y); gsl_vector_set(y, M, 1); for(i= -M; i <= M; i++) { gsl_vector_set(x, i + M, i); } gsl_vector *c = PolyFit(x, y, N); gsl_vector_free(x); gsl_vector_free(y); gsl_vector *fir = gsl_vector_alloc(2 * M + 1); for(i = -M; i <= M; i++) { gsl_vector_set(fir, i + M, gsl_poly_eval (c->data, N + 1, i)); } gsl_vector_free(c); return fir; }
5.平滑滤波的尴尬
当然,这是针对所有FIR滤波器而言,都会有以下几个尴尬,主要在模板值大小上:1.N足够大,就可以获得足够小的NRR(Noice Reduce Rate)2.N过大会使滤波器具有过大的延迟3.N过大会使其主瓣的单边的带宽大大降低,这就有可能在滤波时使有用的信号(或某一个频率组份)也受到损失;因此,在平滑中,N不宜取得过大。致谢:本文的大部分工作都是在CSDN liyuanbhu博主研究基础之上进行的,感谢他的辛苦劳动。
Savitzky-Golay 滤波器详解及C/matlab语言程序设计相关推荐
- gammatone 滤波器详解及其MATLAB代码实现
一.GammaTone 滤波器详解 定义: 外界语音信号进入耳蜗的基底膜后,将依据频率进行分解并产生行波震动,从而刺激听觉感受细胞[1].GammaTone 滤波器是一组用来模拟耳蜗频率分解特点的滤波 ...
- 程序阅读_全面详解LTE:MATLAB建模仿真与实现_自学笔记(1)调制与编码_程序阅读
程序阅读_全面详解LTE:MATLAB建模仿真与实现_自学笔记(1)调制与编码_程序阅读 在粗浅地掌握了LTE知识后,从今天开始对<全面详解LTE:MATLAB建模仿真与实现>一书的学习. ...
- RC滤波器与LC滤波器详解(RC滤波器与LC滤波器区别,RC滤波器和LC滤波器工作原理和经典设计)
RC滤波器与LC滤波器详解(RC滤波器与LC滤波器区别,RC滤波器和LC滤波器工作原理和经典设计) [导读] LC滤波器应用的频率范围为1kHz-1.5GHz.由于受限于其中电感的Q值,频率响应的截至 ...
- C语言位操作详解-4.2.C语言专题第二部分-朱有鹏-专题视频课程
C语言位操作详解-4.2.C语言专题第二部分-11753人已学习 课程介绍 本课程的目标是让大家彻底掌握C语言中位运算符及常规技巧,这些知识在嵌入式代码中用的很多,尤其是设置寄存器以操 ...
- 【一天一门编程语言】Matlab 语言程序设计极简教程
Matlab 语言程序设计极简教程 用 markdown 格式输出答案. 不少于3000字.细分到2级目录. 文章目录 Matlab 语言程序设计极简教程 简介 Matlab 工作空间 Matlab ...
- ise的时钟ip核_Vivado CIC IP核滤波器详解(一)
引言:从本文开始,我们详细介绍Xilinx CIC IP核滤波器相关知识,包括CIC IP核提供的特性.IP核接口描述以及IP核设计指导等相关内容. 1.概述 级联积分器梳状(CIC)滤波器,也称为H ...
- 主成分分析,PCA(Principle Component Analysis)详解 ,结合matlab
https://my.oschina.net/gujianhan/blog/225241 一.简介 PCA(Principal Components Analysis)即主成分分析,是图像处理中经常用 ...
- 基于YOLOv5的目标检测系统详解(附MATLAB GUI版代码)
摘要:本文重点介绍了基于YOLOv5目标检测系统的MATLAB实现,用于智能检测物体种类并记录和保存结果,对各种物体检测结果可视化,提高目标识别的便捷性和准确性.本文详细阐述了目标检测系统的原理,并给 ...
- 大二c语言期末考试题库及详解答案,大学C语言期末考试练习题(带详解答案)...
<大学C语言期末考试练习题(带详解答案)>由会员分享,可在线阅读,更多相关<大学C语言期末考试练习题(带详解答案)(55页珍藏版)>请在金锄头文库上搜索. 1.一. 单项选择题 ...
最新文章
- Cisco ××× 完全配置指南-连载-IPSec
- JavaScript 技术篇-chrome利用ClipboardEvent写入剪切板,chrome写入剪切板不成功解决方法
- Keras处理TIFF图像
- java yii_构建 Java 应用程序
- 赋值给集合_ArrayList集合源码
- antd中tooltip换行_ant design 中实现表格头部可删除和添加
- stm32驱动LTC6912程控放大器程序,PGA可编程增益放大器,可调增益运放电路
- 前端学java还是python_零基础小白 选择学习 java、php、前端 还是 python?
- oracle在cmd中启动数据库实例
- 计算机组成原理:系统总线IO设备
- 14届数独-真题标准数独-Day 9-20220124
- bzoj21012101: [Usaco2010 Dec]Treasure Chest 藏宝箱(滚动数组优化dp)
- [PHP]md5加密
- LEF和GDS匹配问题
- 基于MATLAB 的运动模糊图像复原
- windows XP全公略
- 新浪小编为你介绍留声机的古与今
- Centos7文本处理工具
- 2021年抖音0粉丝无货源直播带货最新技巧
- python写采集程序_Python写的简易采集爬虫(蜘蛛)