gaussian smoothing function

FYI -  here is a gaussian smoothing function, via FFT.
This is an improvement on my previous quick-and-dirty version.

This one applies a Hanning filter (cosine window) to minimize effects
of finite window response, and also returns the output array on the
same grid as the input.  (Because of this, the interface has also
changed --- it is no longer necessary to return the grid.)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
isis> foo=gsmooth();

USAGE: smoothed_y = gsmooth(x_lo, x_hi, y, sigma)

Apply gaussian smoothing filter of width sigma to histogram array y.
  x_lo, x_hi are histogram bin edge coordinate arrays.
  y is histogram values array (x_lo,x_hi,y must have same lengths
METHOD: Convolution is performed via FFT, after applying Hanning filter.
RESTRICTIONS: grid must be uniform.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Two files follow: gsmooth.sl, and hanning.sl:
(local users can find these in ~dph/libisis)


%%%%%%%%%%%%%%%%%%%%%%%%%% gsmooth.sl   %%%%%%%%%%%%%
%
% 2000.10.29 dph
% v.1
% 2001.07.07
% v.2      ------ apply hanning filter; re-align on input x array.
%        NOTE: v.2 changed interface!

% 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?

% EXAMPLE use:
% > ()=evalfile("Sl/gsmooth.sl");
% > ()=evalfile("Sl/hanning.sl");
% >  sigma=0.004;    % approx HEG resolution
% >  % get your data into arrays...  xlo, xhi, counts
% >  scounts=gsmooth(xlo, xhi, counts, sigma);
% >  plot(xlo,xhi,counts,4); ohplot(xlo, xhi, scounts,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% example test case.
% x1 must be defined.  e.g., [0:8191:0.005]+1.  % MEG xlo.
%
%   evalfile(ilib+"/hanning.sl");
%   evalfile(ilib+"/gsmooth-v2.sl");
%   d=(sin(x1*PI*2/0.5)+5.) *exp(-x1/10);
%   d[[4000:5000:1]]=2.0;
%   d[4500]=10.;
%   sd=gsmooth(x1,x2,d,0.008);
%   xrange(0,45);hplot(x1,x2,d,4);ohplot(sx1,sx2,sd,3);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

require("hanning");

define gsmooth()
{
  variable 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;
  variable ymean;

  if (_NARGS != 4)
    {
    message("");
    message("USAGE: smoothed_y = gsmooth(x_lo, x_hi, y, sigma)");
    message("");
    message("Apply gaussian smoothing filter of width sigma to histogram array y.");
    message("  x_lo, x_hi are histogram bin edge coordinate arrays.");
    message("  y is histogram values array (x_lo,x_hi,y must have same lengths");
    message("METHOD: Convolution is performed via FFT, after applying Hanning filter.");
    message("RESTRICTIONS: grid must be uniform.");
    return 1;
    }

%+ pop arguments:
    sigma=();
    y=();
    xhi=();
    xlo=();
%-  

  % 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.*PI)*nsig); % eval kernel

% make array of kernel, same length as y:
%
  karray =  y*0.0;
  karray[ [0:nkbins-1:1] ] = k;

%%+ apply Hanning window (cosine)
   variable f;
   f = hanning(length(y));
   y = y*f;
%%-

%%+ forward fft
%
  (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 ;	                    % convolve: fft product.
%%-

%%+ apply phase shift of half-kernel width in order to return
%
% smoothed array on same x array as input:
% shift by nkbins/2
  variable phase_shift, phase_shift_theta;
  phase_shift_theta = 2.0*PI/length(y)  * [0:length(y)-1:1] * (nkbins/2);
  phase_shift = cos(phase_shift_theta) + sin(phase_shift_theta)*1i;
  c = c*phase_shift;
%%-  

%+ reverse fft
  (sy, isy) = fft1d(Real(c), Imag(c), 1);   % inverse fft;
%%-


%
  variable lz;   % indices of where filter is != 0.
  lz = where( f != 0.0);
  nrm = sum(y) / sum(sy[lz]); % ad hoc normalization; should be known a priori
  sy = sy * nrm;
  sy[lz] = sy[lz] / f[lz] ;   % remove filter.	

  return sy;

}

    

%%%%%%%%%%%% hanning.sl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% the Hanning window is:  1/2 * (1 - cos(2*PI*i/(N-1))), i=0..N-1
%
define hanning(length)
{
 variable result;

 result = 0.5 * (1.0 - cos( 2*PI/(length-1) * [0:length-1:1]) );

 return result;

} 

% modified hanning - ramp up to 1.0 quicker.
% n_end is number <= length/2.
% length is length of array
%
define mhanning(n_end, length)
{
  variable f1, result;

  f1=hanning(n_end*2);
  result = [0:length-1:1]*0.0 + 1.0;
  result[ [0:n_end-1:1] ] = f1[ [0:n_end-1:1] ];
  result[ [length-n_end:length-1:1] ] = f1[ [n_end:n_end*2-1:1] ];

  return result;
}

  






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