scatter plots smooth算法 lowess_smooth rloess c++-程序员宅基地

技术标签: 算法  

/*
 *  c++ implementation of Lowess weighted regression by 
 *  Peter Glaus http://www.cs.man.ac.uk/~glausp/
 *
 *
 *  Based on fortran code by Cleveland downloaded from:
 *  http://netlib.org/go/lowess.f
 *  original author:
* [email protected] Mon Dec 30 16:55 EST 1985
* W. S. Cleveland
* Bell Laboratories
* Murray Hill NJ 07974
 *  
 *  See original documentation below the code for details.
 * 
 */
#include<algorithm>
#include<cmath>
#include<fstream>

using namespace std;

#include "lowess.h"
#include "common.h"

void lowess(const vector<double> &x, const vector<double> &y, double f, long nsteps, vector<double> &ys){//{
   {
   {
   vector<double> rw,res;
   lowess(x,y,f,nsteps,0.,ys,rw,res);
}//}}}
void lowess(const vector<double> &x, const vector<double> &y, double f, long nsteps, double delta, vector<double> &ys, vector<double> &rw, vector<double>&res){ //{
   {
   {
   long n=(long)x.size();
   bool ok=false;
   long nleft,nright, i, j, iter, last, m1, m2, ns;
   double cut, cmad, r, d1, d2, c1, c9, alpha, denom;
   if((n==0)||((long)y.size()!=n)) return;
   ys.resize(n);
   rw.resize(n);
   res.resize(n);
   if(n==1){
      ys[0]=y[0];
      return;
   }
   // ns - at least 2, at most n
   ns = max(min((long)(f*n),n),(long)2);
   for(iter=0;iter<nsteps+1; iter++){
      // robustnes iterations
      nleft = 0;
      nright = ns-1;
      // index of last estimated point
      last = -1;
      // index of current point
      i=0;
      do{
         while(nright<n-1){
            // move <nleft,nright> right, while radius decreases
            d1 = x[i]-x[nleft];
            d2 = x[nright+1] - x[i];
            if(d1<=d2)break;
            nleft++;
            nright++;
         }
         // fit value at x[i]
         lowest(x,y,x[i],ys[i],nleft,nright,res,iter>0,rw,ok);
         if(!ok) ys[i]=y[i];
         if(last<i-1){
            // interpolate skipped points
            if(last<0){
               warning("Lowess: out of range.\n");
            }
            denom = x[i] - x[last];
            for(j=last+1;j<i;j++){
               alpha = (x[j]-x[last])/denom;
               ys[j] = alpha * ys[i] + (1.0-alpha)*ys[last];
            }
         }
         last = i;
         cut = x[last]+delta;
         for(i=last+1;i<n;i++){
            if(x[i]>cut)break;
            if(x[i]==x[last]){
               ys[i]=ys[last];
               last=i;
            }
         }
         i=max(last+1,i-1);
      }while(last<n-1);
      for(i=0;i<n;i++)
         res[i] = y[i]-ys[i];
      if(iter==nsteps)break ;
      for(i=0;i<n;i++)
         rw[i]=abs(res[i]);
      sort(rw.begin(),rw.end());
      m1 = n/2+1;
      m2 = n-m1;
      m1 --;
      cmad = 3.0 *(rw[m1]+rw[m2]);
      c9 = .999*cmad;
      c1 = .001*cmad;
      for(i=0;i<n;i++){
         r = abs(res[i]);
         if(r<=c1) rw[i]=1;
         else if(r>c9) rw[i]=0;
         else rw[i] = (1.0-(r/cmad)*(r/cmad))*(1.0-(r/cmad)*(r/cmad));
      }
   }
}//}}}

void lowest(const vector<double> &x, const vector<double> &y, double xs, double &ys, long nleft, long nright, vector<double> &w, bool userw,  vector<double> &rw, bool &ok){//{
   {
   {
   long n = (long)x.size();
   long nrt, j;
   double a, b, c, h, r, h1, h9, range;
   range = x[n-1]-x[0];
   h = max(xs-x[nleft],x[nright]-xs);
   h9 = 0.999*h;
   h1 = 0.001*h;
   // sum of weights
   a = 0; 
   for(j=nleft;j<n;j++){
      // compute weights (pick up all ties on right)
      w[j]=0.;
      r = abs(x[j]-xs);
      if(r<=h9){
         // small enough for non-zero weight
         if(r>h1) w[j] = (1.0-(r/h)*(r/h)*(r/h))*(1.0-(r/h)*(r/h)*(r/h))*(1.0-(r/h)*(r/h)*(r/h));
         else w[j] = 1.;
         if(userw) w[j] *= rw[j];
         a += w[j];
      }else if(x[j]>xs) break; // get out at first zero wt on right
   }
   nrt = j-1;
   // rightmost pt (may be greater than nright because of ties)
   if(a<=0.) ok = false;
   else{
      // weighted least squares
      ok = true;
      // normalize weights
      for(j=nleft;j<=nrt;j++)
         w[j] /= a;
      if(h>0.){
         // use linear fit
         a = 0.;
         for(j=nleft;j<=nrt;j++)
            a += w[j]*x[j]; // weighted centre of values
         b = xs-a;
         c = 0;
         for(j=nleft;j<=nrt;j++)
            c += w[j]*(x[j]-a)*(x[j]-a);
         if(sqrt(c)>0.001*range){
            // points are spread enough to compute slope
            b /= c;
            for(j=nleft;j<=nrt;j++)
               w[j] *= (1.0+b*(x[j]-a));
         }
      }
      ys = 0;
      for(j=nleft;j<=nrt;j++)
         ys += w[j]*y[j];
   }
}//}}}

/* {
   {
   { Documentation
* [email protected] Mon Dec 30 16:55 EST 1985
* W. S. Cleveland
* Bell Laboratories
* Murray Hill NJ 07974
* 
* outline of this file:
*    lines 1-72   introduction
*        73-177   documentation for lowess
*       178-238   ratfor version of lowess
*       239-301   documentation for lowest
*       302-350   ratfor version of lowest
*       351-end   test driver and fortran version of lowess and lowest
* 
*   a multivariate version is available by "send dloess from a"
* 
*              COMPUTER PROGRAMS FOR LOCALLY WEIGHTED REGRESSION
* 
*             This package consists  of  two  FORTRAN  programs  for
*        smoothing    scatterplots   by   robust   locally   weighted
*        regression, or lowess.   The  principal  routine  is  LOWESS
*        which   computes   the  smoothed  values  using  the  method
*        described in The Elements of Graphing Data, by William S.
*        Cleveland    (Wadsworth,    555 Morego   Street,   Monterey,
*        California 93940).
* 
*             LOWESS calls a support routine, LOWEST, the code for
*        which is included. LOWESS also calls a routine  SORT,  which
*        the user must provide.
* 
*             To reduce the computations, LOWESS  requires  that  the
*        arrays  X  and  Y,  which  are  the  horizontal and vertical
*        coordinates, respectively, of the scatterplot, be such  that
*        X  is  sorted  from  smallest  to  largest.   The  user must
*        therefore use another sort routine which will sort X  and  Y
*        according  to X.
*             To summarize the scatterplot, YS,  the  fitted  values,
*        should  be  plotted  against X.   No  graphics  routines are
*        available in the package and must be supplied by the user.
* 
*             The FORTRAN code for the routines LOWESS and LOWEST has
*        been   generated   from   higher   level   RATFOR   programs
*        (B. W. Kernighan, ``RATFOR:  A Preprocessor for  a  Rational
*        Fortran,''  Software Practice and Experience, Vol. 5 (1975),
*        which are also included.
* 
*             The following are data and output from LOWESS that  can
*        be  used  to check your implementation of the routines.  The
*        notation (10)v means 10 values of v.
* 
* 
* 
* 
*        X values:
*          1  2  3  4  5  (10)6  8  10  12  14  50
* 
*        Y values:
*           18  2  15  6  10  4  16  11  7  3  14  17  20  12  9  13  1  8  5  19
* 
* 
*        YS values with F = .25, NSTEPS = 0, DELTA = 0.0
*         13.659  11.145  8.701  9.722  10.000  (10)11.300  13.000  6.440  5.596
*           5.456  18.998
* 
*        YS values with F = .25, NSTEPS = 0 ,  DELTA = 3.0
*          13.659  12.347  11.034  9.722  10.511  (10)11.300  13.000  6.440  5.596
*            5.456  18.998
* 
*        YS values with F = .25, NSTEPS = 2, DELTA = 0.0
*          14.811  12.115  8.984  9.676  10.000  (10)11.346  13.000  6.734  5.744
*            5.415  18.998
* 
* 
* 
* 
*                                   LOWESS
* 
* 
* 
*        Calling sequence
* 
*        CALL LOWESS(X,Y,N,F,NSTEPS,DELTA,YS,RW,RES)
* 
*        Purpose
* 
*        LOWESS computes the smooth of a scatterplot of Y  against  X
*        using  robust  locally  weighted regression.  Fitted values,
*        YS, are computed at each of the  values  of  the  horizontal
*        axis in X.
* 
*        Argument description
* 
*              X = Input; abscissas of the points on the
*                  scatterplot; the values in X must be ordered
*                  from smallest to largest.
*              Y = Input; ordinates of the points on the
*                  scatterplot.
*              N = Input; dimension of X,Y,YS,RW, and RES.
*              F = Input; specifies the amount of smoothing; F is
*                  the fraction of points used to compute each
*                  fitted value; as F increases the smoothed values
*                  become smoother; choosing F in the range .2 to
*                  .8 usually results in a good fit; if you have no
*                  idea which value to use, try F = .5.
*         NSTEPS = Input; the number of iterations in the robust
*                  fit; if NSTEPS = 0, the nonrobust fit is
*                  returned; setting NSTEPS equal to 2 should serve
*                  most purposes.
*          DELTA = input; nonnegative parameter which may be used
*                  to save computations; if N is less than 100, set
*                  DELTA equal to 0.0; if N is greater than 100 you
*                  should find out how DELTA works by reading the
*                  additional instructions section.
*             YS = Output; fitted values; YS(I) is the fitted value
*                  at X(I); to summarize the scatterplot, YS(I)
*                  should be plotted against X(I).
*             RW = Output; robustness weights; RW(I) is the weight
*                  given to the point (X(I),Y(I)); if NSTEPS = 0,
*                  RW is not used.
*            RES = Output; residuals; RES(I) = Y(I)-YS(I).
* 
* 
*        Other programs called
* 
*               LOWEST
*               SSORT
* 
*        Additional instructions
* 
*        DELTA can be used to save computations.   Very  roughly  the
*        algorithm  is  this:   on the initial fit and on each of the
*        NSTEPS iterations locally weighted regression fitted  values
*        are computed at points in X which are spaced, roughly, DELTA
*        apart; then the fitted values at the  remaining  points  are
*        computed  using  linear  interpolation.   The  first locally
*        weighted regression (l.w.r.) computation is carried  out  at
*        X(1)  and  the  last  is  carried  out at X(N).  Suppose the
*        l.w.r. computation is carried out at  X(I).   If  X(I+1)  is
*        greater  than  or  equal  to  X(I)+DELTA,  the  next  l.w.r.
*        computation is carried out at X(I+1).   If  X(I+1)  is  less
*        than X(I)+DELTA, the next l.w.r.  computation is carried out
*        at the largest X(J) which is greater than or equal  to  X(I)
*        but  is not greater than X(I)+DELTA.  Then the fitted values
*        for X(K) between X(I)  and  X(J),  if  there  are  any,  are
*        computed  by  linear  interpolation  of the fitted values at
*        X(I) and X(J).  If N is less than 100 then DELTA can be  set
*        to  0.0  since  the  computation time will not be too great.
*        For larger N it is typically not necessary to carry out  the
*        l.w.r.  computation for all points, so that much computation
*        time can be saved by taking DELTA to be  greater  than  0.0.
*        If  DELTA =  Range  (X)/k  then,  if  the  values  in X were
*        uniformly  scattered  over  the  range,  the   full   l.w.r.
*        computation  would be carried out at approximately k points.
*        Taking k to be 50 often works well.
* 
*        Method
* 
*        The fitted values are computed by using the nearest neighbor
*        routine  and  robust locally weighted regression of degree 1
*        with the tricube weight function.  A few additional features
*        have  been  added.  Suppose r is FN truncated to an integer.
*        Let  h  be  the  distance  to  the  r-th  nearest   neighbor
*        from X(I).   All  points within h of X(I) are used.  Thus if
*        the r-th nearest neighbor is exactly the  same  distance  as
*        other  points,  more  than r points can possibly be used for
*        the smooth at  X(I).   There  are  two  cases  where  robust
*        locally  weighted regression of degree 0 is actually used at
*        X(I).  One case occurs when  h  is  0.0.   The  second  case
*        occurs  when  the  weighted  standard error of the X(I) with
*        respect to the weights w(j) is  less  than  .001  times  the
*        range  of the X(I), where w(j) is the weight assigned to the
*        j-th point of X (the tricube  weight  times  the  robustness
*        weight)  divided by the sum of all of the weights.  Finally,
*        if the w(j) are all zero for the smooth at X(I), the  fitted
*        value is taken to be Y(I).
* 
* 
* 
* 
*  subroutine lowess(x,y,n,f,nsteps,delta,ys,rw,res)
*  real x(n),y(n),ys(n),rw(n),res(n)
*  logical ok
*  if (n<2){ ys(1) = y(1); return }
*  ns = max0(min0(ifix(f*float(n)),n),2)  # at least two, at most n points
*  for(iter=1; iter<=nsteps+1; iter=iter+1){      # robustness iterations
*         nleft = 1; nright = ns
*         last = 0        # index of prev estimated point
*         i = 1   # index of current point
*         repeat{
*                 while(nright<n){
*  # move nleft, nright to right if radius decreases
*                         d1 = x(i)-x(nleft)
*                         d2 = x(nright+1)-x(i)
*  # if d1<=d2 with x(nright+1)==x(nright), lowest fixes
*                         if (d1<=d2) break
*  # radius will not decrease by move right
*                         nleft = nleft+1
*                         nright = nright+1
*                         }
*                 call lowest(x,y,n,x(i),ys(i),nleft,nright,res,iter>1,rw,ok)
*  # fitted value at x(i)
*                 if (!ok) ys(i) = y(i)
*  # all weights zero - copy over value (all rw==0)
*                 if (last<i-1) { # skipped points -- interpolate
*                         denom = x(i)-x(last)    # non-zero - proof?
*                         for(j=last+1; j<i; j=j+1){
*                                 alpha = (x(j)-x(last))/denom
*                                 ys(j) = alpha*ys(i)+(1.0-alpha)*ys(last)
*                                 }
*                         }
*                 last = i        # last point actually estimated
*                 cut = x(last)+delta     # x coord of close points
*                 for(i=last+1; i<=n; i=i+1){     # find close points
*                         if (x(i)>cut) break     # i one beyond last pt within cut
*                         if(x(i)==x(last)){      # exact match in x
*                                 ys(i) = ys(last)
*                                 last = i
*                                 }
*                         }
*                 i=max0(last+1,i-1)
*  # back 1 point so interpolation within delta, but always go forward
*                 } until(last>=n)
*         do i = 1,n      # residuals
*                 res(i) = y(i)-ys(i)
*         if (iter>nsteps) break  # compute robustness weights except last time
*         do i = 1,n
*                 rw(i) = abs(res(i))
*         call sort(rw,n)
*         m1 = 1+n/2; m2 = n-m1+1
*         cmad = 3.0*(rw(m1)+rw(m2))      # 6 median abs resid
*         c9 = .999*cmad; c1 = .001*cmad
*         do i = 1,n {
*                 r = abs(res(i))
*                 if(r<=c1) rw(i)=1.      # near 0, avoid underflow
*                 else if(r>c9) rw(i)=0.  # near 1, avoid underflow
*                 else rw(i) = (1.0-(r/cmad)**2)**2
*                 }
*         }
*  return
*  end
* 
* 
* 
* 
*                                   LOWEST
* 
* 
* 
*        Calling sequence
* 
*        CALL LOWEST(X,Y,N,XS,YS,NLEFT,NRIGHT,W,USERW,RW,OK)
* 
*        Purpose
* 
*        LOWEST is a support routine for LOWESS and  ordinarily  will
*        not  be  called  by  the  user.   The  fitted  value, YS, is
*        computed  at  the  value,  XS,  of  the   horizontal   axis.
*        Robustness  weights,  RW,  can  be employed in computing the
*        fit.
* 
*        Argument description
* 
* 
*              X = Input; abscissas of the points on the
*                  scatterplot; the values in X must be ordered
*                  from smallest to largest.
*              Y = Input; ordinates of the points on the
*                  scatterplot.
*              N = Input; dimension of X,Y,W, and RW.
*             XS = Input; value of the horizontal axis at which the
*                  smooth is computed.
*             YS = Output; fitted value at XS.
*          NLEFT = Input; index of the first point which should be
*                  considered in computing the fitted value.
*         NRIGHT = Input; index of the last point which should be
*                  considered in computing the fitted value.
*              W = Output; W(I) is the weight for Y(I) used in the
*                  expression for YS, which is the sum from
*                  I = NLEFT to NRIGHT of W(I)*Y(I); W(I) is
*                  defined only at locations NLEFT to NRIGHT.
*          USERW = Input; logical variable; if USERW is .TRUE., a
*                  robust fit is carried out using the weights in
*                  RW; if USERW is .FALSE., the values in RW are
*                  not used.
*             RW = Input; robustness weights.
*             OK = Output; logical variable; if the weights for the
*                  smooth are all 0.0, the fitted value, YS, is not
*                  computed and OK is set equal to .FALSE.; if the
*                  fitted value is computed OK is set equal to
* 
* 
*        Method
* 
*        The smooth at XS is computed using (robust) locally weighted
*        regression of degree 1.  The tricube weight function is used
*        with h equal to the maximum of XS-X(NLEFT) and X(NRIGHT)-XS.
*        Two  cases  where  the  program  reverts to locally weighted
*        regression of degree 0 are described  in  the  documentation
*        for LOWESS.
* 
* 
* 
* 
*  subroutine lowest(x,y,n,xs,ys,nleft,nright,w,userw,rw,ok)
*  real x(n),y(n),w(n),rw(n)
*  logical userw,ok
*  range = x(n)-x(1)
*  h = amax1(xs-x(nleft),x(nright)-xs)
*  h9 = .999*h
*  h1 = .001*h
*  a = 0.0        # sum of weights
*  for(j=nleft; j<=n; j=j+1){     # compute weights (pick up all ties on right)
*         w(j)=0.
*         r = abs(x(j)-xs)
*         if (r<=h9) {    # small enough for non-zero weight
*                 if (r>h1) w(j) = (1.0-(r/h)**3)**3
*                 else      w(j) = 1.
*                 if (userw) w(j) = rw(j)*w(j)
*                 a = a+w(j)
*                 }
*         else if(x(j)>xs)break   # get out at first zero wt on right
*         }
*  nrt=j-1        # rightmost pt (may be greater than nright because of ties)
*  if (a<=0.0) ok = FALSE
*  else { # weighted least squares
*         ok = TRUE
*         do j = nleft,nrt
*                 w(j) = w(j)/a   # make sum of w(j) == 1
*         if (h>0.) {     # use linear fit
*                 a = 0.0
*                 do j = nleft,nrt
*                         a = a+w(j)*x(j) # weighted center of x values
*                 b = xs-a
*                 c = 0.0
*                 do j = nleft,nrt
*                         c = c+w(j)*(x(j)-a)**2
*                 if(sqrt(c)>.001*range) {
*  # points are spread out enough to compute slope
*                         b = b/c
*                         do j = nleft,nrt
*                                 w(j) = w(j)*(1.0+b*(x(j)-a))
*                         }
*                 }
*         ys = 0.0
*         do j = nleft,nrt
*                 ys = ys+w(j)*y(j)
*         }
*  return
*  end
* 
}}}*/




C语音版本

#define FALSE 0
#define TRUE 1

long min (long x, long y)
{
	return((x < y) ? x : y);
}

long max (long x, long y)
{
	return((x > y) ? x : y);
}

static double pow2(double x) { return(x * x); }
static double pow3(double x) { return(x * x * x); }
double fmax(double x, double y) { return (x > y ? x : y); }

int static compar(const void *aa, const void *bb)
{
	const double* a = (double*)aa;
	const double* b = (double*)bb;

	if (*a < *b) return(-1);
	else if (*a > *b) return(1);
	else return(0);
}

static void lowest(const double *x, const double *y, size_t n, double xs, double *ys, long nleft, long nright,
	double *w, int userw, double *rw, int *ok)
{
	double range, h, h1, h9, a, b, c, r;
	long j, nrt;

	range = x[n - 1] - x[0];
	h = fmax(xs - x[nleft], x[nright] - xs);
	h9 = .999 * h;
	h1 = .001 * h;

	/* compute weights (pick up all ties on right) */
	a = 0.0;        /* sum of weights */
	for(j = nleft; j < n; j++)
	{
		w[j]=0.0;
		r = fabs(x[j] - xs);
		if (r <= h9) 
		{    /* small enough for non-zero weight */
			if (r > h1) w[j] = pow3(1.0-pow3(r/h));
			else w[j] = 1.0;
			if (userw) w[j] = rw[j] * w[j];
			a += w[j];
		}
		else if (x[j] > xs) break;  /* get out at first zero wt on right */
	}
	nrt = j - 1;  /* rightmost pt (may be greater than nright because of ties) */
	if (a <= 0.0) *ok = FALSE;
	else { /* weighted least squares */
		*ok = TRUE;

		/* make sum of w[j] == 1 */
		for (j = nleft; j <= nrt; j++) w[j] = w[j] / a;

		if (h > 0.0) {     /* use linear fit */

			/* find weighted center of x values */
			for (j = nleft, a = 0.0; j <= nrt; j++) a += w[j] * x[j];

			b = xs - a;
			for (j = nleft, c = 0.0; j <= nrt; j++) 
				c += w[j] * (x[j] - a) * (x[j] - a);

			if(sqrt(c) > .001 * range) {
				/* points are spread out enough to compute slope */
				b = b/c;
				for (j = nleft; j <= nrt; j++) 
					w[j] = w[j] * (1.0 + b*(x[j] - a));
			}
		}
		for (j = nleft, *ys = 0.0; j <= nrt; j++) *ys += w[j] * y[j];
	}
}

static void sort(double *x, size_t n)
{
	qsort(x, n, sizeof(double), compar);
}

int lowess(const double *x, const double *y, size_t n,double f, size_t nsteps,
	double delta, double *ys, double *rw, double *res)
{
	int iter, ok;
	long i, j, last, m1, m2, nleft, nright, ns;
	double d1, d2, denom, alpha, cut, cmad, c9, c1, r;

	if (n < 2) { ys[0] = y[0]; return(1); }
	ns = max(min((long) (f * n), n), 2);  /* at least two, at most n points */
	for(iter = 1; iter <= nsteps + 1; iter++){      /* robustness iterations */
		nleft = 0; nright = ns - 1;
		last = -1;        /* index of prev estimated point */
		i = 0;   /* index of current point */
		do {
			while(nright < n - 1){
				/* move nleft, nright to right if radius decreases */
				d1 = x[i] - x[nleft];
				d2 = x[nright + 1] - x[i];
				/* if d1 <= d2 with x[nright+1] == x[nright], lowest fixes */
				if (d1 <= d2) break;
				/* radius will not decrease by move right */
				nleft++;
				nright++;
			}
			lowest(x, y,
				n, x[i],
				&ys[i],
				nleft, nright,
				res, (iter > 1), rw, &ok);
			/* fitted value at x[i] */
			if (! ok) ys[i] = y[i];
			/* all weights zero - copy over value (all rw==0) */
			if (last < i - 1) { /* skipped points -- interpolate */
				denom = x[i] - x[last];    /* non-zero - proof? */
				for(j = last + 1; j < i; j = j + 1){
					alpha = (x[j] - x[last]) / denom;
					ys[j] = alpha * ys[i] + (1.0 - alpha) * ys[last];
				}
			}
			last = i;        /* last point actually estimated */
			cut = x[last] + delta;     /* x coord of close points */
			for(i=last + 1; i < n; i++) {     /* find close points */
				if (x[i] > cut) break;     /* i one beyond last pt within cut */
				if(x[i] == x[last]) {      /* exact match in x */
					ys[i] = ys[last];
					last = i;
				}
			}
			i = max(last + 1,i - 1);
			/* back 1 point so interpolation within delta, but always go forward */
		} while(last < n - 1);
		for (i = 0; i < n; i++)      /* residuals */
			res[i] = y[i] - ys[i];
		if (iter > nsteps) break; /* compute robustness weights except last time */
		for (i = 0; i < n; i++) 
			rw[i] = fabs(res[i]);
		sort(rw,n);
		m1 = 1 + n / 2; m2 = n - m1 + 1;
		cmad = 3.0 * (rw[m1] + rw[m2]);      /* 6 median abs resid */
		c9 = .999 * cmad; c1 = .001 * cmad;
		for (i = 0; i < n; i++) {
			r = fabs(res[i]);
			if(r <= c1) rw[i] = 1.0;      /* near 0, avoid underflow */
			else if(r > c9) rw[i] = 0.0;  /* near 1, avoid underflow */
			else rw[i] = pow2(1.0 - pow2(r / cmad));
		}
	}
	return(0);
}

测试main

void test_smooth()
{
	const double in[] = {
		-55.1047,
		-56.3673,
		-56.314,
		-55.8626,
		-56.3733,
		-55.8694,
		-54.4476,
		-56.1106,
		-58.3752,
		-56.4632,
		-57.5418,
		-57.259,
		-54.999,
		-56.5903,
		-56.5675,
		-57.2702,
		-56.5198,
		-59.275,
		-58.0706,
		-55.4509,
		-58.1618,
		-57.2596,
		-55.4555,
		-55.7893,
		-55.7953,
		-55.8368,
		-57.3642,
		-57.5559,
		-56.0514,
		-56.8001,
		-58.4702,
		-55.175,
		-53.7347,
		-54.4519,
		-54.5773,
		-56.9627,
		-57.1959,
		-55.6704,
		-55.0481,
		-57.9426,
		-57.3462,
		-55.6331,
		-55.629,
		-55.3975,
		-56.4719,
		-58.1078,
		-56.1705,
		-54.758,
		-58.8711,
		-57.9153,
		-56.4004,
		-55.1685,
		-55.7985,
		-54.3574,
		-55.6311,
		-55.4626,
		-56.6099,
		-55.4795,
		-54.0479,
		-56.069,
		-56.2238,
		-56.643,
		-57.3297,
		-57.2569,
		-55.3871,
		-55.8629,
		-55.827,
		-55.3129,
		-56.7753,
		-54.7421,
		-53.2383,
		-53.1972,
		-54.2125,
		-55.1294,
		-55.3516,
		-54.4107,
		-56.1692,
		-55.6607,
		-54.1289,
		-56.2226,
		-54.9853,
		-54.5406,
		-55.8668,
		-54.4344,
		-51.34,
		-53.4457,
		-55.3933,
		-56.4022,
		-57.494,
		-57.042,
		-53.8239,
		-54.7248,
		-55.1078,
		-54.9422,
		-57.6964,
		-57.2593,
		-54.7605,
		-56.342,
		-57.4363,
		-53.8504,
		-52.5132,
		-54.1004,
		-55.4099,
		-55.062,
		-54.2594,
		-52.8564,
		-52.14,
		-51.7081,
		-52.2992,
		-52.3724,
		-51.8007,
		-55.592,
		-55.7873,
		-53.5016,
		-54.1608,
		-53.7607,
		-52.8233,
		-54.0887,
		-54.6511,
		-54.4701,
		-54.7901,
		-56.7217,
		-55.1668,
		-54.6029,
		-55.2335,
		-53.67,
		-53.9694,
		-56.1003,
		-55.6258,
		-55.8682,
		-54.7942,
		-54.6461,
		-56.4503,
		-57.3815,
		-55.4552,
		-56.4655,
		-55.3854,
		-54.1829,
		-53.3174,
		-54.3715,
		-53.5423,
		-54.4937,
		-56.8845,
		-54.4562,
		-53.0783,
		-54.8609,
		-52.7257,
		-53.1482,
		-55.1311,
		-54.7786,
		-54.3794,
		-55.2594,
		-52.3897,
		-50.3529,
		-50.6989,
		-50.8013,
		-50.2843,
		-51.4467,
		-52.3954,
		-51.4057,
		-51.6931,
		-53.6986,
		-52.1103,
		-49.9167,
		-51.4783,
		-53.218,
		-53.8505,
		-52.805,
		-51.589,
		-51.8991,
		-53.1859,
		-50.7663,
		-51.6103,
		-52.6432,
		-50.0238,
		-52.5902,
		-54.4426,
		-51.1014,
		-51.337,
		-52.8024,
		-53.7283,
		-53.3006,
		-54.6558,
		-56.1147,
		-53.3179,
		-52.8044,
		-52.334,
		-51.9468,
		-51.2366,
		-56.9915,
		-54.5127,
		-53.0841,
		-54.4758,
		-53.2447,
		-54.7147,
		-54.3793,
		-52.4646,
		-52.7251,
		-53.2872,
		-51.1914,
		-52.1654,
		-53.3249,
		-52.2418,
		-50.3992,
		-51.7577,
		-50.617,
		-50.6632,
		-54.7326,
		-52.635,
		-51.3294,
		-54.1903,
		-53.3178,
		-51.235,
		-53.1232,
		-52.5514,
		-51.5221,
		-49.9557,
		-52.2744,
		-53.2325,
		-51.3947,
		-51.9394,
		-52.0155,
		-51.9813,
		-52.9384,
		-51.6752,
		-51.4657,
		-53.9598,
		-56.0678,
		-55.6356,
		-54.9773,
		-52.1095,
		-49.7851,
		-51.2385,
		-52.6269,
		-53.3314,
		-55.0205,
		-55.7239,
		-53.4701,
		-52.5249,
		-53.1064,
		-55.33,
		-53.1046,
		-52.853,
		-53.7369,
		-54.7797,
		-54.5858,
		-54.2175,
		-53.2216,
		-50.8936,
		-36.8913,
		-35.1439,
		-37.0516,
		-50.6355,
		-52.1987,
		-53.0451,
		-53.1897,
		-52.8646,
		-52.6694,
		-52.9763,
		-53.23,
		-54.3301,
		-55.2726,
		-54.0729,
		-51.3799,
		-50.7924,
		-51.3911,
		-51.1238,
		-50.1222,
		-51.9869,
		-51.6642,
		-50.5145,
		-50.098,
		-49.673,
		-51.0346,
		-48.4426,
		-46.532,
		-49.659,
		-49.8172,
		-47.0652,
		-48.479,
		-50.3125,
		-50.6249,
		-52.0916,
		-49.325,
		-46.4799,
		-48.3024,
		-51.8701,
		-48.8237,
		-50.4471,
		-50.5064,
		-49.4765,
		-51.0378,
		-49.8951,
		-50.867,
		-51.7528,
		-49.7907,
		-50.9414,
		-50.196,
		-50.7166,
		-48.2638,
		-48.1027,
		-49.4721,
		-51.5115,
		-49.6891,
		-49.1679,
		-50.4271,
		-50.3478,
		-50.5238,
		-49.163,
		-48.1769,
		-48.8715,
		-51.5548,
		-48.3888,
		-47.5323,
		-50.1061,
		-49.1536,
		-49.2668,
		-49.7307,
		-51.1017,
		-54.1429,
		-50.4325,
		-49.1318,
		-48.6643,
		-50.4365,
		-49.6167,
		-48.396,
		-49.2512,
		-50.9879,
		-49.5467,
		-50.9555,
		-53.2533,
		-50.8848,
		-50.4579,
		-50.1226,
		-49.8508,
		-49.3174,
		-50.3957,
		-47.6939,
		-47.1738,
		-49.8836,
		-49.2091,
		-49.4535,
		-53.6354,
		-52.5986,
		-51.6961,
		-52.7263,
		-48.666,
		-49.2609,
		-52.923,
		-52.6971,
		-51.1352,
		-50.2645,
		-48.4548,
		-48.1124,
		-47.8769,
		-48.4172,
		-46.8714,
		-48.2583,
		-50.2179,
		-48.3703,
		-49.6104,
		-49.8433,
		-46.8929,
		-47.5154,
		-50.6053,
		-51.3085,
		-51.9856,
		-50.2469,
		-47.5982,
		-49.3081,
		-51.3249,
		-48.7754,
		-49.3255,
		-50.9454,
		-49.4825,
		-50.0666,
		-49.8912,
		-48.8952,
		-48.1874,
		-49.0473,
		-48.5272,
		-47.6251,
		-50.1978,
		-51.7487,
		-48.6737,
		-49.9187,
		-50.7832,
		-48.591,
		-47.5307,
		-53.948,
		-49.7888,
		-47.7759,
		-49.4646,
		-49.2951,
		-48.0835,
		-50.9995,
		-47.7416,
		-47.4029,
		-51.7832,
		-47.9682,
		-46.4668,
		-49.6323,
		-50.7472,
		-48.5351,
		-48.8773,
		-49.049,
		-48.3435,
		-50.9687,
		-49.1747,
		-46.6598,
		-48.7942,
		-49.6085,
		-47.7136,
		-46.3717,
		-47.9584,
		-48.8272,
		-47.049,
		-48.0484,
		-48.5147,
		-47.7235,
		-48.9985,
		-48.3041,
		-48.7325,
		-52.1014,
		-48.6559,
		-45.9734,
		-47.9816,
		-49.9001,
		-49.744,
		-49.0705,
		-49.5089,
		-50.6053,
		-51.4918,
		-50.2007,
		-49.1247,
		-50.8951,
		-53.99,
		-50.9961,
		-51.2368,
		-54.649,
		-50.5483,
		-49.6662,
		-49.7394,
		-48.4737,
		-50.7644,
		-52.0322,
		-51.6659,
		-48.8891,
		-50.7512,
		-51.6192,
		-49.0519,
		-48.6193,
		-49.1657,
		-49.9413,
		-50.675,
		-50.7688,
		-47.5235,
		-48.7736,
		-51.2266,
		-50.0237,
		-48.7437,
		-51.924,
		-52.801,
		-49.2737,
		-46.3321,
		-48.3736,
		-48.1676,
		-46.3604,
		-48.1548,
		-51.7357,
		-48.9502,
		-48.1962,
		-48.5177,
		-49.363,
		-48.0272,
		-45.426,
		-48.0143,
		-48.0975,
		-45.1166,
		-46.3444,
		-48.1079,
		-47.5219,
		-47.5311,
		-47.5127,
		-45.4789,
		-47.1243,
		-48.7736,
		-46.554,
		-48.3435,
		-50.5012,
		-46.9389,
		-47.8544,
		-48.8036,
		-50.0142,
		-50.8305,
		-51.3919,
		-50.5138,
		-47.5832,
		-47.4375,
		-47.9406,
		-48.4136,
		-47.9574,
		-46.4125,
		-46.3805,
		-49.6796,
		-48.7645,
		-46.8555,
		-48.4917,
		-48.5139,
		-45.8423,
		-47.3436,
		-49.4883,
		-47.0694,
		-46.6695,
		-48.0118,
		-47.1591,
		-49.1952,
		-50.3417,
		-48.5627,
		-47.191,
		-48.3246,
		-48.2971,
		-47.0113,
		-48.5018,
		-49.5539,
		-48.1113,
		-48.9499,
		-48.0196,
		-44.7132,
		-46.8052,
		-49.3817,
		-48.9602,
		-51.7208,
		-47.954,
		-45.4842,
		-49.5872,
		-47.2174,
		-44.829,
		-46.5841,
		-48.7515,
		-47.0837,
		-47.6698,
		-49.6554,
		-49.0501,
		-48.0787,
		-47.2196,
		-48.6574,
		-46.4617,
		-47.7125,
		-47.2467,
		-49.2226,
		-50.1231,
		-47.5745,
		-45.2543,
		-46.22,
		-46.6847,
		-45.2459,
		-45.3651,
		-48.327,
		-45.4815,
		-45.0398,
		-45.7934,
		-45.2835,
		-44.7138,
		-45.9321,
		-45.2813,
		-47.0122,
		-46.104,
		-45.367,
		-45.9876,
		-47.9313,
		-48.014,
		-45.8913,
		-45.8209,
		-44.539,
		-43.8343,
		-47.3054,
		-48.6662,
		-47.7547,
		-46.8956,
		-48.9327,
		-47.4592,
		-48.1918,
		-47.9374,
		-45.3557,
		-45.4929,
		-45.3678,
		-43.5012,
		-43.1875,
		-45.2978,
		-46.5465,
		-49.1348,
		-49.1225,
		-46.0337,
		-46.0285,
		-47.0877,
		-44.4196,
		-44.1308,
		-46.6495,
		-45.2522,
		-44.6563,
		-45.8002,
		-44.6618,
		-44.5779,
		-44.4855,
		-44.2913,
		-43.1099,
		-43.9997,
		-45.2453,
		-45.0873,
		-46.2542,
		-47.2774,
		-45.2886,
		-44.5012,
		-45.9774,
		-44.1381,
		-45.4811,
		-48.5527,
		-47.2975,
		-44.8606,
		-46.5022,
		-46.7482,
		-46.028,
		-48.0085,
		-47.6895,
		-45.0938,
		-46.893,
		-47.9799,
		-46.101,
		-46.2139,
		-48.2228,
		-47.5895,
		-45.3641,
		-45.9702,
		-45.4339,
		-47.2054,
		-50.1577,
		-47.9846,
		-48.8655,
		-48.461,
		-48.6598,
		-48.6151,
		-49.7885,
		-47.47,
		-49.6225,
		-47.6268,
		-47.3095,
		-48.6953,
		-47.3902,
		-47.3501,
		-46.8645,
		-49.2446,
		-46.3305,
		-44.4973,
		-46.9743,
		-48.1871,
		-48.1099,
		-46.9727,
		-45.9809,
		-47.0705,
		-48.4405,
		-45.9353,
		-45.4907,
		-46.1851,
		-44.4699,
		-47.7013,
		-46.1764,
		-46.9609,
		-46.7163,
		-45.2347,
		-45.2775,
		-45.7384,
		-45.4096,
		-46.1625,
		-43.9906,
		-42.2237,
		-43.7476,
		-44.4925,
		-42.0441,
		-43.5432,
		-44.4195,
		-44.1445,
		-45.2037,
		-43.0051,
		-41.7508,
		-43.4584,
		-45.7652,
		-43.0487,
		-44.2161,
		-45.7777,
		-44.9954,
		-45.5829,
		-45.2463,
		-46.069,
		-46.7692,
		-44.4243,
		-45.4386,
		-45.8462,
		-45.5292,
		-42.6122,
		-42.72,
		-45.0977,
		-47.1575,
		-41.8999,
		-42.8585,
		-44.6372,
		-45.2421,
		-45.8997,
		-43.8037,
		-43.0463,
		-43.7606,
		-49.1001,
		-43.6061,
		-42.4586,
		-43.3683,
		-43.391,
		-44.8932,
		-46.0093,
		-46.4953,
		-46.0914,
		-45.0499,
		-44.6887,
		-44.5065,
		-44.6259,
		-43.2564,
		-41.6824,
		-44.2345,
		-46.3672,
		-44.4248,
		-44.7575,
		-45.4479,
		-46.4974,
		-46.4438,
		-48.1768,
		-47.8746,
		-48.1585,
		-46.0249,
		-44.2233,
		-45.0271,
		-47.7074,
		-45.1241,
		-45.602,
		-48.034,
		-45.7877,
		-45.7414,
		-48.8029,
		-46.1515,
		-44.7503,
		-46.0245,
		-48.1391,
		-42.7574,
		-45.0704,
		-44.4313,
		-42.3392,
		-42.7127,
		-42.63,
		-42.1699,
		-43.1941,
		-44.8777,
		-43.7871,
		-44.7024,
		-43.1134,
		-41.052,
		-41.5106,
		-42.7823,
		-43.9133,
		-46.6444,
		-44.1329,
		-42.4254,
		-41.9902,
		-43.1927,
		-42.4992,
		-43.2412,
		-45.32,
		-43.1171,
		-41.8517,
		-45.7062,
		-44.5682,
		-43.9681,
		-42.4479,
		-42.7704,
		-45.2259,
		-46.2546,
		-43.2546,
		-42.3056,
		-42.9318,
		-41.2086,
		-40.1974,
		-38.86,
		-41.7293,
		-45.3347,
		-43.45,
		-43.1411,
		-42.5701,
		-42.2307,
		-42.1508,
		-40.1295,
		-39.5435,
		-44.5262,
		-44.3852,
		-41.6828,
		-42.4575,
		-41.5466,
		-41.2296,
		-40.6723,
		-40.7443,
		-41.0065,
		-44.0477,
		-44.4363,
		-41.7784,
		-43.8286,
		-44.3334,
		-40.4965,
		-40.7184,
		-42.1522,
		-40.372,
		-40.0213,
		-42.1974,
		-43.9423,
		-41.4528,
		-41.3604,
		-40.0896,
		-40.5994,
		-43.1084,
		-44.5182,
		-40.5199,
		-43.2448,
		-43.9948,
		-42.4505,
		-43.8804,
		-41.658,
		-41.7391,
		-43.9621,
		-42.5052,
		-41.941,
		-43.953,
		-44.4372,
		-42.9089,
		-41.4266,
		-44.5048,
		-44.0139,
		-43.1635,
		-42.9775,
		-42.5989,
		-44.2387,
		-43.6555,
		-42.2137,
		-40.8761,
		-41.2583,
		-41.6775,
		-39.2895,
		-38.7132,
		-39.5674,
		-41.6928,
		-38.4184,
		-36.5258,
		-36.6337,
		-36.45,
		-35.4038,
		-34.3154,
		-33.8786,
		-34.3,
		-32.2024,
		-31.6128,
		-29.9677,
		-28.604,
		-26.7203,
		-21.491,
		-20.0435,
		-20.0388,
		-19.4482,
		-19.675,
		-15.5392,
		-15.7468,
		-15.7612,
		-14.1639,
		-17.6221,
		-20.3647,
		-15.2342,
		-14.4646,
		-15.9993,
		-16.7893,
		-20.0614,
		-18.6413,
		-15.156,
		-15.7714,
		-18.9955,
		-11.9265,
		-13.1928,
		-17.4033,
		-15.9239,
		-15.6211,
		-15.1626,
		-17.2121,
		-15.1454,
		-14.846,
		-17.4043,
		-14.9008,
		-17.8761,
		-13.6937,
		-12.7696,
		-16.4788,
		-17.8299,
		-14.1835,
		-15.9506,
		-15.5977,
		-16.4307,
		-16.6612,
		-17.2893,
		-17.1491,
		-15.2785,
		-14.299,
		-11.5101,
		-12.0057,
		-14.3695,
		-12.3953,
		-14.5376,
		-18.6463,
		-17.6179,
		-17.2854,
		-14.838,
		-15.901,
		-15.4921,
		-16.3385,
		-14.7699,
		-11.5637,
		-13.0916,
		-16.2336,
		-16.4366,
		-15.3438,
		-16.4631,
		-16.434,
		-15.4712,
		-14.356,
		-16.113,
		-16.3932,
		-14.813,
		-10.4706,
		-14.0698,
		-18.1615,
		-16.0609,
		-15.0685,
		-18.28,
		-15.6291,
		-15.2425,
		-18.0123,
		-16.8178,
		-13.9477,
		-18.1202,
		-20.0955,
		-18.5818,
		-18.6161,
		-17.0272,
		-15.1053,
		-14.8645,
		-14.2902,
		-14.2992,
		-18.4018,
		-18.3067,
		-14.4854,
		-17.361,
		-16.466,
		-12.033,
		-12.1618,
		-18.0388,
		-14.6222,
		-14.4357,
		-16.3651,
		-14.5199,
		-16.5932,
		-18.0895,
		-18.2942,
		-14.4461,
		-14.9826,
		-13.7885,
		-11.4138,
		-14.7876,
		-17.2821,
		-15.1468,
		-14.2192,
		-14.4969,
		-14.6106,
		-19.7936,
		-18.5169,
		-15.4286,
		-17.5611,
		-17.1634,
		-13.5942,
		-13.6943,
		-16.9909,
		-16.9429,
		-16.4109,
		-18.8415,
		-16.5409,
		-15.0603,
		-15.8583,
		-15.0508,
		-14.7035,
		-20.1458,
		-14.2932,
		-11.0877,
		-14.5908,
		-18.6891,
		-16.1547,
		-15.6604,
		-17.4981,
		-15.1965,
		-16.2621,
		-15.3162,
		-15.8825,
		-18.2088,
		-17.8679,
		-12.9174,
		-13.0332,
		-14.5191,
		-15.1663,
		-17.884,
		-19.5843,
		-16.3794,
		-15.5378,
		-15.8095,
		-16.8979,
		-16.1351,
		-16.554,
		-14.7715,
		-11.7863,
		-15.7083,
		-16.0304,
		-16.8404,
		-19.9122,
		-15.3532,
		-17.1025,
		-16.169,
		-18.5948,
		-18.5845,
		-17.9948,
		-20.0339,
		-16.0732,
		-15.9746,
		-18.4749,
		-17.5057,
		-15.7567,
		-18.9827,
		-18.5605,
		-19.4898,
		-19.8482,
		-19.9323,
		-19.4093,
		-20.9551,
		-19.3766,
		-18.472,
		-16.8604,
		-15.9094,
		-16.6785,
		-18.2195,
		-20.4397,
		-17.7038,
		-16.4312,
		-20.627,
		-20.5637,
		-17.8994,
		-19.3752,
		-15.7811,
		-15.8661,
		-21.2333,
		-17.3369,
		-17.6729,
		-17.944,
		-16.6844,
		-15.3104,
		-18.4577,
		-17.4173,
		-14.2345,
		-16.6316,
		-15.1417,
		-12.4047,
		-16.2736,
		-18.1997,
		-13.7877,
		-14.0002,
		-15.3601,
		-15.24,
		-16.6319,
		-17.412,
		-15.6287,
		-16.9564,
		-18.5734,
		-12.4446,
		-13.2975,
		-16.535,
		-16.7369,
		-16.7826,
		-16.7616,
		-17.2345,
		-15.6208,
		-15.6083,
		-14.7413,
		-15.5053,
		-19.3031,
		-14.9355,
		-12.4761,
		-16.9753,
		-19.1808,
		-14.589,
		-15.7518,
		-17.9452,
		-17.4412,
		-18.4867,
		-19.3394,
		-15.8374,
		-15.3029,
		-16.1687,
		-13.7927,
		-12.9113,
		-14.1068,
		-14.4755,
		-16.5459,
		-20.7771,
		-18.5065,
		-16.205,
		-14.6572,
		-14.6384,
		-15.8938,
		-17.6287,
		-14.2392,
		-12.21,
		-14.5071,
		-18.1666,
		-16.0671,
		-14.9861,
		-16.4561,
		-18.072,
		-15.2579,
		-13.6526,
		-14.9221,
		-15.7395,
		-13.5946,
		-10.9806,
		-14.4211,
		-17.8838,
		-15.2315,
		-16.0815,
		-17.462,
		-16.6821,
		-15.6187,
		-19.4224,
		-14.3781,
		-13.0367,
		-18.8141,
		-19.094,
		-17.2449,
		-19.502,
		-18.4997,
		-16.9683,
		-18.2979,
		-16.559,
		-14.2958,
		-17.2245,
		-16.5611,
		-15.965,
		-16.2773,
		-16.0848,
		-13.4482,
		-13.628,
		-17.3424,
		-15.2759,
		-15.4864,
		-15.5957,
		-15.4528,
		-17.081,
		-18.0077,
		-16.725,
		-16.0681,
		-16.6705,
		-14.7472,
		-11.1078,
		-15.5698,
		-18.9295,
		-14.5944,
		-15.3533,
		-14.4935,
		-14.1263,
		-18.3201,
		-16.4641,
		-14.0456,
		-16.3689,
		-15.9044,
		-13.9083,
		-15.1376,
		-18.0881,
		-15.4433,
		-18.0392,
		-19.8115,
		-17.2612,
		-16.1393,
		-14.6346,
		-12.639,
		-13.3451,
		-20.3823,
		-15.5033,
		-11.6252,
		-17.2243,
		-19.6212,
		-15.889,
		-16.2187,
		-17.2933,
		-15.3106,
		-15.2818,
		-15.3002,
		-15.3482,
		-16.0188,
		-15.113,
		-13.0713,
		-12.5087,
		-15.7176,
		-17.3474,
		-16.7773,
		-19.4358,
		-16.6256,
		-15.8939,
		-17.3975,
		-16.3073,
		-15.8907,
		-17.3589,
		-15.003,
		-11.9169,
		-16.0895,
		-14.965,
		-15.6599,
		-18.6817,
		-17.8752,
		-18.2858,
		-17.1836,
		-18.8071,
		-17.1642,
		-16.1271,
		-16.7472,
		-19.1897,
		-18.4565,
		-20.2236,
		-18.0552,
		-18.4603,
		-20.0452,
		-18.7502,
		-17.8809,
		-22.238,
		-29.1337,
		-30.655,
		-30.9884,
		-31.0855,
		-26.6656,
		-30.7148,
		-31.103,
		-29.2949,
		-26.9712,
		-21.8321,
		-17.1371,
		-19.202,
		-20.7523,
		-18.956,
		-18.5722,
		-22.2649,
		-17.2581,
		-17.3095,
		-17.8291,
		-16.7546,
		-18.3583,
		-19.4196,
		-17.1012,
		-19.0665,
		-19.6574,
		-18.7763,
		-15.6548,
		-17.41,
		-14.4814,
		-12.649,
		-16.4868,
		-19.6048,
		-15.9429,
		-14.6092,
		-16.2591,
		-16.6707,
		-16.7363,
		-17.5889,
		-15.5971,
		-17.7072,
		-16.9493,
		-12.379,
		-14.5237,
		-17.3588,
		-18.0528,
		-18.7503,
		-17.4407,
		-17.773,
		-16.1955,
		-16.316,
		-17.5686,
		-16.1528,
		-17.5352,
		-15.1959,
		-12.5491,
		-16.2168,
		-17.951,
		-16.2831,
		-16.8028,
		-21.7191,
		-17.5696,
		-18.5079,
		-17.8364,
		-16.8529,
		-17.3568,
		-16.4887,
		-14.5074,
		-13.6032,
		-15.0561,
		-14.8592,
		-17.1061,
		-18.2881,
		-17.3267,
		-17.4816,
		-15.2771,
		-15.9633,
		-17.1671,
		-16.8132,
		-15.4977,
		-13.4716,
		-14.6176,
		-15.4932,
		-15.5906,
		-16.0732,
		-16.6312,
		-19.8594,
		-17.7166,
		-14.7387,
		-16.1423,
		-17.7272,
		-16.0497,
		-12.0825,
		-14.8194,
		-18.6578,
		-16.2647,
		-16.952,
		-21.209,
		-17.2636,
		-15.4765,
		-21.519,
		-18.9253,
		-15.6929,
		-18.4423,
		-22.0175,
		-20.0796,
		-19.7565,
		-18.4389,
		-17.4147,
		-16.1838,
		-15.7579,
		-15.7313,
		-17.9307,
		-17.6456,
		-15.472,
		-16.2165,
		-18.3629,
		-13.9219,
		-13.3908,
		-18.0883,
		-15.8902,
		-15.8117,
		-17.6531,
		-15.399,
		-17.1985,
		-20.5278,
		-18.2017,
		-16.0071,
		-16.6188,
		-15.3262,
		-12.4708,
		-16.0418,
		-18.4806,
		-15.2803,
		-14.1438,
		-15.4045,
		-17.4351,
		-21.2766,
		-18.956,
		-14.8291,
		-15.6915,
		-17.9646,
		-13.796,
		-14.2603,
		-17.9241,
		-17.1372,
		-20.2816,
		-18.5349,
		-16.7273,
		-15.9286,
		-17.0292,
		-15.6254,
		-16.2442,
		-21.4171,
		-15.2941,
		-11.6143,
		-16.4831,
		-19.1565,
		-17.3051,
		-16.2204,
		-16.7853,
		-14.3267,
		-16.5206,
		-17.4987,
		-15.8609,
		-16.1717,
		-18.7535,
		-13.7566,
		-13.515,
		-16.0579,
		-16.2419,
		-17.9404,
		-20.2519,
		-18.8312,
		-16.4888,
		-16.1905,
		-18.4802,
		-15.4975,
		-16.8566,
		-15.1004,
		-12.0728,
		-16.1559,
		-16.5879,
		-17.5269,
		-17.6385,
		-17.0088,
		-18.2969,
		-17.6654,
		-18.2338,
		-18.5142,
		-18.5781,
		-23.0018,
		-18.7451,
		-18.7929,
		-20.3262,
		-18.2918,
		-20.992,
		-20.8563,
		-21.0641,
		-19.8925,
		-20.8524,
		-19.6837,
		-20.3027,
		-20.8755,
		-18.7651,
		-19.745,
		-19.0357,
		-22.9735,
		-19.6601,
		-17.5248,
		-22.4799,
		-22.2368,
		-20.741,
		-21.4475,
		-19.1367,
		-20.0864,
		-20.6374,
		-17.0022,
		-17.0454,
		-20.2269,
		-20.2276,
		-18.1309,
		-18.3089,
		-19.0382,
		-17.7965,
		-17.7277,
		-19.0534,
		-16.6147,
		-16.9542,
		-15.08,
		-13.8843,
		-18.0411,
		-20.2552,
		-17.1335,
		-15.8885,
		-16.4164,
		-19.8657,
		-21.2539,
		-21.6384,
		-17.9591,
		-16.4973,
		-19.4645,
		-13.7794,
		-14.5086,
		-17.499,
		-18.7878,
		-22.569,
		-18.4554,
		-19.112,
		-17.7439,
		-16.5406,
		-19.7502,
		-17.2812,
		-19.5686,
		-16.5935,
		-14.6281,
		-16.7502,
		-19.8096,
		-16.3894,
		-17.3978,
		-17.9746,
		-16.3131,
		-16.844,
		-18.1621,
		-16.3371,
		-16.2464,
		-18.0401,
		-14.3329,
		-13.5673,
		-16.2798,
		-15.4207,
		-17.4112,
		-20.4657,
		-19.6506,
		-18.5436,
		-18.5071,
		-18.5784,
		-18.1828,
		-18.302,
		-15.7905,
		-14.5212,
		-15.8785,
		-17.87,
		-17.5357,
		-17.2979,
		-19.5085,
		-15.4981,
		-16.5418,
		-16.4873,
		-16.6953,
		-18.9113,
		-17.337,
		-12.8396,
		-15.5374,
		-18.1843,
		-18.6142,
		-17.9819,
		-20.3665,
		-18.874,
		-17.8146,
		-21.8009,
		-19.4206,
		-15.1749,
		-16.941,
		-21.3388,
		-19.2361,
		-18.5137,
		-19.5836,
		-17.8294,
		-18.0813,
		-16.2072,
		-15.2342,
		-17.3446,
		-17.9416,
		-17.0639,
		-19.0464,
		-18.1247,
		-14.0078,
		-14.1669,
		-17.2219,
		-16.7482,
		-15.523,
		-17.009,
		-17.2274,
		-17.2415,
		-20.7098,
		-22.7188,
		-15.6795,
		-16.1834,
		-15.076,
		-12.8413,
		-17.1034,
		-21.6745,
		-16.3876,
		-15.0064,
		-15.6645,
		-15.5831,
		-19.8088,
		-18.1582,
		-15.6864,
		-16.8537,
		-16.605,
		-14.3675,
		-14.9494,
		-18.2275,
		-17.9394,
		-17.9904,
		-21.1667,
		-18.6842,
		-17.2572,
		-16.2678,
		-15.2376,
		-15.7674,
		-20.0147,
		-15.5518,
		-12.6766,
		-16.4039,
		-15.9753,
		-16.8163,
		-16.7973,
		-17.6159,
		-17.4915,
		-17.6545,
		-17.1689,
		-16.1077,
		-18.1998,
		-18.469,
		-13.9715,
		-13.9777,
		-16.6906,
		-16.3311,
		-17.1617,
		-19.8228,
		-17.2765,
		-18.7073,
		-20.429,
		-19.0216,
		-16.8083,
		-20.0403,
		-18.0369,
		-14.1958,
		-18.8201,
		-19.6189,
		-18.9291,
		-21.8515,
		-21.4555,
		-19.4145,
		-21.1048,
		-24.6117,
		-30.1201,
		-30.2087,
		-32.4339,
		-34.0687,
		-34.3296,
		-35.9902,
		-35.9791,
		-36.9887,
		-38.0751,
		-37.0478,
		-38.1653,
		-39.9887,
		-41.1763,
		-40.8435,
		-40.5692,
		-42.2582,
		-43.4551,
		-43.1281,
		-41.4753,
		-41.3605,
		-44.7594,
		-48.0518,
		-45.8695,
		-44.5699,
		-44.9164,
		-46.7043,
		-46.9795,
		-45.987,
		-46.8277,
		-46.9546,
		-48.517,
		-48.9181,
		-48.5357,
		-47.1731,
		-46.1186,
		-46.5777,
		-49.208,
		-50.2139,
		-47.5279,
		-47.6892,
		-45.6928,
		-42.3576,
		-45.6602,
		-50.0009,
		-44.0409,
		-44.8954,
		-46.0335,
		-46.3665,
		-47.0359,
		-46.5136,
		-43.8663,
		-45.1994,
		-48.9513,
		-43.6399,
		-44.5852,
		-48.9535,
		-48.7395,
		-50.6469,
		-47.9366,
		-50.2694,
		-53.1197,
		-51.219,
		-48.9197,
		-49.6145,
		-49.3445,
		-46.2711,
		-46.8294,
		-48.7328,
		-48.689,
		-46.8953,
		-46.455,
		-47.4365,
		-49.5511,
		-47.2848,
		-48.2091,
		-48.48,
		-49.6225,
		-49.5303,
		-45.9227,
		-45.6378,
		-47.6079,
		-48.3332,
		-51.1252,
		-50.0722,
		-47.7211,
		-46.3724,
		-47.9191,
		-46.9733,
		-46.4178,
		-47.9862,
		-48.7368,
		-47.3079,
		-48.573,
		-47.8612,
		-47.7978,
		-47.9014,
		-47.465,
		-48.1392,
		-51.7731,
		-50.9328,
		-48.3029,
		-48.4842,
		-51.4898,
		-48.723,
		-48.3427,
		-50.284,
		-50.1632,
		-49.0694,
		-49.4466,
		-49.4901,
		-51.7125,
		-51.408,
		-49.0947,
		-46.7358,
		-48.8949,
		-50.5354,
		-51.6669,
		-48.0811,
		-46.9125,
		-47.2465,
		-48.6119,
		-48.0179,
		-47.1497,
		-49.5535,
		-49.6195,
		-51.0516,
		-50.4786,
		-49.7187,
		-47.9761,
		-49.0253,
		-50.6974,
		-49.7176,
		-48.4694,
		-49.1365,
		-48.0496,
		-48.3849,
		-49.6228,
		-48.9914,
		-46.2694,
		-47.6119,
		-49.1836,
		-50.8548,
		-50.1212,
		-49.1347,
		-47.36,
		-49.579,
		-47.6367,
		-48.0236,
		-50.94,
		-51.1781,
		-50.3151,
		-47.9266,
		-49.14,
		-48.8484,
		-47.8906,
		-51.3558,
		-51.0253,
		-51.8225,
		-52.3138,
		-50.3064,
		-50.2402,
		-50.6688,
		-47.7897,
		-48.0829,
		-50.0277,
		-49.8401,
		-46.2139,
		-48.474,
		-48.5406,
		-46.5788,
		-46.0183,
		-48.064,
		-47.0858,
		-48.4572,
		-51.3744,
		-49.071,
		-48.7529,
		-49.3962,
		-47.3005,
		-47.7479,
		-49.1209,
		-45.7782,
		-46.6839,
		-49.1914,
		-48.0888,
		-48.3552,
		-48.7767,
		-50.5968,
		-50.9137,
		-49.9867,
		-47.3503,
		-46.2074,
		-50.2556,
		-49.8635,
		-46.9128,
		-48.4134,
		-50.8621,
		-48.2746,
		-47.2561,
		-48.3277,
		-48.462,
		-49.5798,
		-51.1221,
		-49.8699,
		-49.8763,
		-52.7402,
		-51.8437,
		-50.9413,
		-49.4195,
		-48.6313,
		-48.7089,
		-48.26,
		-49.1006,
		-51.7097,
		-51.2647,
		-50.0222,
		-46.266,
		-47.6016,
		-51.4242,
		-51.2863,
		-50.4255,
		-50.777,
		-50.6262,
		-52.4123,
		-51.1417,
		-48.6119,
		-49.1773,
		-49.7257,
		-51.0008,
		-51.1137,
		-49.4485,
		-47.3073,
		-46.4001,
		-49.1763,
		-49.6873,
		-48.8188,
		-47.7588,
		-48.8807,
		-49.5315,
		-50.5,
		-49.1801,
		-46.205,
		-49.0411,
		-48.2937,
		-47.7029,
		-49.864,
		-49.8018,
		-51.6115,
		-50.3661,
		-49.6347,
		-48.2495,
		-51.1527,
		-50.0531,
		-50.046,
		-51.5537,
		-52.1198,
		-52.4282,
		-51.7337,
		-50.2507,
		-50.2803,
		-49.4137,
		-49.3984,
		-50.0768,
		-48.9663,
		-49.2549,
		-49.3198,
		-48.9314,
		-51.3619,
		-50.9843,
		-49.4912,
		-50.9405,
		-48.5417,
		-46.2278,
		-48.7857,
		-49.9857,
		-48.5979,
		-49.3343,
		-49.8103,
		-51.2711,
		-50.6175,
		-52.9629,
		-46.3871,
		-47.3597,
		-50.0049,
		-51.3863,
		-50.9092,
		-49.2139,
		-49.7167,
		-48.2345,
		-49.7085,
		-47.9198,
		-47.2056,
		-48.779,
		-49.5685,
		-47.9345,
		-48.9486,
		-48.4122,
		-46.8779,
		-48.9406,
		-47.3671,
		-49.5441,
		-50.6269,
		-48.1014,
		-47.5496,
		-50.8901,
		-50.5588,
		-50.3705,
		-49.6213,
		-50.9238,
		-51.3464,
		-51.3144,
		-51.2102,
		-51.0348,
		-48.9954,
		-50.0406,
		-51.2378,
		-48.9083,
		-50.871,
		-49.1628,
		-49.3872,
		-47.8336,
		-45.5733,
		-45.3397,
		-48.3252,
		-49.8814,
		-49.5617,
		-52.5514,
		-50.8642,
		-47.4364,
		-46.9441,
		-49.2649,
		-48.6666,
		-48.7245,
		-52.2812,
		-50.4663,
		-48.9467,
		-50.4195,
		-49.3005,
		-48.0456,
		-49.7586,
		-51.4669,
		-49.4233,
		-49.3954,
		-51.918,
		-52.5363,
		-53.245,
		-50.1197,
		-49.5715,
		-48.6225,
		-51.4814,
		-48.6501,
		-48.435,
		-49.0636,
		-46.99,
		-47.1601,
		-51.3511,
		-51.7712,
		-50.9188,
		-49.9227,
		-47.5559,
		-47.5108,
		-47.5403,
		-48.5864,
		-48.7457,
		-52.6943,
		-47.607,
		-45.5889,
		-46.5265,
		-47.079,
		-46.4308,
		-46.6699,
		-46.537,
		-46.3403,
		-50.2296,
		-49.3748,
		-48.8292,
		-48.3511,
		-49.7581,
		-48.8874,
		-49.5481,
		-51.5764,
		-53.0174,
		-51.1775,
		-49.355,
		-48.1056,
		-50.0939,
		-52.3718,
		-49.6871,
		-48.6252,
		-48.5005,
		-49.1256,
		-47.7784,
		-49.501,
		-49.6147,
		-50.4504,
		-52.848,
		-52.5858,
		-51.5435,
		-50.0252,
		-50.4699,
		-49.7307,
		-50.7538,
		-52.7033,
		-51.7882,
		-52.9639,
		-53.4579,
		-50.473,
		-50.4015,
		-52.9672,
		-52.2223,
		-51.081,
		-52.4087,
		-53.6466,
		-53.7278,
		-51.9755,
		-50.2002,
		-48.3145,
		-50.6177,
		-51.2229,
		-48.8809,
		-52.0095,
		-53.5939,
		-50.5413,
		-52.9807,
		-53.6694,
		-53.6789,
		-54.1299,
		-51.5615,
		-51.0408,
		-54.0533,
		-51.3977,
		-49.306,
		-51.6357,
		-53.7669,
		-55.5275,
		-54.4678,
		-51.1255,
		-50.2686,
		-48.2089,
		-50.3321,
		-48.6478,
		-47.8024,
		-51.0462,
		-51.0639,
		-49.8515,
		-52.1705,
		-50.1129,
		-49.5335,
		-51.425,
		-52.8038,
		-50.2589,
		-50.7481,
		-48.5987,
		-47.1937,
		-48.4121,
		-50.5755,
		-51.7702,
		-50.9743,
		-51.1422,
		-50.6376,
		-51.6163,
		-50.8467,
		-52.6728,
		-50.7194,
		-50.9636,
		-53.1127,
		-53.6698,
		-51.9594,
		-49.8107,
		-48.4892,
		-46.9902,
		-49.1573,
		-49.527,
		-49.2974,
		-53.1054,
		-51.2926,
		-50.9592,
		-50.8261,
		-48.1979,
		-48.3261,
		-49.4646,
		-48.531,
		-49.1327,
		-50.5197,
		-50.4211,
		-52.3129,
		-52.3579,
		-49.2716,
		-47.8853,
		-48.2628,
		-46.8968,
		-47.9587,
		-50.7874,
		-50.925,
		-49.1492,
		-52.0662,
		-51.2057,
		-49.733,
		-49.6531,
		-50.086,
		-51.8516,
		-52.064,
		-49.5581,
		-48.727,
		-52.7099,
		-51.1081,
		-48.84,
		-51.0247,
		-54.5521,
		-55.9998,
		-51.3357,
		-51.0436,
		-53.4228,
		-51.0122,
		-50.6261,
		-53.7108,
		-54.3617,
		-52.8832,
		-50.9914,
		-48.9337,
		-48.592,
		-49.4738,
		-50.0585,
		-50.9642,
		-52.1417,
		-52.0936,
		-52.8795,
		-52.6958,
		-49.9302,
		-47.6595,
		-52.1132,
		-51.3054,
		-50.0328,
		-50.9853,
		-50.9056,
		-52.5706,
		-52.3182,
		-49.0835,
		-48.2247,
		-52.7862,
		-51.8768,
		-49.5388,
		-51.7117,
		-53.3644,
		-52.5562,
		-53.1913,
		-53.2061,
		-52.064,
		-51.4112,
		-52.2875,
		-51.7046,
		-53.4237,
		-51.6789,
		-48.9602,
		-49.615,
		-54.3065,
		-52.7595,
		-51.7539,
		-52.4489,
		-51.0018,
		-48.8141,
		-49.9922,
		-50.1697,
		-51.5571,
		-55.0326,
		-49.2505,
		-47.7099,
		-50.4369,
		-49.4718,
		-47.7357,
		-49.81,
		-49.7595,
		-49.5549,
		-52.5567,
		-50.5915,
		-51.071,
		-50.3581,
		-48.8385,
		-47.9719,
		-50.531,
		-53.754,
		-54.5275,
		-51.6008,
		-50.5125,
		-50.2016,
		-51.6389,
		-52.0413,
		-52.9119,
		-52.7238,
		-51.7382,
		-51.692,
		-50.3941,
		-52.276,
		-53.055,
		-51.7107,
		-51.5907,
		-56.1152,
		-54.4785,
		-53.1961,
		-52.6983,
		-51.7241,
		-55.7689,
		-51.4337,
		-52.7248,
		-53.1615,
		-56.9085,
		-52.0995,
		-54.9093,
		-56.176,
		-55.7428,
		-55.3085,
		-55.296,
		-55.3061,
		-54.5521,
		-38.7349,
		-36.6404,
		-35.6112,
		-36.3573,
		-38.6825,
		-40.9994,
		-54.838,
		-55.4042,
		-55.9963,
		-56.5734,
		-55.3016,
		-52.379,
		-54.4772,
		-57.5687,
		-55.7754,
		-56.1422,
		-58.4314,
		-56.3081,
		-58.3597,
		-54.8967,
		-55.1681,
		-54.94,
		-56.9699,
		-54.5725,
		-53.3058,
		-55.4576,
		-57.0444,
		-54.5734,
		-55.4882,
		-58.8719,
		-57.8147,
		-57.161,
		-55.9259,
		-54.1616,
		-54.9256,
		-57.9064,
		-55.5892,
		-54.4938,
		-56.6013,
		-53.3273,
		-53.5378,
		-56.8796,
		-54.122,
		-54.2593,
		-56.193,
		-53.9324,
		-55.5943,
		-57.948,
		-55.9327,
		-55.7052,
		-54.8921,
		-55.0404,
		-54.565,
		-55.2447,
		-53.395,
		-54.9005,
		-56.3495,
		-57.184,
		-56.2635,
		-57.2707,
		-57.4239,
		-56.9224,
		-56.0371,
		-56.0452,
		-57.16,
		-58.3847,
		-58.6806,
		-56.2668,
		-58.4003,
		-56.7104,
		-56.2934,
		-57.4744,
		-56.0358,
		-54.3808,
		-54.3828,
		-54.9765,
		-55.1621,
		-55.7418,
		-55.4705,
		-58.0448,
		-56.9212,
		-57.1734,
		-55.65,
		-54.5887,
		-57.506,
		-60.9133,
		-57.8444,
		-57.463,
		-60.6389,
		-57.9056,
		-56.315,
		-57.7275,
		-59.5114,
		-58.9604,
		-57.4103,
		-56.0218,
		-55.2313,
		-55.5317,
		-57.8286,
		-55.9575,
		-57.8296,
		-58.639,
		-55.4288,
		-55.5388,
		-54.5277,
		-55.7242,
		-56.5171,
		-57.9004,
		-58.5774,
		-58.0525,
		-57.3392,
		-56.6642,
		-57.3167,
		-57.4948,
		-57.6944,
		-58.4729,
		-56.7523,
		-56.0368,
		-56.9025,
		-56.2476,
		-55.9178,
		-56.8256,
		-57.2025,
		-56.0975,
		-56.4654,
		-55.9343,
		-58.3103,
		-56.9723,
		-56.2622,
		-58.8845,
		-57.8157,
		-58.9594,
		-58.1085,
		-57.076,
		-58.3772,
		-61.9133,
		-59.8652,
		-58.0047,
		-60.263,
		-58.0358,
		-58.024,
		-61.0247,
		-59.0908,
		-58.1108,
		-58.3913,
		-57.9262,
		-58.2836,
		-59.0939,
		-58.8803,
		-58.1376,
		-57.91,
		-58.8683,
		-58.4782,
		-59.3043,
		-58.0344,
		-57.5405,
		-57.2232,
		-56.4142,
		-55.2308,
		-58.6367,
		-60.2894,
		-58.3026,
		-57.6476,
		-58.7072,
		-57.2483,
		-57.0139,
		-58.3794,
		-58.8506,
		-58.9216,
		-61.2998,
		-60.9501,
		-60.7788,
		-57.7742,
		-56.8946,
		-56.6763,
		-58.6238,
		-59.6427,
		-58.9964,
		-60.5465,
		-58.0314,
		-57.6388,
		-59.3954,
		-58.3974,
		-57.6562,
		-58.0348,
		-58.047,
		-60.6587,
		-60.8311,
		-57.935,
		-58.2454,
		-58.5662,
		-58.1951,
		-58.0662,
		-59.569,
		-59.0606,
		-58.507,
		-57.3905,
		-58.0025,
		-59.7966,
		-58.2268,
		-58.0629,
		-57.8326,
		-57.1805,
		-59.3286,
		-59.5327,
		-59.4585,
		-61.4747,
		-59.7711,
		-58.6101,
		-58.3887,
		-58.4864,
		-58.0447,
		-56.8286,
		-55.9196,
		-56.3552,
		-59.9929,
		-59.1182,
		-60.4975,
		-59.3041,
		-57.4506,
		-56.5174,
		-57.1703,
		-60.0005,
		-59.6632,
		-59.0451,
		-59.9379,
		-57.636,
		-56.5565,
		-58.4742,
		-58.3485,
		-55.38,
		-59.3836,
		-57.7059,
		-57.7295,
		-57.0076,
		-56.6843,
		-58.024,
		-58.3442,
		-58.6589,
		-58.0776,
		-57.2852,
		-59.8936,
		-58.368,
		-58.8886,
		-61.1485,
		-57.9197,
		-59.1586,
		-60.8125,
		-57.1896,
		-56.3304,
		-60.1898,
		-56.4633,
};
	double *out = (double *)malloc(sizeof(in));
	FILE *fp = fopen("smooth.txt","wb");

	printf("\n=================================\n");
	if(1)
	{
		//const double f = 0.25;
		double f = 0.01;
		const size_t nsteps = 3;
		const size_t n = sizeof(in)/sizeof(in[0]);
		const double delta = 0.3;

		double *ys  = (double *) malloc(sizeof(double)*n);
		double *rw  = (double *) malloc(sizeof(double)*n);
		double *res = (double *) malloc(sizeof(double)*n);
		const double *y = &in[0];
		double *x = (double *) malloc(sizeof(double)*n);
		
		double start_freq = -125000000;
		double rbw = 100000;
		for (int i = 0; i < n; i++)
		{
			x[i] = start_freq+(rbw*i);
		}

		lowess((const double*)x, (const double*)y, n,f, nsteps,delta, ys, rw, res);
		memcpy(out,in,sizeof(in));

		int center_point = 250;
		int smooth_count = 10;
		for (int i = 0; i < smooth_count; i++)
		{
			if (abs(abs(in[(center_point-5)+i]) - abs(ys[(center_point-5)+i])) > 3)
			{
				out[(center_point-5)+i] = ys[(center_point-5)+i];
			}
		}
		center_point = 2500 - 250;
		for (int i = 0; i < smooth_count; i++)
		{
			if (abs(abs(in[(center_point-5)+i]) - abs(ys[(center_point-5)+i])) > 3)
			{
				out[(center_point-5)+i] = ys[(center_point-5)+i];
			}
		}

		for (int i = 0; i < n; i++)
		{
			fprintf(fp,"%.4f,%.4f\r\n",in[i],out[i]);
		}
		fclose(fp);
		
		fp = fopen("smooth_2.csv","wb");
		const int column = 50;
		double *result[column] = {NULL};
		
		for (int i = 0; i < column; i++)
		{
			f = 0.001 + (i*0.005);
			result[i] = (double *) malloc(sizeof(double)*n);
			lowess((const double*)x, (const double*)y, n,f, nsteps,delta,result[i], rw, res);
		}
		for (int j = 0; j < n; j++)
		{
			for (int i = 0; i < column; i++)
			{
				fprintf(fp,"%.4f",result[i][j]);
				if (i != column-1)
				{
					fprintf(fp,",");
				}	
			}
			fprintf(fp,"\r\n");
		}
		for (int i = 0; i < column; i++)
		{
			free(result[i]);
		}

		free(x);
		free(ys);
		free(rw);
		free(res);
	}
	fclose(fp);
	printf("create ok\n");
	getchar();
}


版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/cupidove/article/details/52260432

智能推荐

win7 + centos7 双系统启动_centos7 windows 同时启动-程序员宅基地

文章浏览阅读1.2w次,点赞5次,收藏16次。本文介绍的是在已有windows系统(默认安装在C盘)基础上安装centos7,设置双系统启动。难点: 1、linux安装程序无法识别NTFS,windows系统无法读写ext3 2、U盘启动盘无法放入大于4G的文件 3、Gurb2 添加启动项分区准备1、下载分区助手,从原有分区中切割分区并删除分区作为安装linux的空间 2、因为linux安装程序无法识别NTFS,U盘启动盘(FAT32格_centos7 windows 同时启动

HTML文档模板 | html/template-程序员宅基地

文章浏览阅读1.7k次。import "html/template"概述索引示例概观模板包(html/template)实现了数据驱动的模板,以便在代码注入过程中安全地生成HTML输出。它提供了与包文本/模板相同的接口,只要输出是HTML,就应该使用它来代替文本/模板。这里的文档侧重于包的安全特性。有关如何自行编写模板的信息,请参阅文本/模板的文档。介绍该软件包包装文本/模板,以便您可以共享其模板API以安全地解析和执行..._html template

FastReport教程:如何在报表中使用多个数据库_fastreport 数据源-程序员宅基地

文章浏览阅读2.6k次。下载FastReport.Net最新版本有时,我们必须以不同的格式处理来自不同来源的数据。对于分析师和报表开发人员来说,这可能是一个令人头疼的问题。毕竟,你必须以某种方式组合数据。幸运的是,在FastReport.Net的报表中,您可以创建许多数据连接。而且,数据源可以完全不同 - 文本文件,数据库。多亏了这一点,我们将能够在一份报表中整合数据。 在本文中,我们将介绍在报表中创建两个数据源以及..._fastreport 数据源

学术海报Poster-- 模板分享_学术海报poster模板下载-程序员宅基地

文章浏览阅读4.4k次,点赞18次,收藏13次。读研期间,发表的论文被录用,一般会通过口述演讲或者Poster海报的形式向参与者展示你的论文科研成果,其中受众面积最大的一般是Poster海报分享的形式。对于论文录用者来说,它也是最简单的一种参会形式,而拥有一份精美的海报模板,对于广大的研究生来说,能省时省力不少,科研工作成果好很重要,但是,用精美的海报展示您的科研成果,让更多的读者了解到你的科研内容/成果,同样也非常重要。我在读研期间,就苦于寻找一份精美的海报模板而花费大量时间,现将这100份模板海报分享给你,希望对你的科研之路也能有一些帮助_学术海报poster模板下载

The k-th Largest Group POJ - 2985 treap+并查集_newman likes playing with cats. he possesses lots -程序员宅基地

文章浏览阅读306次。Newman likes playing with cats. He possesses lots of cats in his home. Because the number of cats is really huge, Newman wants to group some of the cats. To do that, he first offers a number to each o_newman likes playing with cats. he possesses lots of cats in his

C#向Excel报表中插入图片的2种方法_c#insertpictoexcel 类数据集-程序员宅基地

文章浏览阅读7.2k次。这几天做向Excel插入数据,其中有插入图片的需求,经试验,下面2种方法都可以插入图片,但各有不同的用处。现将这2种方法共享出来, 希望需要的朋友进行参考,代码中已经有详细注释了。注意:使用之前需要引用COM:Microsoft Office 11.0 Object Library如果引用列表中没有,需要自行添加 C:/Program Files/Microsoft Office/OFFIC..._c#insertpictoexcel 类数据集

随便推点

避免出现bitmap内存限制OUT OF MEMORY的一种方法-程序员宅基地

文章浏览阅读261次。在编写Android程序的时候,我们总是难免会碰到OOM(OUT OF MEMORY)的错误,那么这个错误究竟是怎么来的呢,可以先看一下这篇文章ANDROID BITMAP内存限制OOM,OUT OF MEMORY。 这里,我使用Gallery来举例,在模拟器中,不会出现OOM错误,但是,一旦把程序运行到真机里,图片文件一多,必然会出现OOM,我们通过做一些额外的处理来避免。1.创建一..._取消限制outofmemory限制

Bailian2966 时区转换【时区计算】-程序员宅基地

文章浏览阅读399次。2966:时区转换总时间限制: 1000ms 内存限制: 65536kB描述直到19世纪,时间校准是一个纯粹的地方现象。每一个村庄当太阳升到最高点的时候把他们的时钟调到中午12点。一个钟表制造商人家或者村里主表的时间被认为是官方时间,市民们把自家的钟表和这个时间对齐。每周一些热心的市民会带着时间标准的表,游走大街小巷为其他市民对表。在城市之间旅游的话,在到达新地方的时候需要把怀表校准。但是,..._2966

android之SharedPreferes_在android中sharedpreferences的文件可以跨类读取吗-程序员宅基地

文章浏览阅读710次。SharedPreferences是Android平台上一个轻量级的存储类,主要是保存一些常用的配置比如窗口状态,一般在Activity中 重载窗口状态onSaveInstanceState保存一般使用SharedPreferences完成,它提供了Android平台常规的Long长 整形、Int整形、String字符串型的保存,它是什么样的处理方式呢?SharedPreferences类似过去W_在android中sharedpreferences的文件可以跨类读取吗

win7和linux(centos)双系统安装成win7单系统_centoslinux能改win7系统吗?-程序员宅基地

文章浏览阅读453次。首先说明一下环境:安装了win7旗舰版和centos双系统。现在的需求是删除linux系统,然后只剩下win7系统,之后再把win7重新安装。步骤:1,先使用大白菜制作一个u盘启动工具,并放入下载好的win7系统镜像2,在win7的磁盘管理(计算机右键->管理->磁盘管理)3,在磁盘管理中删除linux的分区4,重新启动系统,这时有些系统就不能进入win7了5,_centoslinux能改win7系统吗?

layout_gravity和gravity的作用和区别(上)-程序员宅基地

文章浏览阅读2k次,点赞6次,收藏7次。layout_gravity和gravity的作用和区别(上)问题描述: 最近在学习Android布局时需要将位于LinearLayout布局中的控件放在布局底部,这时会使用到android:layout_gravity=“bottom”,但是发现这样是行不通的,最后查资料找出原因,防止忘记在这记一下。gravity: 是对view控件本身来说的,是用来设置view本身的内容应该显示在vi..._layout_gravity

[Perforce系列—] 1. Perforce 的使用和常用命令-程序员宅基地

文章浏览阅读2.8w次,点赞7次,收藏20次。常见使用 p4 的方式就是使用Client 端,但是有的时候遇到这样的状况:1. Client 端的操作方式的局限, 有些操作无法实现 (一般这样的状况不多)2. 使用P4 和其他一些工具进行整合, 比如p4 整合自动部署3. 使用代码的方式来使用p4, 比如使用Perl语言来与P4 进行交互。。。。以上的状况, 直接使用P4 Client 的话, 就没办法达成我们的要求了。这个时候自然就要使用到P4的命令行。_perforce

推荐文章

热门文章

相关标签