speana2.c

实现ipvlbi数据记录

       //int mode=1;   // linear frequency scale
       int mode, pmode; 
       int scount[4];
       float fsampl,tsampl,fact;
	   double f1khz,f2khz;
       int *idata;
       double *cdata;
       float *xdata, *ydata, *fdat;
       double sum,dc;
       float pows, *spows;
       float pmax,pmin,pmaxdbm,pmindbm,cofst,df;
       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="tds.data";
       float pmaxdat[4],fdatmax[4];
	   int nmax=4;
       ip_vlbi_file_t X;

	   if(argc == 1){
          printf("Usage: speana2 filename [smode [sekibun1 [sekibun [soffset [comment [f1khz [f2khz]]]]]]]\n");
          printf("\n");
          printf("         where  filename ---- Data file name\n");
          printf("                              if 0 is given, tds.data will be used\n");
          printf("                smode -------- process mode\n");
          printf("                              =0 : 4096 points FFT (default)\n");
          printf("                                   (about 1 kHz resolution @4MHz sampling)\n");
          printf("                              =1 : 4194304 points FFT (about 1Hz resolution)\n");
          printf("                sekibun1 --- unit Averaging period (integer seconds)\n");
          printf("                             default is 1\n");
          printf("                sekibun ---- total averaging period\n");
          printf("                             default is 300\n");
          printf("                soffset --- Start time offset in sec (integer seconds)\n");
          printf("                             default is 0\n");
          printf("                comment --- Comment appeared on a graph\n");
          printf("                            if omitted, you will be asked later.\n");
          printf("                f1khz   ----  serach frequency range (lower boundary) (kHz)\n");
          printf("                              for power Maximum (default is 0)\n");
          printf("                f2khz   ----  serach frequency range (upper boundary) (kHz)\n");
          printf("                              for power Maximum (default is video frequency)\n");
		  exit(0);
	   }

       strcpy(fname,*(argv+1));
	   if (fname[0]=='0'){
          strcpy(fname,dflt_fname);
	   }
       X.fname=fname;
	
	   if (argc>=3){
	      smode = atoi(*(argv + 2 ));
	   } else {
		  smode=0;
       }

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

	   if (argc>=5){
	      sekibun = atoi(*(argv + 4 ));
	   } else {
          sekibun=300;     
       }

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

       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;
       }

        /* get comment */

	   if (argc>=7){
          strcpy(iplab,*(argv+6));
       } else {
          /* enter board ID */
          /*
          printf("Enter IP-VLBI BOARD ID -> ");
          read1line(iplab);      // 1 line read from stdin
          printf("Board ID is %s\n",iplab);      
          strcpy(lab2,"IP-VLBI Board ID : ");
          strcat(lab2,iplab);
          */

          /* enter any comment */

          printf("Enter any comment -> ");
          read1line(iplab);      // 1 line read from stdin
          //printf("Board ID is %s\n",iplab);      
          //strcat(lab2,iplab);
       }
	
	   if (argc>=8){
	      f1khz = atof(*(argv + 7 ));
	   } else {
	   	  f1khz=0.0;
	   }
	   if (argc>=9){
	   	  f2khz = atof(*(argv + 8 ));
	   } else {
	   	  f2khz=0.0;
	   }
	printf("f1khz,f2khz %f %f\n",f1khz,f2khz);
       /* Data header check to get sampling parameters */
       if(checkheader(&X) < 0){
            goto exitend;
       }
		
	
	   amode=0; mode=1;   // lenear scale and spectrum mode
	   pmode=0;

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

	    usampl=4000;
        //if(smode ==1){ usampl=10000; }
        if(smode ==1) usampl=4000000; 
    //    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;}
        nfft=4096;
	    if(smode==1) nfft=4194304;
	
	printf("usampl,nfft %d %d\n",usampl,nfft);
        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));
        npmax=nfft/2;

        fdat   = (float *)malloc(npmax*sizeof(float));
        //autocor=(float *)malloc(4*nfft*sizeof(float));
        df=fsampl/(float)nfft;
        for(n=0; n<npmax; n++) fdat[n]=(float)n*df/1000.0; 
      

        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;
       }


        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);
        
	    if( (stream  = fopen( "speana2.txt", "w" )) == NULL ){
            //perror( "擖椡僼傽僀儖偑奐偗傑偣傫" );
            //exit( 1 );
	        // return -1;
            printf("File Open Error! \n");
		    goto exitend2;
        }

	   //sekibun1=30;
	   //sekibun1=1;
	   nloop1=0;
	
	   fprintf(stream,"%d %d %d\n", nfft, sfreq, adbit);
	   fprintf(stream,"%d %d\n",sekibun1,sekibun);
	
       nloop=0;
	   sekibun1st=1;
loop_main:
        for(k=0; k<4*nfft; k++){spows[k]=0.0;}
        for(k=0; k<4; k++){scount[k]=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);
        }
    	if(sekibun1st==1){
        	sprintf(ltime1st,"Time(UT) %02d:%02d:%02d\0",X.hh,X.mm,X.ss);
            strcpy(lab2,ltime1st);
            strcat(lab2,"  ");
            strcat(lab2,iplab);
            sekibun1st=0;
	    }
        for(i=0; i< imax; i++){
             bytesread = fread( ibuf, sizeof( char ), numb, X.stream );
        	if(bytesread <= 0) goto exitend;
             //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=(float)((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++;
       nloop1++;
       if(nloop1 < sekibun1) goto loop1;
	   nloop1=0;

       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);
            	
            	if(kplot==0) fprintf(stream,"%s\n",lab2);
                for(k=0; k<nfft; k++) xdata[k]=spows[k*4+m];
            	
            	if(kplot==0){
            	   for(k=0; k<npmax; k++) {
                	fprintf(stream,"%d %f %f\n",k+1,fdat[k],xdata[k]);
                   }
            	}
                if(amode > 0) {
                     plotautocor(xdata, nfft, sfreq, adbit, lab, lab2, amode);
                } else {
                     plotspe_a(xdata, nfft, sfreq, adbit, lab, lab2, mode, f1khz, f2khz, pmaxdat, fdatmax, nmax);
                	if(kplot==0){
                        printf("1st Max %fdbm at %fkHz\n",pmaxdat[0],fdatmax[0]);
                        printf("2nd Max %fdbm at %fkHz\n",pmaxdat[1],fdatmax[1]);
                        printf("3rd Max %fdbm at %fkHz\n",pmaxdat[2],fdatmax[2]);
                        printf("4th Max %fdbm at %fkHz\n",pmaxdat[3],fdatmax[3]);
                	}
                }

            }
            /* panel undivided to date write */
            cpgsubp(1,1);
            /* date etc write */
            cpgiden();
            /* close the graphics device */
            cpgclos();
            bypass:
        }
	if(nloop < sekibun) {
		sekibun1st=1;
		goto loop_main;
	}
 exitend: ;
 	   i=fclose( X.stream ) ;
	   i=fclose(stream); 
       free(ibuf);
       free(cdata);
       free(idata);
       free(spows);
       free(xdata);
       free(ydata);
       free(fdat);

 exitend2: ;

       return 0;
}