我有一个关于 sgolay
的问题Matlab R2013a 中的函数。我的数据库有 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')
最佳答案
在一个固有的嘈杂过程中求导。因此,如果数据中已经存在一些噪声,那么当您采用高阶导数时,噪声确实会被放大。 Savitzky-Golay 是将平滑和微分结合到一个操作中的一种非常有用的方法。这是一种通用方法,它计算任意阶的导数。不过,这也是需要权衡的。对于具有一定结构的数据还存在其他特殊方法。
关于你的申请,我没有任何具体的答案。很大程度上取决于数据的性质(采样率、噪声比等)。如果使用过多的平滑,则会弄脏数据或产生锯齿。如果您使用高阶多项式系数 K
过度拟合数据,也会出现同样的情况。在您的演示代码中,您还应该绘制 sin
函数的解析导数。然后使用不同量的输入噪声和平滑滤波器。如果您可以近似了解真实数据的各个方面,那么具有已知确切答案的工具可能会有所帮助。在实践中,我尝试使用尽可能少的平滑,以产生不太嘈杂的导数。通常,这意味着三阶多项式 (K = 3
) 和尽可能小的窗口大小 F
。
所以,是的,许多人建议您 use your eyes 来调整这些参数。然而,最近也有一些关于自动选择系数的研究: On the Selection of Optimum Savitzky-Golay Filters (2013)。 Savitzky-Golay 还有其他替代方案,例如基于 this paper 的 regularization ,但您可能需要自己在 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 。
关于matlab - sgolay函数在Matlab R2013a中如何工作?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/23943080/