Savitzky-Golay 滤波器的核心代码主要集中在计算投影矩阵B
并使用这个矩阵对输入信号进行滤波。这部分核心代码包括计算B
矩阵、处理边界效应和进行实际滤波操作。以下是对核心代码的一点解释:
① 计算 Savitzky-Golay 投影矩阵B
B = sgolay(order, framelen, weights); % 计算 Savitzky-Golay 投影矩阵
Savitzky-Golay 投影矩阵B
是通过多项式拟合计算出来的。这个矩阵用于对输入信号进行平滑处理。函数 sgolay
负责生成这个矩阵。
② 处理数据的维度
if isempty(dim)
[x, nshifts] = shiftdim(x); % 沿第一个非单维度工作
else
perm = [dim, 1:dim-1, dim+1:ndims(x)];
x = permute(x, perm); % 将 DIM 置于第一维度
end
if size(x, 1) < framelen, error(message('signal:sgolayfilt:InvalidDimensionsTooSmall')), end
这段代码确保信号x
的长度至少为 framelen
,并将数据沿指定维度进行处理。
③ 计算滤波结果
前端过渡
ybegin = B(end:-1:(framelen-1)/2+2,:) * x(framelen:-1:1,:);
这部分代码计算信号的前端过渡部分。通过倒序取出 framelen
个数据点并与矩阵B
相乘,计算出平滑后的前端部分。
稳态输出
ycenter = filter(B((framelen-1)/2+1,:), 1, x);
这部分代码计算信号的稳态输出部分。通过应用 filter
函数进行卷积操作,实现对信号中心部分的平滑处理。
后端过渡
yend = B((framelen-1)/2:-1:1,:) * x(end:-1:end-(framelen-1),:);
这部分代码计算信号的后端过渡部分。通过倒序取出信号末尾的 framelen
个数据点并与矩阵B
相乘,计算出平滑后的后端部分。
合并结果
y = [ybegin; ycenter(framelen:end,:); yend];
将前端过渡部分、稳态输出部分和后端过渡部分合并,得到完整的平滑信号。
④ 恢复原始数据的形状
if isempty(dim)
y = shiftdim(y, -nshifts); % 恢复原始维度
else
y = ipermute(y, perm); % 恢复原始维度顺序
end
将平滑后的信号恢复到与输入信号相同的维度顺序。