Matlab中Savitzky-Golay filtering（最小二乘平滑滤波）函数sgolayfilt的使用方法

语法规则

y = sgolayfilt(x,order,framelen)
y = sgolayfilt(x,order,framelen,weights)
y = sgolayfilt(x,order,framelen,weights,dim)

语法描述

y = sgolayfilt(x,order,framelen)：对数据向量x使用Savitzky-Golay FIR平滑滤波器。如果x是一个矩阵，则sgolayfilt对每一列进行操作。多项式阶次order必须小于框长度framelen，因此framelen必须是奇数。如果order=framelen-1，则滤波器不进行平滑。

y = sgolayfilt(x,order,framelen,weights)：指定了一个权重向量weights，它为正实数，用于最小二乘估计的最小化。如果该值没有指定，或定义为空[]，它默认为单位矩阵。

y = sgolayfilt(x,order,framelen,weights,dim)：指定了维数dim。如果dim没有被指定，则sgolayfilt对第一个不为1的维进行操作；指定为1则对列向量进行操作，指定2则对行向量进行操作。

举例

稳态和瞬态Savitzky-Golay滤波器

生成一个随机信号并使用sgolayfilt平滑。指定多项式阶数为3，框长度为11.画出原始的和平滑后的信号。

order = 3;
framelen = 11;

lx = 34;
x = randn(lx,1);

sgf = sgolayfilt(x,order,framelen);

plot(x,':')
hold on
plot(sgf,'.-')
legend('signal','sgolay')

sgolayfilt函数通过对sgolay的输出B的中间行进行卷积来执行大部分滤波。结果是过滤的信号的稳态部分。生成并画出这个部分。

m = (framelen-1)/2;

B = sgolay(order,framelen);

steady = conv(x,B(m+1,:),'same');

plot(steady)
legend('signal','sgolay','steady')

靠近信号边缘的样本不能放置在对称窗口的中心，并且必须以不同的方式进行处理。

为了确定启动瞬态，矩阵将B的第一个(framelen-1)/2行乘以信号的第一个framelen长度的采样。

ybeg = B(1:m,:)*x(1:framelen);

为了确定终端瞬态，矩阵将B的最后的(framelen-1)/2行乘以信号的最后的framelen长度的采样。

yend = B(framelen-m+1:framelen,:)*x(lx-framelen+1:lx);

连接瞬态和稳态部分以产生完整的信号。

cmplt = steady;
cmplt(1:m) = ybeg;
cmplt(lx-m+1:lx) = yend;

plot(cmplt)
legend('signal','sgolay','steady','complete')
hold off

添加权重会破坏B的对称性，并且需要额外步骤来提供合适的解决方法。

对声音信号的Savitzky-Golay滤波

用三阶Savitzky-Golay平滑mtlb信号，数据框长度为41.

load mtlb
smtlb = sgolayfilt(mtlb,3,41);

subplot(2,1,1)
plot(1:2000, mtlb(1:2000))
axis([0 2000 -4 4])
title('mtlb')
grid

subplot(2,1,2)
plot(1:2000,smtlb(1:2000))
axis([0 2000 -4 4])
title('smtlb')
grid