fftw smoothing

I'm sending this to isis-users so somebody can improve it ;-)

It's a Fourier, Gaussian-kernel smoother, useful for getting rid of
that ragged structure in your under-binned data. It uses the isis
hooks to fftw in a quick-and-dirty way.  I haven't done any apodizing
to take care of edge affects, nor have I carefully considered
coordinate offset or kernel-centering on the data to put the output on
the same grid as input (hence the x-array fudge and output
x-coordinates).  But it works well enough to make pretty pictures
(don't do quantitative measurements on smoothed data!).

-- Dave


% 2000.10.29 dph
% v.1

% gaussian smoothing
%
% convolve y by gaussian kernel, width sigma, in x's units
% Use fft1d
% Assumes uniform grid.
% Does not worry about integral form of y.  Should it?

% not checked for quantitative accuracy; checked for visualization.
% not checked for edge affects; probably should subtract mean and
%     apply cosine bell to ends 
%    (works OK for counts, which fade away at array ends; probably
%    won't work for fluxed data)

% example use:
%
% > ()=evalfile("Sl/gsmooth.sl");
% >  sigma=0.004;    % approx HEG resolution
% >  % get your data into arrays...  xlo, xhi, counts
% >  (sxlo,sxhi,scounts)=gsmooth(xlo, xhi, counts, sigma);
% >  hplot(sxlo, sxhi, scounts);


define gsmooth(xlo,xhi,y,sigma)
{
  variable k, nk, karray, dx, nsig;
  variable nkbins;
  variable rffty, iffty, rfftk, ifftk, cffty, cfftk, c, sy, isy, nrm;
  variable xlo_out, xhi_out;

  % make kernel
  dx =  xhi[0]-xlo[0];
  nsig = sigma/dx;          % # bins in sigma.
  nkbins = int(4.*nsig)+1;          % # bins to use in the gaussian kernel; +-2sigma.

  k=exp( -([0:nkbins-1:1]-nkbins/2.)^2/ nsig^2/2.) / (sqrt(2.*3.1415926)*nsig); % eval kernel

  karray =  y*0.0;
  karray[ [0:nkbins-1:1] ] = k;

  (rffty, iffty) = fft1d(y, y*0., -1);               % fft of data
  (rfftk, ifftk) = fft1d(karray, karray*0., -1);     % fft of kernel

  cffty = rffty + iffty * 1i;     % convert to Slang complex types.
  cfftk = rfftk + ifftk * 1i;

  c = cffty * cfftk ;	                    % fft product.
  (sy, isy) = fft1d(Real(c), Imag(c), 1);   % inverse fft;

  nrm = sum(y)/sum(sy);       % ad hoc normalization; should be known a priori
  sy = sy * nrm;

  xlo_out = xlo - nkbins/2. * dx;            % fudge the x-array 
  xhi_out = xlo_out + dx;

  return xlo_out, xhi_out, sy;

}

    

  

----
You received this message because you are
subscribed to the isis-users list.
To unsubscribe, send a message to
isis-users-request_at_email.domain.hiddenwith the first line of the message as:
unsubscribe