I am trying to smooth a time series signal with a Gaussian filter and then differentiate the signal (this is for an application for edge detection). A nice property of convolution is: $$ \frac{d}{dx} \left( s(t) * g(t) \right) = \frac{ds(t)}{dt} * g(t) = s(t) * \frac{dg(t)}{dt} $$ When $g(t)$ is a Gaussian kernel, $dg(t)/dt$ is analytically computable and so I don't need to use any finite differences to compute $ds(t)/dt$. But in doing some numerical tests on the above equation, I found some curious differences between the right-most equality in the equation above. The following matlab code illustrates the problem.
It defines a simple sine wave signal $s(t) = sin(t)$, then computes a 211 point Gaussian kernel along with its derivative. It then convolves $s(t)$ with $dg/dt$ and then convolves $ds/dt$ with $g(t)$ and compares the solutions. I would expect the solutions to be identical based on the known result above, but surprisingly I find small differences between the two solutions. Can anyone shed any light on what is going on?
h = 0.001; % step size for signal
K = 211; % window size of gaussian
sigma = 0.25; % stddev of gaussian window
% define signal and its derivative
x = -pi:h:pi;
signal = sin(x)';
dsignal = cos(x)';
% make Gaussian window, size K
L = K - 1;
n = [-1:2/L:1]';
g = exp(-(1/2)*(n/sigma).^2);
% now its derivative
dg = -(1/sigma^2/(L/2)) * n .* g / h;
% normalize so convolution has unit gain
w = g / sum(g);
dw = dg / sum(g);
% compute sol1 = d/dx signal(x) * G(x)
sol1 = conv(dsignal, w, 'same');
% compute sol2 = signal(x) * dG(x)/dx
sol2 = conv(signal, dw, 'same');
% plot solutions - they are visually identical
plot(x,sol1,x,sol2+1e-2)
% ...but this tells another story
r = sol1 - sol2;
rmsdiff = sqrt(dot(r,r) / length(r))
maxdiff = max(r)
mindiff = min(r)
The output of the matlab code is:
rmsdiff =
9.266306205125738e-04
maxdiff =
9.940615138074316e-04
mindiff =
-0.008096374978515
where I'd expect all of these to be zero. A plot of the two solutions, offset by a small epsilon is below.
So the two solutions look the same, but in fact have small differences. Since everything in the script uses analytic results (no finite differencing), I would expect the two solutions to be identical. Can anyone explain why they are not?
bfloat
, but it takes forever on my potato. If you think you have the patience... $\endgroup$norm()
, seems to be more precise in pointing the finger. But, sure, differences in this case can be strange. $\endgroup$