sgolay 函数在 Matlab R2013a 中如何工作？

Question

我对 Matlab R2013a 中的 sgolay 函数有疑问。我的数据库有 165 个光谱和 2884 个变量，我想取它们的一阶和二阶导数。如何将输入 K 和 F 定义为 sgolay？

下面是一个例子：

sgolay 用于平滑噪声正弦曲线，并将得到的一阶和二阶导数与使用diff 计算的一阶和二阶导数进行比较。请注意使用 diff 如何放大噪音并产生无用的结果。

K = 4;                 % Order of polynomial fit
F = 21;                % Window length
[b,g] = sgolay(K,F);   % Calculate S-G coefficients

dx = .2;
xLim = 200;
x = 0:dx:xLim-1;

y = 5*sin(0.4*pi*x)+randn(size(x));  % Sinusoid with noise

HalfWin  = ((F+1)/2) -1;
for n = (F+1)/2:996-(F+1)/2,
  % Zero-th derivative (smoothing only)
  SG0(n) =   dot(g(:,1), y(n - HalfWin: n + HalfWin));

  % 1st differential
  SG1(n) =   dot(g(:,2), y(n - HalfWin: n + HalfWin));

  % 2nd differential
  SG2(n) = 2*dot(g(:,3)', y(n - HalfWin: n + HalfWin))';
end

SG1 = SG1/dx;         % Turn differential into derivative
SG2 = SG2/(dx*dx);    % and into 2nd derivative

% Scale the "diff" results
DiffD1 = (diff(y(1:length(SG0)+1)))/ dx;    
DiffD2 = (diff(diff(y(1:length(SG0)+2)))) / (dx*dx);

subplot(3,1,1);
plot([y(1:length(SG0))', SG0'])
legend('Noisy Sinusoid','S-G Smoothed sinusoid')

subplot(3, 1, 2);
plot([DiffD1',SG1'])
legend('Diff-generated 1st-derivative', 'S-G Smoothed 1st-derivative')

subplot(3, 1, 3);
plot([DiffD2',SG2'])
legend('Diff-generated 2nd-derivative', 'S-G Smoothed 2nd-derivative')

Answer 1

在一个固有的嘈杂过程中取导数。因此，如果您的数据中已经存在一些噪声，实际上，当您采用高阶导数时，它会被放大。 Savitzky-Golay 是将平滑和微分结合到一个操作中的一种非常有用的方法。这是一种通用方法，它计算任意顺序的导数。不过，也有取舍。对于具有特定结构的数据，还有其他特殊方法。

就您的申请而言，我没有任何具体的答案。很大程度上取决于数据的性质（采样率、噪声比等）。如果你使用过多的平滑，你会弄脏你的数据或产生锯齿。如果您使用高阶多项式系数K 过度拟合数据，情况也是如此。在您的演示代码中，您还应该绘制sin 函数的解析导数。然后使用不同数量的输入噪声和平滑滤波器。如果您可以近似真实数据的各个方面，那么具有已知确切答案的此类工具可能会有所帮助。在实践中，我尝试使用尽可能少的平滑来生成不太嘈杂的导数。这通常意味着三阶多项式 (K = 3) 和窗口大小 F 尽可能小。

所以是的，很多人建议您use your eyes 调整这些参数。然而，最近也有一些关于自动选择系数的研究：On the Selection of Optimum Savitzky-Golay Filters (2013)。 Savitzky-Golay 还有其他替代方案，例如基于regularization 的this paper，但您可能需要自己在Matlab 中实现它们。

顺便说一句，不久前我写了一个小替换 sgolay。和你一样，我只需要第二个输出，微分过滤器，G，所以这就是它计算的全部。这个函数也更快（大约 2-4 倍）：

function G=sgolayfilt(k,f)
%SGOLAYFILT  Savitzky-Golay differentiation filters
s = vander(0.5*(1-f):0.5*(f-1));
S = s(:,f:-1:f-k);
[~,R] = qr(S,0);
G = S/R/R';

带有输入验证的此函数的完整版本是available on my GitHub。