Return to previous page

Contents of file 'sgfilter/sgfilter.c':




Deprecated: Function split() is deprecated in /storage/content/45/2011745/jonsson.eu/public_html/php/htmlicise.php on line 46
    1   /*13:*/
    2   #line 975 "./sgfilter.w"
    3   
    4   /*14:*/
    5   #line 1002 "./sgfilter.w"
    6   
    7   #include <math.h> 
    8   #include <stdlib.h> 
    9   #include <stdarg.h> 
   10   #include <stdio.h> 
   11   #include <string.h> 
   12   #include <ctype.h> 
   13   #include <time.h> 
   14   #include <sys/time.h> 
   15   #include "sgfilter.h"
   16   
   17   /*:14*/
   18   #line 976 "./sgfilter.w"
   19   
   20   /*15:*/
   21   #line 1022 "./sgfilter.w"
   22   
   23   extern char*optarg;
   24   char*progname;
   25   
   26   /*:15*/
   27   #line 977 "./sgfilter.w"
   28   
   29   /*22:*/
   30   #line 1210 "./sgfilter.w"
   31   
   32   /*23:*/
   33   #line 1250 "./sgfilter.w"
   34   
   35   void log_printf(const char*function_name,const char*format,...){
   36   va_list args;
   37   time_t time0;
   38   struct tm lt;
   39   struct timeval tv;
   40   char logentry[1024];
   41   
   42   gettimeofday(&tv,NULL);
   43   time(&time0);
   44   lt= *localtime(&time0);
   45   sprintf(logentry,"%02u%02u%02u %02u:%02u:%02u.%03d ",
   46   lt.tm_year-100,lt.tm_mon+1,lt.tm_mday,
   47   lt.tm_hour,lt.tm_min,lt.tm_sec,tv.tv_usec/1000);
   48   sprintf(logentry+strlen(logentry),"(%s) ",function_name);
   49   va_start(args,format);
   50   vsprintf(logentry+strlen(logentry),format,args);
   51   va_end(args);
   52   sprintf(logentry+strlen(logentry),"\n");
   53   fprintf(stdout,"%s",logentry);
   54   return;
   55   }
   56   
   57   /*:23*/
   58   #line 1211 "./sgfilter.w"
   59   
   60   /*24:*/
   61   #line 1276 "./sgfilter.w"
   62   
   63   int*ivector(long nl,long nh){
   64   int*v;
   65   v= (int*)malloc((size_t)((nh-nl+2)*sizeof(int)));
   66   if(!v){
   67   log("Error: Allocation failure.");
   68   exit(1);
   69   }
   70   return v-nl+1;
   71   }
   72   
   73   /*:24*/
   74   #line 1212 "./sgfilter.w"
   75   
   76   /*25:*/
   77   #line 1290 "./sgfilter.w"
   78   
   79   double*dvector(long nl,long nh){
   80   double*v;
   81   long k;
   82   v= (double*)malloc((size_t)((nh-nl+2)*sizeof(double)));
   83   if(!v){
   84   log("Error: Allocation failure.");
   85   exit(1);
   86   }
   87   for(k= nl;k<=nh;k++)v[k]= 0.0;
   88   return v-nl+1;
   89   }
   90   
   91   /*:25*/
   92   #line 1213 "./sgfilter.w"
   93   
   94   /*26:*/
   95   #line 1307 "./sgfilter.w"
   96   
   97   double**dmatrix(long nrl,long nrh,long ncl,long nch)
   98   {
   99   long i,nrow= nrh-nrl+1,ncol= nch-ncl+1;
  100   double**m;
  101   m= (double**)malloc((size_t)((nrow+1)*sizeof(double*)));
  102   if(!m){
  103   log("Allocation failure 1 occurred.");
  104   exit(1);
  105   }
  106   m+= 1;
  107   m-= nrl;
  108   m[nrl]= (double*)malloc((size_t)((nrow*ncol+1)*sizeof(double)));
  109   if(!m[nrl]){
  110   log("Allocation failure 2 occurred.");
  111   exit(1);
  112   }
  113   m[nrl]+= 1;
  114   m[nrl]-= ncl;
  115   for(i= nrl+1;i<=nrh;i++)m[i]= m[i-1]+ncol;
  116   return m;
  117   }
  118   
  119   /*:26*/
  120   #line 1214 "./sgfilter.w"
  121   
  122   /*27:*/
  123   #line 1333 "./sgfilter.w"
  124   
  125   void free_ivector(int*v,long nl,long nh){
  126   free((char*)(v+nl-1));
  127   }
  128   
  129   /*:27*/
  130   #line 1215 "./sgfilter.w"
  131   
  132   /*28:*/
  133   #line 1341 "./sgfilter.w"
  134   
  135   void free_dvector(double*v,long nl,long nh){
  136   free((char*)(v+nl-1));
  137   }
  138   
  139   /*:28*/
  140   #line 1216 "./sgfilter.w"
  141   
  142   /*29:*/
  143   #line 1350 "./sgfilter.w"
  144   
  145   void free_dmatrix(double**m,long nrl,long nrh,long ncl,long nch){
  146   free((char*)(m[nrl]+ncl-1));
  147   free((char*)(m+nrl-1));
  148   }
  149   
  150   /*:29*/
  151   #line 1217 "./sgfilter.w"
  152   
  153   /*30:*/
  154   #line 1372 "./sgfilter.w"
  155   
  156   void lubksb(double**a,int n,int*indx,double b[])
  157   {
  158   int i,ii= 0,ip,j;
  159   double sum;
  160   
  161   for(i= 1;i<=n;i++){
  162   ip= indx[i];
  163   sum= b[ip];
  164   b[ip]= b[i];
  165   if(ii)
  166   for(j= ii;j<=i-1;j++)sum-= a[i][j]*b[j];
  167   else if(sum)ii= i;
  168   b[i]= sum;
  169   }
  170   for(i= n;i>=1;i--){
  171   sum= b[i];
  172   for(j= i+1;j<=n;j++)sum-= a[i][j]*b[j];
  173   b[i]= sum/a[i][i];
  174   }
  175   }
  176   
  177   /*:30*/
  178   #line 1218 "./sgfilter.w"
  179   
  180   /*31:*/
  181   #line 1407 "./sgfilter.w"
  182   
  183   void ludcmp(double**a,int n,int*indx,double*d)
  184   {
  185   int i,imax= 0,j,k;
  186   double big,dum,sum,temp;
  187   double*vv;
  188   
  189   vv= dvector(1,n);
  190   *d= 1.0;
  191   for(i= 1;i<=n;i++){
  192   big= 0.0;
  193   for(j= 1;j<=n;j++)
  194   if((temp= fabs(a[i][j]))> big)big= temp;
  195   if(big==0.0){
  196   log("Error: Singular matrix found in routine ludcmp()");
  197   exit(1);
  198   }
  199   vv[i]= 1.0/big;
  200   }
  201   for(j= 1;j<=n;j++){
  202   for(i= 1;i<j;i++){
  203   sum= a[i][j];
  204   for(k= 1;k<i;k++)sum-= a[i][k]*a[k][j];
  205   a[i][j]= sum;
  206   }
  207   big= 0.0;
  208   for(i= j;i<=n;i++){
  209   sum= a[i][j];
  210   for(k= 1;k<j;k++)
  211   sum-= a[i][k]*a[k][j];
  212   a[i][j]= sum;
  213   if((dum= vv[i]*fabs(sum))>=big){
  214   big= dum;
  215   imax= i;
  216   }
  217   }
  218   if(j!=imax){
  219   for(k= 1;k<=n;k++){
  220   dum= a[imax][k];
  221   a[imax][k]= a[j][k];
  222   a[j][k]= dum;
  223   }
  224   *d= -(*d);
  225   vv[imax]= vv[j];
  226   }
  227   indx[j]= imax;
  228   if(a[j][j]==0.0)a[j][j]= EPSILON;
  229   if(j!=n){
  230   dum= 1.0/(a[j][j]);
  231   for(i= j+1;i<=n;i++)a[i][j]*= dum;
  232   }
  233   }
  234   free_dvector(vv,1,n);
  235   }
  236   
  237   /*:31*/
  238   #line 1219 "./sgfilter.w"
  239   
  240   /*32:*/
  241   #line 1471 "./sgfilter.w"
  242   
  243   #include <math.h> 
  244   #define SWAP(a,b) tempr= (a);(a)= (b);(b)= tempr
  245   void four1(double data[],unsigned long nn,int isign)
  246   {
  247   unsigned long n,mmax,m,j,istep,i;
  248   double wtemp,wr,wpr,wpi,wi,theta;
  249   double tempr,tempi;
  250   
  251   n= nn<<1;
  252   j= 1;
  253   for(i= 1;i<n;i+= 2){
  254   if(j> i){
  255   SWAP(data[j],data[i]);
  256   SWAP(data[j+1],data[i+1]);
  257   }
  258   m= n>>1;
  259   while(m>=2&&j> m){
  260   j-= m;
  261   m>>= 1;
  262   }
  263   j+= m;
  264   }
  265   mmax= 2;
  266   while(n> mmax){
  267   istep= mmax<<1;
  268   theta= isign*(6.28318530717959/mmax);
  269   wtemp= sin(0.5*theta);
  270   wpr= -2.0*wtemp*wtemp;
  271   wpi= sin(theta);
  272   wr= 1.0;
  273   wi= 0.0;
  274   for(m= 1;m<mmax;m+= 2){
  275   for(i= m;i<=n;i+= istep){
  276   j= i+mmax;
  277   tempr= wr*data[j]-wi*data[j+1];
  278   tempi= wr*data[j+1]+wi*data[j];
  279   data[j]= data[i]-tempr;
  280   data[j+1]= data[i+1]-tempi;
  281   data[i]+= tempr;
  282   data[i+1]+= tempi;
  283   }
  284   wr= (wtemp= wr)*wpr-wi*wpi+wr;
  285   wi= wi*wpr+wtemp*wpi+wi;
  286   }
  287   mmax= istep;
  288   }
  289   }
  290   #undef SWAP
  291   
  292   
  293   /*:32*/
  294   #line 1220 "./sgfilter.w"
  295   
  296   /*33:*/
  297   #line 1532 "./sgfilter.w"
  298   
  299   void twofft(double data1[],double data2[],double fft1[],double fft2[],
  300   unsigned long n)
  301   {
  302   void four1(double data[],unsigned long nn,int isign);
  303   unsigned long nn3,nn2,jj,j;
  304   double rep,rem,aip,aim;
  305   
  306   nn3= 1+(nn2= 2+n+n);
  307   for(j= 1,jj= 2;j<=n;j++,jj+= 2){
  308   fft1[jj-1]= data1[j];
  309   fft1[jj]= data2[j];
  310   }
  311   four1(fft1,n,1);
  312   fft2[1]= fft1[2];
  313   fft1[2]= fft2[2]= 0.0;
  314   for(j= 3;j<=n+1;j+= 2){
  315   rep= 0.5*(fft1[j]+fft1[nn2-j]);
  316   rem= 0.5*(fft1[j]-fft1[nn2-j]);
  317   aip= 0.5*(fft1[j+1]+fft1[nn3-j]);
  318   aim= 0.5*(fft1[j+1]-fft1[nn3-j]);
  319   fft1[j]= rep;
  320   fft1[j+1]= aim;
  321   fft1[nn2-j]= rep;
  322   fft1[nn3-j]= -aim;
  323   fft2[j]= aip;
  324   fft2[j+1]= -rem;
  325   fft2[nn2-j]= aip;
  326   fft2[nn3-j]= rem;
  327   }
  328   }
  329   
  330   /*:33*/
  331   #line 1221 "./sgfilter.w"
  332   
  333   /*34:*/
  334   #line 1576 "./sgfilter.w"
  335   
  336   void realft(double data[],unsigned long n,int isign)
  337   {
  338   void four1(double data[],unsigned long nn,int isign);
  339   unsigned long i,i1,i2,i3,i4,np3;
  340   double c1= 0.5,c2,h1r,h1i,h2r,h2i;
  341   double wr,wi,wpr,wpi,wtemp,theta;
  342   
  343   theta= 3.141592653589793/(double)(n>>1);
  344   if(isign==1){
  345   c2= -0.5;
  346   four1(data,n>>1,1);
  347   }else{
  348   c2= 0.5;
  349   theta= -theta;
  350   }
  351   wtemp= sin(0.5*theta);
  352   wpr= -2.0*wtemp*wtemp;
  353   wpi= sin(theta);
  354   wr= 1.0+wpr;
  355   wi= wpi;
  356   np3= n+3;
  357   for(i= 2;i<=(n>>2);i++){
  358   i4= 1+(i3= np3-(i2= 1+(i1= i+i-1)));
  359   h1r= c1*(data[i1]+data[i3]);
  360   h1i= c1*(data[i2]-data[i4]);
  361   h2r= -c2*(data[i2]+data[i4]);
  362   h2i= c2*(data[i1]-data[i3]);
  363   data[i1]= h1r+wr*h2r-wi*h2i;
  364   data[i2]= h1i+wr*h2i+wi*h2r;
  365   data[i3]= h1r-wr*h2r+wi*h2i;
  366   data[i4]= -h1i+wr*h2i+wi*h2r;
  367   wr= (wtemp= wr)*wpr-wi*wpi+wr;
  368   wi= wi*wpr+wtemp*wpi+wi;
  369   }
  370   if(isign==1){
  371   data[1]= (h1r= data[1])+data[2];
  372   data[2]= h1r-data[2];
  373   }else{
  374   data[1]= c1*((h1r= data[1])+data[2]);
  375   data[2]= c1*(h1r-data[2]);
  376   four1(data,n>>1,-1);
  377   }
  378   }
  379   
  380   /*:34*/
  381   #line 1222 "./sgfilter.w"
  382   
  383   /*35:*/
  384   #line 1638 "./sgfilter.w"
  385   
  386   char convlv(double data[],unsigned long n,double respns[],unsigned long m,
  387   int isign,double ans[])
  388   {
  389   void realft(double data[],unsigned long n,int isign);
  390   void twofft(double data1[],double data2[],double fft1[],double fft2[],
  391   unsigned long n);
  392   unsigned long i,no2;
  393   double dum,mag2,*fft;
  394   
  395   fft= dvector(1,n<<1);
  396   for(i= 1;i<=(m-1)/2;i++)
  397   respns[n+1-i]= respns[m+1-i];
  398   for(i= (m+3)/2;i<=n-(m-1)/2;i++)
  399   respns[i]= 0.0;
  400   twofft(data,respns,fft,ans,n);
  401   no2= n>>1;
  402   for(i= 2;i<=n+2;i+= 2){
  403   if(isign==1){
  404   ans[i-1]= (fft[i-1]*(dum= ans[i-1])-fft[i]*ans[i])/no2;
  405   ans[i]= (fft[i]*dum+fft[i-1]*ans[i])/no2;
  406   }else if(isign==-1){
  407   if((mag2= ans[i-1]*ans[i-1]+ans[i]*ans[i])==0.0){
  408   log("Attempt of deconvolving at zero response in convlv().");
  409   return(1);
  410   }
  411   ans[i-1]= (fft[i-1]*(dum= ans[i-1])+fft[i]*ans[i])/mag2/no2;
  412   ans[i]= (fft[i]*dum-fft[i-1]*ans[i])/mag2/no2;
  413   }else{
  414   log("No meaning for isign in convlv().");
  415   return(1);
  416   }
  417   }
  418   ans[2]= ans[n+1];
  419   realft(ans,n,-1);
  420   free_dvector(fft,1,n<<1);
  421   return(0);
  422   }
  423   
  424   /*:35*/
  425   #line 1223 "./sgfilter.w"
  426   
  427   /*36:*/
  428   #line 1695 "./sgfilter.w"
  429   
  430   char sgcoeff(double c[],int np,int nl,int nr,int ld,int m)
  431   {
  432   void lubksb(double**a,int n,int*indx,double b[]);
  433   void ludcmp(double**a,int n,int*indx,double*d);
  434   int imj,ipj,j,k,kk,mm,*indx;
  435   double d,fac,sum,**a,*b;
  436   
  437   if(np<nl+nr+1||nl<0||nr<0||ld> m||nl+nr<m){
  438   log("Inconsistent arguments detected in routine sgcoeff.");
  439   return(1);
  440   }
  441   indx= ivector(1,m+1);
  442   a= dmatrix(1,m+1,1,m+1);
  443   b= dvector(1,m+1);
  444   for(ipj= 0;ipj<=(m<<1);ipj++){
  445   sum= (ipj?0.0:1.0);
  446   for(k= 1;k<=nr;k++)sum+= pow((double)k,(double)ipj);
  447   for(k= 1;k<=nl;k++)sum+= pow((double)-k,(double)ipj);
  448   mm= (ipj<2*m-ipj?ipj:2*m-ipj);
  449   for(imj= -mm;imj<=mm;imj+= 2)a[1+(ipj+imj)/2][1+(ipj-imj)/2]= sum;
  450   }
  451   ludcmp(a,m+1,indx,&d);
  452   for(j= 1;j<=m+1;j++)b[j]= 0.0;
  453   b[ld+1]= 1.0;
  454   lubksb(a,m+1,indx,b);
  455   for(kk= 1;kk<=np;kk++)c[kk]= 0.0;
  456   for(k= -nl;k<=nr;k++){
  457   sum= b[1];
  458   fac= 1.0;
  459   for(mm= 1;mm<=m;mm++)sum+= b[mm+1]*(fac*= k);
  460   kk= ((np-k)%np)+1;
  461   c[kk]= sum;
  462   }
  463   free_dvector(b,1,m+1);
  464   free_dmatrix(a,1,m+1,1,m+1);
  465   free_ivector(indx,1,m+1);
  466   return(0);
  467   }
  468   
  469   /*:36*/
  470   #line 1224 "./sgfilter.w"
  471   
  472   /*37:*/
  473   #line 1775 "./sgfilter.w"
  474   
  475   char sgfilter(double yr[],double yf[],int mm,int nl,int nr,int ld,int m){
  476   int np= nl+1+nr;
  477   double*c;
  478   char retval;
  479   
  480   #if CONVOLVE_WITH_NR_CONVLV 
  481   c= dvector(1,mm);
  482   retval= sgcoeff(c,np,nl,nr,ld,m);
  483   if(retval==0)
  484   convlv(yr,mm,c,np,1,yf);
  485   free_dvector(c,1,mm);
  486   #else 
  487   int j;
  488   long int k;
  489   c= dvector(1,nl+nr+1);
  490   retval= sgcoeff(c,np,nl,nr,ld,m);
  491   if(retval==0){
  492   for(k= 1;k<=nl;k++){
  493   for(yf[k]= 0.0,j= -nl;j<=nr;j++){
  494   if(k+j>=1){
  495   yf[k]+= c[(j>=0?j+1:nr+nl+2+j)]*yr[k+j];
  496   }
  497   }
  498   }
  499   for(k= nl+1;k<=mm-nr;k++){
  500   for(yf[k]= 0.0,j= -nl;j<=nr;j++){
  501   yf[k]+= c[(j>=0?j+1:nr+nl+2+j)]*yr[k+j];
  502   }
  503   }
  504   for(k= mm-nr+1;k<=mm;k++){
  505   for(yf[k]= 0.0,j= -nl;j<=nr;j++){
  506   if(k+j<=mm){
  507   yf[k]+= c[(j>=0?j+1:nr+nl+2+j)]*yr[k+j];
  508   }
  509   }
  510   }
  511   }
  512   free_dvector(c,1,nr+nl+1);
  513   #endif
  514   return(retval);
  515   }
  516   
  517   /*:37*/
  518   #line 1225 "./sgfilter.w"
  519   
  520   /*38:*/
  521   #line 1830 "./sgfilter.w"
  522   
  523   /*39:*/
  524   #line 1844 "./sgfilter.w"
  525   
  526   short pathcharacter(int ch){
  527   return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
  528   (ch=='-')||(ch=='+'));
  529   }
  530   
  531   /*:39*/
  532   #line 1831 "./sgfilter.w"
  533   
  534   char*strip_away_path(char filename[]){
  535   int j,k= 0;
  536   while(pathcharacter(filename[k]))k++;
  537   j= (--k);
  538   while(isalnum((int)(filename[j])))j--;
  539   j++;
  540   return(&filename[j]);
  541   }
  542   
  543   /*:38*/
  544   #line 1226 "./sgfilter.w"
  545   
  546   /*40:*/
  547   #line 1854 "./sgfilter.w"
  548   
  549   /*41:*/
  550   #line 1908 "./sgfilter.w"
  551   
  552   void hl(const char*format,...){
  553   va_list args;
  554   char line[1024];
  555   
  556   va_start(args,format);
  557   vsprintf(line,format,args);
  558   va_end(args);
  559   sprintf(line+strlen(line),"\n");
  560   fprintf(stdout,"%s",line);
  561   return;
  562   }
  563   
  564   /*:41*/
  565   #line 1855 "./sgfilter.w"
  566   
  567   void showsomehelp(void){
  568   hl("Usage: %s [options]",progname);
  569   hl("Options:");
  570   hl(" -h, --help");
  571   hl("    Display this help message and exit clean.");
  572   hl(" -i, --inputfile <str>");
  573   hl("    Specifies the file name from which unfiltered data is to be read.");
  574   hl("    The input file should describe the input as two columns contain-");
  575   hl("    ing $x$- and $y$-coordinates of the samples.");
  576   hl(" -o, --outputfile <str>");
  577   hl("    Specifies the file name to which filtered data is to be written,");
  578   hl("    again in a two-column format containing $x$- and $y$-coordinates");
  579   hl("    of the filtered samples. If this option is omitted, the generated");
  580   hl("    filtered data will instead be written to the console (terminal).");
  581   hl(" -nl <nl>");
  582   hl("    Specifies the number of samples nl to use to the 'left' of the");
  583   hl("    basis sample in the regression window (kernel). The total number");
  584   hl("    of samples in the window will be nL+nR+1.");
  585   hl(" -nr <nr>");
  586   hl("    Specifies the number of samples nr to use to the 'right' of the");
  587   hl("    basis sample in the regression window (kernel). The total number");
  588   hl("    of samples in the window will be nL+nR+1.");
  589   hl(" -m <m>");
  590   hl("    Specifies the order m of the polynomial to use in the regression");
  591   hl("    analysis leading to the Savitzky-Golay coefficients. Typical");
  592   hl("    values are between m=2 and m=6. Beware of too high values, which");
  593   hl("    easily makes the regression too sensitive, with an oscillatory");
  594   hl("    result.");
  595   hl(" -ld <ld>");
  596   hl("    Specifies the order of the derivative to extract from the ");
  597   hl("    Savitzky--Golay smoothing algorithm. For regular Savitzky-Golay");
  598   hl("    smoothing of the input data as such, use ld=0. For the Savitzky-");
  599   hl("    Golay smoothing and extraction of derivatives, set ld to the");
  600   hl("    order of the desired derivative and make sure that you correctly");
  601   hl("    interpret the scaling parameters as described in 'Numerical");
  602   hl("    Recipes in C', 2nd Edn (Cambridge University Press, New York,");
  603   hl("    1994).");
  604   hl(" -v, --verbose");
  605   hl("    Toggle verbose mode. (Default: Off.)  This option should always");
  606   hl("    be omitted whenever no  output file has been specified (that is");
  607   hl("    to say, omit any --verbose or -v option whenever --outputfile or");
  608   hl("    -o has been omitted), as the verbose logging otherwise will");
  609   hl("    contaminate the filtered data stream written to the console");
  610   hl("    (terminal).");
  611   }
  612   
  613   /*:40*/
  614   #line 1227 "./sgfilter.w"
  615   
  616   /*42:*/
  617   #line 1925 "./sgfilter.w"
  618   
  619   long int num_coordinate_pairs(FILE*file){
  620   double tmp;
  621   int tmpch;
  622   long int mm= 0;
  623   fseek(file,0L,SEEK_SET);
  624   while((tmpch= getc(file))!=EOF){
  625   ungetc(tmpch,file);
  626   fscanf(file,"%lf",&tmp);
  627   fscanf(file,"%lf",&tmp);
  628   mm++;
  629   tmpch= getc(file);
  630   while((tmpch!=EOF)&&(!isdigit(tmpch)))tmpch= getc(file);
  631   if(tmpch!=EOF)ungetc(tmpch,file);
  632   }
  633   fseek(file,0L,SEEK_SET);
  634   return(mm);
  635   }
  636   
  637   /*:42*/
  638   #line 1228 "./sgfilter.w"
  639   
  640   
  641   /*:22*/
  642   #line 978 "./sgfilter.w"
  643   
  644   
  645   int main(int argc,char*argv[])
  646   {
  647   /*16:*/
  648   #line 1046 "./sgfilter.w"
  649   
  650   int no_arg;
  651   int nl= DEFAULT_NL;
  652   int nr= DEFAULT_NR;
  653   int ld= DEFAULT_LD;
  654   int m= DEFAULT_M;
  655   long int k,mm= 0;
  656   double*x,*yr,*yf;
  657   char input_filename[NCHMAX]= "",output_filename[NCHMAX]= "";
  658   char verbose= 0;
  659   FILE*file;
  660   
  661   /*:16*/
  662   #line 982 "./sgfilter.w"
  663   
  664   /*17:*/
  665   #line 1064 "./sgfilter.w"
  666   
  667   progname= strip_away_path(argv[0]);
  668   no_arg= argc;
  669   while(--argc){
  670   if(!strcmp(argv[no_arg-argc],"-o")||
  671   !strcmp(argv[no_arg-argc],"--outputfile")){
  672   --argc;
  673   strcpy(output_filename,argv[no_arg-argc]);
  674   }else if(!strcmp(argv[no_arg-argc],"-i")||
  675   !strcmp(argv[no_arg-argc],"--inputfile")){
  676   --argc;
  677   strcpy(input_filename,argv[no_arg-argc]);
  678   }else if(!strcmp(argv[no_arg-argc],"-h")||
  679   !strcmp(argv[no_arg-argc],"--help")){
  680   showsomehelp();
  681   exit(0);
  682   }else if(!strcmp(argv[no_arg-argc],"-v")||
  683   !strcmp(argv[no_arg-argc],"--verbose")){
  684   verbose= (verbose?0:1);
  685   }else if(!strcmp(argv[no_arg-argc],"-nl")){
  686   --argc;
  687   if(!sscanf(argv[no_arg-argc],"%d",&nl)){
  688   log("Error in '-nl' option.");
  689   exit(1);
  690   }
  691   }else if(!strcmp(argv[no_arg-argc],"-nr")){
  692   --argc;
  693   if(!sscanf(argv[no_arg-argc],"%d",&nr)){
  694   log("Error in '-nr' option.");
  695   exit(1);
  696   }
  697   }else if(!strcmp(argv[no_arg-argc],"-ld")){
  698   --argc;
  699   if(!sscanf(argv[no_arg-argc],"%d",&ld)){
  700   log("Error in '-ld' option.");
  701   exit(1);
  702   }
  703   }else if(!strcmp(argv[no_arg-argc],"-m")){
  704   --argc;
  705   if(!sscanf(argv[no_arg-argc],"%d",&m)){
  706   log("Error in '-m' option.");
  707   exit(1);
  708   }
  709   }else{
  710   log("Unrecognized option '%s'.",argv[no_arg-argc]);
  711   showsomehelp();
  712   exit(1);
  713   }
  714   }
  715   if(verbose)fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
  716   
  717   /*:17*/
  718   #line 983 "./sgfilter.w"
  719   
  720   /*18:*/
  721   #line 1121 "./sgfilter.w"
  722   
  723   if(!strcmp(input_filename,"")){
  724   log("No input file specified! (Please use the '--inputfile' option.)");
  725   log("Execute '%s --help' for help.",progname);
  726   exit(1);
  727   }
  728   if((file= fopen(input_filename,"r"))==NULL){
  729   log("Could not open %s for loading raw data!",input_filename);
  730   exit(1);
  731   }
  732   mm= num_coordinate_pairs(file);
  733   if(mm<nl+nr+1){
  734   log("Error: The number M=%ld of data points must be at least nl+nr+1=%d",
  735   mm,nl+nr+1);
  736   log("Please check your -nl or -nr options.");
  737   exit(1);
  738   }
  739   if(verbose){
  740   log("Loading %ld unfiltered samples from %s...",mm,input_filename);
  741   log(" ... allocating memory for storage ...");
  742   }
  743   x= dvector(1,mm);
  744   yr= dvector(1,mm);
  745   #if CONVOLVE_WITH_NR_CONVLV
  746   yf= dvector(1,2*mm);
  747   #else
  748   yf= dvector(1,mm);
  749   #endif
  750   if(verbose)
  751   log(" ... scanning %s for input data ...",input_filename);
  752   for(k= 1;k<=mm;k++){
  753   fscanf(file,"%lf",&x[k]);
  754   fscanf(file,"%lf",&yr[k]);
  755   }
  756   fclose(file);
  757   if(verbose)
  758   log(" ... done. Input now residing in RAM.");
  759   
  760   /*:18*/
  761   #line 984 "./sgfilter.w"
  762   
  763   /*19:*/
  764   #line 1162 "./sgfilter.w"
  765   
  766   if(!sgfilter(yr,yf,mm,nl,nr,ld,m)){
  767   if(verbose)
  768   log("Successfully performed Savitzky-Golay filtering.");
  769   }else{
  770   if(verbose)
  771   log("Error: Could not perform Savitzky-Golay filtering.");
  772   }
  773   
  774   /*:19*/
  775   #line 985 "./sgfilter.w"
  776   
  777   /*20:*/
  778   #line 1178 "./sgfilter.w"
  779   
  780   if(!strcmp(output_filename,"")){
  781   if(verbose)
  782   log("Writing %ld filtered samples to console...",mm);
  783   }else{
  784   if(verbose)
  785   log("Writing %ld filtered samples to %s...",mm,output_filename);
  786   if((file= freopen(output_filename,"w",stdout))==NULL){
  787   log("Error: Unable to redirect stdout stream to file %s.",
  788   output_filename);
  789   exit(1);
  790   }
  791   }
  792   for(k= 1;k<=mm;k++)fprintf(stdout,"%1.8f %1.8f\n",x[k],yf[k]);
  793   #ifdef UNIX_LIKE_OS
  794   freopen("/dev/tty","a",stdout);
  795   #endif
  796   if(verbose)log(" ... done.");
  797   
  798   /*:20*/
  799   #line 986 "./sgfilter.w"
  800   
  801   /*21:*/
  802   #line 1199 "./sgfilter.w"
  803   
  804   free_dvector(x,1,mm);
  805   free_dvector(yr,1,mm);
  806   #if CONVOLVE_WITH_NR_CONVLV
  807   free_dvector(yf,1,2*mm);
  808   #else
  809   free_dvector(yf,1,mm);
  810   #endif
  811   
  812   /*:21*/
  813   #line 987 "./sgfilter.w"
  814   
  815   return(EXIT_SUCCESS);
  816   }
  817   
  818   /*:13*/
  819   

Return to previous page

Generated by ::viewsrc::

Last modified Saturday 17 Dec 2011