2
$\begingroup$

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. Fig1

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?

$\endgroup$
  • $\begingroup$ The difference shows 0.00x difference, maximum, at the edges, so I suspect it's some numerical artifact due to the sin/cos based on the step size, or normalization, or such. I am just trying to make this in wxMaxima with bfloat, but it takes forever on my potato. If you think you have the patience... $\endgroup$ – a concerned citizen Mar 19 '18 at 6:45
  • $\begingroup$ I would expect any errors introduced by sin/cos or normalization would be on the order of the machine epsilon, ~1e-16, not 0.07.... so it is quite puzzling. $\endgroup$ – vibe Mar 19 '18 at 14:53
  • $\begingroup$ You also have the derivative, calculated as 1/sigma^2/L, which is ~0.076, then the convolution itself, which may introduce a bit, too, due to the relative large order and length of signal, then the step, 0.001, which is decimal, not binary, and may be represented with approximations, float or double, whichever it is. Also, I simply used the difference, not norm(), seems to be more precise in pointing the finger. But, sure, differences in this case can be strange. $\endgroup$ – a concerned citizen Mar 19 '18 at 15:39
  • $\begingroup$ Yes norm may not be the best. I just calculated the rms difference, which is on the order of 1e-4. It seems to be of that same order regardless of the step size $h$. $\endgroup$ – vibe Mar 19 '18 at 16:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.