Contents of file 'sgfilter/sgfilter.c':
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
Generated by ::viewsrc::