Re: Gaussian or other smoothing functions.

On May 28, 2:02 pm, Tim <tim...@[EMAIL PROTECTED]
> wrote:
> Hi All,
>
> I've looked in the do***entation and online, and I can't seem to find
> any way to do Gaussian smoothing.
>
> Normally, I'd simply write a function for it, but the problem is that
> this type of smoothing of course requires as many dimensions as there
> are data points. That is, each point of data is scaled as a function
> of the height of all other points. Since (I think) GNUPLOT reads one
> line at a time from data files, I don't know if it's possible.
>
> That said, GNUPLOT is so powerful, I can't see it not being able to do
> something like this. Could somebody please point me in the right
> direction?
>
> Cheers,
>
> Tim.


I guess we're all friends again :) Thanks to all of you... even
goog...@[EMAIL PROTECTED]
 :P

I wanted to give this follow-up post. I'm sure I'm not the only one
who has wanted to do some gaussian smoothing. I wrote this last week.
I haven't had a chance to comment it much, but I think it's straight
forward.  It's written in C, and should be ****table (don't forget -lm
for the math library).

It takes command line arguments of the standard deviation and kernel
radius (in that order). The standard deviation can be any floating
point value, but the kernel radius should be integral. The kernel
radius should be 3-5 times the standard deviation in order to preserve
the normality of the convolution, but keep in mind that the kernel
radius will be shaved off of the end of the data, since there is no
way to preserve normality at the boundary. One could easily program in
some boundary conditions (reflectance... or trickier... periodicity).

Anyhow, here is the code. I hope it helps someone:

// gsmooth.c
//
// Author: Timothy A.V. Teatro <timothy.teatro@[EMAIL PROTECTED]
>
// Date  : May 30, 2008
//
// Description:
//  This program reads data from standard in and convoludes it with a
//  discretized Gaussian lineshape. The output is sent to stdout and
//  error re****ting is minimal as the program is intended to be read
//  directly by GNUPLOT.
//
//  EDIT: This version is simplified for general use; to be posted on
//   comp.graphics.apps.gnuplot
//
//  Usage: ./gsmooth sigma krad < data.dat
//   sigma is the standard deviation of the Gaussian kernel:
//      g(x) = exp( -x^2/(2*sigma^2) ) / (2pi)^1/2 sigma
//   and krad is the radius of the kernel, in datapoints. That is,
//   over what radius the "averaging" is going to occur. In order to
//   preserve normality, krad should not be less than three times the
//   standard deviation.
//
//   Usage within gnuplot is shown in the following example:
//>    #set term post monochrome
//>    #set output "plotroughsmooth.ps"
//>
//>    #unset key
//>
//>    plot 'in.dat' using 1:2 title "Noise" w l, \
//>         '< gsmooth 2.00 7 < in.dat' using 1:2 title "Smoothed:
sig=2.00" w l
//
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double ReadDouble(FILE*file,char*flag);
double smoothY(double Y[], double kernel[], int krad);

int main(int argc, char *argv[]) {

   double sigma, *kernel, scale;
   double *X, *Y;
   int krad, i;
   char End;
   //
   // Parse Commandline arguments. We get the Standard Deviation
   //    and the kernel size.
   //
   sigma = atof(argv[1]);
   krad  = atoi(argv[2]);
   //
   // Now initialize the kernel
   //
   // Create the space in memory:
   //
   kernel=(double *)malloc( sizeof(double[2*krad+1]) );
   //
   // Define the scaling factor. Since exp(0)=1, we also assign the
value
   //    of scale to the center point (x=mu) in the kernel.
   //
   scale = 1/(sqrt(2*M_PI)*sigma);
   kernel[krad] = scale;
   //
   // Now define the rest of the kernal, symmetrically about the
center.
   //
   for (i=0; i<=krad-1; i++) {
      kernel[i]=scale*exp( -(i-krad)*(i-krad)/(2*sigma*sigma) );
      kernel[krad*2-i]=kernel[i];
//ec      printf("kernel[%d] = %f\n", i, kernel[i]);
//ec      printf("kernel[%d] = %f\n\n", krad*2-i, kernel[krad*2-i]);
   }
//ec   for (i=0; i<=krad*2; i++) {
//ec      printf("%d   %f\n", i-krad, kernel[i]);
//ec   }
   //
   // Now that the kernel is defined, we can start reading in data.
   //    We must read ahead of our calculation by krad entries.
   //
   //
   // Start by allocating space for the coordinates.
   //
   X = (double *) malloc( sizeof(double[2*krad+1]) );
   Y = (double *) malloc( sizeof(double[2*krad+1]) );
   //
   // Read off the first kernel diameter.
   //
   for (i=0; i<=2*krad; i++) {
      if (End == EOF) break;
      X[i] = ReadDouble(stdin, &End);
      Y[i] = ReadDouble(stdin, &End);
   }
   do {
      printf("%f\t%f\n", X[krad], smoothY(Y, kernel, krad));
      for (i=0; i<=2*krad-1; i++) {
         X[i]=X[i+1];
         Y[i]=Y[i+1];
      }
      X[2*krad] = ReadDouble(stdin, &End);
      Y[2*krad] = ReadDouble(stdin, &End);
   } while (End != EOF);
   return 0;

}

double ReadDouble(FILE*file,char*flag)
{
   char number[15];
   double Num;
   int i=0;
   char ch;
   do
   {
      ch = fgetc(file);
      if(ch!=' ' && ch!='\t' && ch!='\n' && ch!=EOF)
      {
         do
         {
            number[i]=ch;
            i++;
            ch = fgetc(file);
         }while(ch!=' ' && ch!='\t' && ch!='\n' && ch!=EOF);
         break;
      }
   } while (ch != EOF);
   *flag=ch;
   number[i]='\0';
   Num=atof(number);
   return Num;
}



double smoothY(double Y[], double kernel[], int krad) {

   int i;
   double G=0.;

   for (i=0; i<=krad * 2; i++) {
      G += Y[i] * kernel[i];
   }

   return G;
}