Merge pull request #2 from ymaan4/write_plot_data

README.md

@@ -6,8 +6,8 @@ periodic RFI from GMRT time-domain data. Over the time, rfiClean has evolved
 to mitigate any spiky (in time or frequency) RFI as well, and from any SIGPROC
 filterbank format data file. It is written primarily in C. For handling the
 filterbank format I/O, rfiClean uses several modules from Duncan Lorimer's SIGPROC
-(https://sourceforge.net/p/sigproc/wiki/Home/; thank you Dunc) most of which are
-included here in the src/ext/ folder.
+(https://sourceforge.net/p/sigproc/wiki/Home/; thank you Dunc) which are included
+here in the src/ext/ folder (some of these codes are also modified suitably).
 
 #### Dependencies:
 FFTW3
@@ -17,8 +17,12 @@ PGPLOT
 
 
 ### Bash-script based parallel processing
-For faster processing, use the bash script in bin/ to use rfiClean on different parts of a single data file simultaneously and then combine the output products at the end --- parallel processing in a rather crude way.
+For faster processing, use the bash script in bin/ to use rfiClean on different parts of a single data file simultaneously and then combine the output products at the end --- parallel processing in a rather crude but efficient way.
 
+### Diagnostic plot
+rfiClean produces a diagnostic plot showing which Fourier frequencies are
+mitigated from the data and how frequently. The plot data are also output
+in a file.
 
 
 Yogesh Maan  <ymaan4[@]gmail.com>

include/header.h

@@ -15,6 +15,6 @@ double analog_power[2];
 
 /* added frequency table for use with non-contiguous data */
 //double frequency_table[4096]; /* note limited number of channels */
-double frequency_table[8192]; /* note limited number of channels */
+double frequency_table[32768]; /* note limited number of channels */
 //double frequency_table[16384]; /* note limited number of channels */
 long int npuls; /* added for binary pulse profile format */

include/rficlean.h

@@ -27,7 +27,7 @@ void bcast_header();
 void rficlean_data(FILE *input, FILE *output);
 float vmax(float *vec, int n);
 float vmin(float *vec, int n);
-void plot_data(char outdev[]);
+void plot_data(char outdev[], int wpout);
 
 /* include these from sigproc-4.3 */
 void slaCldj ( int iy, int im, int id, double *djm, int *j );

src/ext/read_header.c

@@ -14,7 +14,13 @@ void get_string(FILE *inputfile, int *nbytes, char string[]) /* includefile */
   strcpy(string,"ERROR");
   fread(&nchar, sizeof(int), 1, inputfile);
   if (feof(inputfile)) exit(0);
-  if (nchar>80 || nchar<1) return;
+  if (nchar>80 || nchar<1){
+    printf("\nnchar: %d ; more than allowed limit of 80.\n",nchar);
+    *nbytes=sizeof(int);
+    fread(string, 79, 1, inputfile);
+    printf("Partial string: '%s'\n",string);
+    return;
+  }
   *nbytes=sizeof(int);
   fread(string, nchar, 1, inputfile);
   string[nchar]='\0';

src/plot_data.c

@@ -23,12 +23,17 @@ float vmax(float *vec, int n)
   return (maxvalue);
 }
 
+struct headrecord
+{
+   long int bl,totalsamp;
+   int nchn,numints;
+   double sampt,fchan1,chanwidth,startmjd;
+};
 
-
-void plot_data(char outdev[])
+void plot_data(char outdev[], int wpout)
 {
-  char plotdev[100];
-  int i,j,k,npt, i1,i2,j1,j2;
+  char plotdev[100], plotout[100]={ '\0' };
+  int i,j,k,n, npt, i1,i2,j1,j2;
   long int ntot;
   float *xarr, *yarr, *fcstat;
   double *tmparr;
@@ -38,6 +43,7 @@ void plot_data(char outdev[])
   float tt, ft, th, fh;   /* thin and fat thicknesses and heights */
   float lm, rm, tm, bm;   /* LRTB margins */
   int ii, mincol, maxcol, numcol;
+  struct headrecord head;
 
   npt = maxhist*nsub;
   ntot = nchans*naddt;
@@ -48,8 +54,41 @@ void plot_data(char outdev[])
   //sprintf(plotdev, "%s/CPS", plotdev);
   strcpy(plotdev, outdev);
 
+  /*-----------------------------------------------------------*/
+  /* write the plot data out, if asked for */
+  if (wpout>0){
+      n = strlen(outdev);
+      //printf("\n %d %c %c %s\n",n,outdev[n-2],outdev[n-1],outdev);
+      if(outdev[n-2]=='p' && outdev[n-1]=='s'){
+        strncpy(plotout, outdev, n-2);
+        strcat(plotout, "pdat");
+      } else {
+        strcpy(plotout,"abc.pdat");
+      }
+      printf ("\nPlot data will be written in: %s\n", plotout);
+
+      FILE *optr;
+      optr = fopen(plotout,"wb");
+      //fprintf(optr,"%ld  %d  %ld  %f ",naddt,nchans,ntot,tsamp);
+      //fprintf(optr,"%f  %f   %f ",fch1,foff,(fch1+nchans*foff));
+      head.bl = naddt;
+      head.totalsamp = totsamp;
+      head.nchn = nchans;
+      head.numints = nints;
+      head.sampt = tsamp;
+      head.fchan1 = fch1;
+      head.chanwidth = foff;      
+      head.startmjd = tstart;
+      fwrite(&head,sizeof(head),1,optr);
+      fwrite(chanstat,sizeof(float),nchans,optr);
+      fwrite(fftstat,sizeof(float),ntot,optr);
+      fclose(optr);
+  }
+  /*-----------------------------------------------------------*/
+
 
   /* Open and prep the device */
+  strcat(plotdev,"/CPS");
   cpgopen(plotdev);
   cpgpap(10.25, 8.5 / 11.0);
   cpgpage();
@@ -237,7 +276,7 @@ void plot_data(char outdev[])
         cpgswin(xl, xh, yl, yh);
         cpgbox("BCNLST", 0.0, 0, "BNST", 0.0, 0);
         //cpgbox("", 0.0, 0, "CST", 0.0, 0);
-        cpgmtxt("B", 2.6, 0.5, 0.5, "Fluc. Frequency (Hz)");
+        cpgmtxt("B", 2.6, 0.5, 0.5, "Fourier Frequency (Hz)");
         cpgmtxt("L", 2.1, 0.5, 0.5, "Radio Frequency (MHz)");
         /* Label */
         left = lm + 3.0 * ft + 4.0 * tt;
@@ -324,7 +363,7 @@ void plot_data(char outdev[])
         cpgbin(naddt/2, xarr, yarr, 1);
         //for (i=0;i<naddt/2;i++) printf("%f  %f\n",xarr[i],yarr[i]);
         cpgbox("CMLST", 0.0, 0, "", 0.0, 0);
-        cpgmtxt("T", 1.8, 0.5, 0.5, "Fluc. Frequency (Hz)");
+        cpgmtxt("T", 1.8, 0.5, 0.5, "Fourier Frequency (Hz)");
      /*=======================================================*/
 
      /*==== Chan-reject stats ==========================*/

src/rficlean.c

@@ -105,11 +105,12 @@ void print_version(char *program, char *argument)
 
 void main (int argc, char *argv[])
 {
-  int i, nc, headersize, headerless=0,gm=0;
+  int i, ibits, nc, headersize, headerless=0,gm=0;
   char string[80];
 
   /* set up default global variables */
   obits=headerless=naddt=nsamp=0;
+  ibits=0;
   fthresh = 6.0;
   forcefthresh = 1000.0;
   rthresh = 4.0;
@@ -159,6 +160,9 @@ void main (int argc, char *argv[])
 	/* get and open file for output */
         strcpy(outfile,argv[++i]);
 	output=fopen(outfile,"wb");
+      } else if (strings_equal(argv[i],"-ibits")) {
+        i++;
+        ibits=atoi(argv[i]);
       } else if (strings_equal(argv[i],"-T")) {
 	i++;
 	nsamp=atoi(argv[i]); 
@@ -246,6 +250,7 @@ void main (int argc, char *argv[])
     printf ("\n Reading the GMRT header & time-stamp file...\n");
     read_gmheader(gminfofile, gmhdrfile);
     if (naddt <= 1) naddt=4096;
+    if (ibits > 0) nbits=ibits;
     if (obits == 0) obits=nbits;
     if( nsamp>0) tsamp = tsamp*nsamp ; 
     if (!headerless) {
@@ -270,6 +275,8 @@ void main (int argc, char *argv[])
     case 2:
       nchans=1;
       break;
+    case 6:
+      break;
     default:
       error_message("ERROR: Input data to rficlean is not in filterbank format");
       break;
@@ -292,5 +299,6 @@ void main (int argc, char *argv[])
   /* finally clean and output the data */
   rficlean_data(input,output);
 
+  printf("\n");
   exit(0);
 }

src/rficlean_data.c

@@ -19,7 +19,7 @@ void rficlean_data(FILE *input, FILE *output)
   float *fblock,min,max,hpower,realtime,fsaved[2], *ts0dm, atemp;
   unsigned short *sblock;
   unsigned char  *cblock;
-  int nsaved=0,opened=0;
+  int nsaved=0,opened=0, wpout=10;
   long int ns,nsblk,nout,iter,i,j,k, jt1, jt2, iblock;
   long int itemp, isum, nsize, istart, ii, jj, kk, n0;
 
@@ -34,7 +34,7 @@ void rficlean_data(FILE *input, FILE *output)
   nsub = 1;
   nsubchans = nchans/nsub ;
   strcpy(plotdevice, psfile);
-  strcat(plotdevice,"/CPS");
+  ////strcat(plotdevice,"/CPS");
   nsize = naddt;
   if(nsize < nchans) nsize = nchans;
   nsize = 2*nsize;
@@ -244,9 +244,12 @@ void rficlean_data(FILE *input, FILE *output)
   printf (" Data rfiClean-ed and written out! \n\n");
   printf (" Now making diagnostic plots ... ");
 
-  plot_data(plotdevice);
+  plot_data(plotdevice,wpout);
   printf (" Done! \n ");
 
+  fclose(input);
+  fclose(output);
+
   fftw_destroy_plan ( fplan );
   fftw_destroy_plan ( bplan );
   fftw_free ( in );