Re: Convolving a non-uniformly sampled signal with a Gaussian

[Paul Kienzle]

Here's a function to compute the convolution of a linear spline R(Q) 
with a gaussian of width sigma calculated at Qo:

function [Ro,R] = gaussconv(Q,R,Qo,sigma)
  z = Qo - Q;
  G = exp(-z.*z/2/sigma^2);
  GR = G.*R;
  Gerf = erf(-z/sqrt(2)/sigma);
  m = diff(R)./diff(Q);
  b = R(2:end) - m .* Q(2:end);
  Ro = sum(0.5*(m*Qo+b).*diff(Gerf) - sigma/sqrt(2*pi)*m.*diff(G));

  ## Assuming the spline goes to 0 at Inf, we can approximate
  ## the convolution of the tail by multiplying the area of the
  ## tail of the gaussian by the height of the endpoint of the
  ## spline.  Similarly for -Inf.
  ## Alternatively, we can normalize by the area under the gaussian.
  ## In either case, we should be going out far enough in z that
  ## the contribution of the tails is negligible.
  if 0
    ## normalization
    Ro = 2 * Ro / (Gerf(end)-Gerf(1));
  end

%!test
%! Qo = 4; sigma = 1;
%! Q = [0,1,3,4,5,8];
%! R = sinc((Qo-Q)/2);
%! %plot(Q,R)
%! Qfine = linspace(0,8,10000);
%! Rfine = interp1(Q,R,Qfine);
%! Gfine = exp(-(Qfine-Qo).^2/2/sigma^2)/sqrt(2*pi*sigma^2);
%! %hold on; plot(Qfine,Rfine); hold off;
%! xfine = trapz(Qfine,Rfine.*Gfine);
%! x=gaussconv(Q,R,Qo,sigma);
%! assert(x,xfine,1e-8);


I have an oct-file for computing this for a set of points Q_i,sigma_i.  
Let me know if you are interested.

- Paul


On Jul 7, 2005, at 9:43 AM, Søren Hauberg wrote:

> Hi
> I have a signal that I need to convolve with a Gaussian (or another 
> similar filter). The problem is that the signal is not uniformly 
> sampled, meaning I can't use conv directly. My first thought was to 
> interpolate the signal, perform convolution, and then resample, but I 
> don't want to do this as the function that generated isn't continuos.
>
> Does anybody know a easy way to do this, or do I have to implement 
> this from scratch.
>
> /Søren
>
>
>
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
>
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------




-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------