speana_cont.c

实现ipvlbi数据记录

       int smode,usampl,nfft,kplot;
	unsigned char bit64[8];
       unsigned char *ibuf;
	int nbyte=5000;
       FILE *stream;
       fpos_t pos;
       int bytesread;
       int bits, bytes, numb;
       //int mode=0;   // log frequency scale
       //int mode=1;   // linear frequency scale
       int mode, pmode; 
       int scount[4];
       float fsampl,tsampl,fact;
       int *idata;
       double *cdata;
       float *xdata, *ydata;
       double sum,dc;
       float pows, *spows;
       float pmax,pmin,pmaxdbm,pmindbm,cofst;
       char fname[60];
       char timelab[30];
       char lab[60], lab2[60], iplab[40], ltime1st[60];
       char cbuf[60];
       char labfsampl[10];
	//char outfile[] = "./tds.data";
	char* dflt_fname="sampler.dat";
       ip_vlbi_file_t X;

	if(argc == 1){
          printf("Usage: speana filename [mode [sekibun [pmode [comment [soffset]]]]]\n");
          printf("\n");
          printf("         where  filename ---- Data file name\n");
          printf("                              if 0 is given, tds.data will be used\n");
          printf("                mode -------- Frequenxy axis mode\n");
          printf("                              =0 : log scale\n");
          printf("                              =1 : linear scale(default)\n");
          printf("                              =-1 : Auto correlation function plot\n");
          printf("                                     instead of spectrum\n");
          printf("                              =-N : N points auto correlation function plot\n");
          printf("                sekibun --- Averaging period (integer seconds)\n");
          printf("                             default is 1\n");
          printf("                pmode ---- Plot device selection\n");
          printf("                           =0 :  XTERM and PostScript file (pgplot.ps) out\n");
          printf("                           =1 :  PostScript file (pgplot.ps) out only\n");
          printf("                           =2 :  XTERM only\n");
          printf("                          default is 0\n");
          printf("                comment --- Comment appeared on a graph\n");
          printf("                soffset --- Start time offset in sec (integer seconds)\n");
          printf("                             default is 0\n");
		 exit(0);
	}

       strcpy(fname,*(argv+1));
	if (fname[0]=='0'){
          strcpy(fname,dflt_fname);
	}
       X.fname=fname;
	if (argc>=3)
	{
	   mode = atoi(*(argv + 2 ));
          if(mode <0 ){
            amode=-mode;
            mode=0;
          } else {
            amode=0;
          }
	}
       else
       {
          mode=1;     // default is log frequency scale mode
          amode=0;
       }

	if (argc>=4)
	{
	   sekibun = atoi(*(argv + 3 ));
	}
       else
       {
          sekibun=1;     // default integration time (sec)
       }

	if (argc>=5)
	{
	   pmode = atoi(*(argv + 4 ));
	}
       else
       {
          pmode=0;     // default plot mode
       }


        /* get comment */

	if (argc>=6){
          strcpy(iplab,*(argv+5));
       } else {
       	strcpy(iplab,"  ");
       }
	if (argc>=7)
	{
	   startoffset = (long int)atoi(*(argv + 6 ));
	}
       else
       {
          startoffset=0;     // default startoffset
       }

	
	
loopentry:

       printf("Data File is %s\n",fname);
       if( (X.stream  = fopen( X.fname, "rb" )) == NULL )
       {
         //perror( "擖椡僼傽僀儖偑奐偗傑偣傫" );
         //exit( 1 );
	  // return -1;
         printf("File Open Error! \n");
		goto exitend2;
        }

       /* Data header check to get sampling parameters */
       if(checkheader(&X) < 0){
			i=fclose( X.stream ) ;
            goto exitend1;
       }

nextstep:

       /* parameter set */
        sfreq=X.sfreq;
        adbit=X.adbit;
        numch=X.numch;
        fsampl=1000.0*(float)sfreq;   // sampling frequency (Hz)
        tsampl=1.0/fsampl;            // sampling period (sec)
        bits=((long)fsampl)*adbit*numch;  // # of bits in 1 sec data except header
        bytes=bits/8;                     // # of bytes in 1 sec data except header
        X.fsample=fsampl;
        //sekibun=1;                       // Averaging period for Spectrum (sec)

        smode=1;                         // low resolution mode

        if(smode ==1){ usampl=10000; }
        else {
           usampl=50000;
           if(sfreq == 40){ usampl=40000; }
           if(sfreq == 100){ usampl=100000; }
           if(sfreq == 200){ usampl=200000; }
        }
        imax=(sfreq*1000)/usampl;
        numb=usampl*adbit*numch/8;   // # of bytes in unit # of samples

        printf("imax numb  %d %d\n",imax,numb);    //debug

        nfft=1024;
        if(usampl < 524288){nfft= 524288;}
        if(usampl < 262144){nfft= 262144;}
        if(usampl < 131072){nfft= 131072;}
        if(usampl < 65536){nfft= 65536;}
        if(usampl < 32768){nfft= 32768;}
        if(usampl < 16384){nfft= 16384;}

        ibuf = (unsigned char *)malloc(numb);
        idata = (int *)malloc(4*usampl*sizeof(int));
        cdata = (double *)malloc(2*nfft*sizeof(double));
        spows = (float *)malloc(4*nfft*sizeof(float));
        xdata = (float *)malloc(nfft*sizeof(float));
        ydata = (float *)malloc(nfft*sizeof(float));
        //autocor=(float *)malloc(4*nfft*sizeof(float));

        fact=((float)usampl)*((float)usampl);
        
        /* sampling frequency in text for print */
        if(sfreq == 40){strcpy(labfsampl,"40kHz"); }
        if(sfreq == 100){strcpy(labfsampl,"100kHz"); }
        if(sfreq == 200){strcpy(labfsampl,"200kHz"); }
        if(sfreq == 500){strcpy(labfsampl,"500kHz"); }
        if(sfreq == 1000){strcpy(labfsampl,"1MHz"); }
        if(sfreq == 2000){strcpy(labfsampl,"2MHz"); }
        if(sfreq == 4000){strcpy(labfsampl,"4MHz"); }
        if(sfreq == 8000){strcpy(labfsampl,"8MHz"); }
        if(sfreq == 18000){strcpy(labfsampl,"16MHz"); }


       /* foward startffset seconds */
       if(fwdNsec(ibuf, startoffset, numb, imax, &X) <0){
              goto exitend;
       }

        for(k=0; k<4*nfft; k++){spows[k]=0.0;}
        for(k=0; k<4; k++){scount[k]=0;}

        printf("  Time %02d:%02d:%02d   Sec in Day = %d\n",X.hh,X.mm,X.ss,X.seconds);
        sprintf(ltime1st,"Time(UT) %02d:%02d:%02d\0",X.hh,X.mm,X.ss);
        strcpy(lab2,ltime1st);
        strcat(lab2,"  ");
        strcat(lab2,iplab);
        

       nloop=0;
loop1:
        if(nloop >0) {
           if(checkheader(&X) < 0){ goto exitend; }
            printf("  Time %02d:%02d:%02d   Sec in Day = %d\n",X.hh,X.mm,X.ss,X.seconds);
        }

        for(i=0; i< imax; i++)
        {
             bytesread = fread( ibuf, sizeof( char ), numb, X.stream );
             //printf("bytes read 2 %d\n",bytesread); //debug
             parseADdata_2(ibuf,numb,adbit,idata, &kpos);
             //printf("numb adbit kpos %d %d %d\n",numb,adbit,kpos); //debug
             //printf("ibuf %d %d %d %d\n",ibuf[0],ibuf[1],ibuf[2],ibuf[3]); //debug
             //printf("idata %d %d %d %d\n",idata[0],idata[1],idata[2],idata[3]); //debug
             //for(k=0; k<kpos; k++){printf("k idata %d %d\n",k,idata[k]);} // debug
             for(m=0; m<numch; m++)
             {   
                kaisu=0;
                sum=0.0;
                for(k=0; k<2*nfft; k++){cdata[k]=0.0;} // initialize
                for(k=0; k<kpos/numch; k++)
                {
                   //if(m==0){printf("k idata %d %d\n",k,idata[k*numch+m]);} //debug
                   cdata[2*k]=(double)idata[k*numch+m];        
                   // cdata[2*k+1]=0.0;  
                   sum=sum+cdata[2*k];
                   kaisu++;
                }
             
                dc=sum/(double)kaisu;    // DC component
     
                //printf("m DC = %d %f\n",m,dc);  //debug
                
                for(k=0; k<kpos/numch; k++)
                {
                  //if(m==0){printf("cdata %f",cdata[2*k]);} //debug
                  cdata[2*k]-=dc;
                  //if(m==0){printf(" %f",cdata[2*k]);} //debug
                }
                //hanning(cdata,nfft);  // hanning filter
                //cdft(2*nfft,-1,cdata);
                if(amode ==0) {    // spectrum mode
                    hanning(cdata,usampl);  // hanning filter
                    cdft(2*nfft,1,cdata);
                } else {
                    cdft(2*nfft,-1,cdata);
                }
                for(k=0; k<nfft; k++)
                {
                   pows=(cdata[2*k]*cdata[2*k]+cdata[2*k+1]*cdata[2*k+1])/(double)fact;
                   spows[k*4+m]=spows[k*4+m]+pows;
                }

                scount[m]++;
             }
       }
       nloop++;
       if(nloop < sekibun){ goto loop1;}
       fgetpos(X.stream,&pos);  //debug
	i=fclose( X.stream ) ;
       printf("%d bytes read\n",pos); 

       if(amode > 0) {   // auto correlation function calc
          for(m=0; m<numch; m++){
             for(k=0; k<nfft; k++){
                spows[k*4+m]/=scount[m]; // nomalization
                cdata[2*k]=spows[k*4+m];
                cdata[2*k+1]=0.0;
             }
             if(m==1 || m==3) cdata[0]=1.0;  // to avoid math error
             cdft(2*nfft,1,cdata);
             for(k=0; k<nfft; k++){
               // if(k<10) printf("%d %e\n",k,cdata[2*k]);  // debug
               // if(k>nfft-10) printf("%d %e\n",k,cdata[2*k]);  // debug
                if(k==0) pmax=cdata[0];
                if(cdata[2*k] > pmax) pmax=cdata[2*k];
                spows[k*4+m]=cdata[2*k];
             }
             for(k=0; k<nfft; k++){  // normalize
                spows[k*4+m]/=pmax;
               // if(k<10) printf("   %d %e\n",k,spows[k*4+m]);  // debug
               // if(k>nfft-10) printf("   %d %e\n",k,spows[k*4+m]);  // debug
             }
             
          }
         
       } else {



       	/* convert power to dB  then dBm */       
       	cofst=OFFSET1;       //1 bit sampling full power convert to 0 dB
       	if( adbit==8){cofst=OFFSET8; }    // 8bit sampling offset
       	if( adbit==4){cofst=OFFSET4; }    // 4bit sampling offset
       	if( adbit==2){cofst=OFFSET2; }    // 4bit sampling offset


       	for(m=0; m<numch; m++){
          		for(k=0; k<nfft; k++){
              		spows[k*4+m]/=scount[m]; // nomalization
              		spows[k*4+m]=10.0*log10(spows[k*4+m]) + cofst; 
             		if(k==0){
               		 pmaxdbm=spows[k*4+m];
               		 pmindbm=spows[k*4+m];
              	} else {
                		if(spows[k*4+m] > pmaxdbm){ pmaxdbm=spows[k*4+m]; }
                		if(spows[k*4+m] < pmindbm){ pmindbm=spows[k*4+m]; }
              	}
           	}
           	pmax=pmaxdbm-cofst;
           	pmin=pmindbm-cofst;
           	printf("CH# %d: Maximum data (dB, dBm) is %f %f\n", m+1,pmax,pmaxdbm);      
        	}

        }

        /* plot spectrum (or auto cor) POSTSCRIPT OUT and XTERM*/

        for(kplot=0; kplot<2; kplot++)
        {
            if(pmode==1 && kplot==1) goto bypass;
            if(pmode==2 && kplot==0) goto bypass;
            //if (cpgopen("?") < 1)
            if(kplot==0) {
               if (cpgopen("/PS") < 1) goto exitend;
            }
            if(kplot==1) {
               if (cpgopen("/XTERM") < 1) goto exitend;
            }

            cpgsubp(2,2);  // divide plot pannel into 2X2

            for(m=0; m<numch; m++)
            {
                sprintf(lab,"CH#:%1d %1dbit %s sampling, Tavg(sec)=%3d\0",m+1,
                       adbit,labfsampl,sekibun);
                //printf("%s\n",lab);
                for(k=0; k<nfft; k++) 
                {
                   xdata[k]=spows[k*4+m];
                }
                if(amode > 0) {
                     plotautocor(xdata, nfft, sfreq, adbit, lab, lab2, amode);
                } else {
                     plotspe_a(xdata, nfft, sfreq, adbit, lab, lab2, mode);
                }

            }
            /* panel undivided to date write */
            cpgsubp(1,1);
            /* date etc write */
            cpgiden();
            /* close the graphics device */
            cpgclos();
            bypass:
        }

 exitend: ;
       
       free(ibuf);
       free(cdata);
       free(idata);
       free(spows);
       free(xdata);
       free(ydata);
	
 exitend1:;
	
	  goto loopentry;

 exitend2: ;

       return 0;
}