speana_cont.c

实现ipvlbi数据记录

/**************************************************************/
/*                                                            */
/*   Program Name :  speana_cont                              */
/*                                                            */
/*      Version:  0.00  2002-09-15                            */
/*                0.10  2002-09-15                            */
/*                                                            */
/*      Copyright (c) 2002 T.Kondo/CRL All Right Reserved     */
/*                                                            */
/**************************************************************/
/*
 俬俹亅倁俴俛俬儃乕僪庢摼僨乕僞偺僗儁僋僩儖夝愅
        乮僆僼儔僀儞張棟乯
         侾昩僨乕僞僼傽僀儖傪巊偭偰偺弨幚帪娫張棟
         僨乕僞偼 sampling2 偺 mode=6 傪巊偭偰廂廤偡傞偙偲

          T.Kondo 2002.5.20
               V0.10 2002.9.15

    憱傜偣曽 speana_cont filename [mode [sekibun [pmode [comment [soffset]]]]]
         偙偙偱 dirname -- 侾昩僨乕僞僼傽僀儖偺僼傽僀儖柤
                            捠忢 sampler.dat
                            0 偲偡傞偲僇儗儞僩僨傿儗僋僩儕偺丂sampler.dat 
               mode -- 廃攇悢幉儌乕僪
                       0: 儘僌僗働乕儖
                       1: 儕僯傾僗働乕儖乮僨僼僅儖僩乯 
                      -value 帺屓憡娭娭悢昞帵乮value揰奼戝乯
                      -1: 帺屓憡娭娭悢僼儖儗儞僕
                      
               sekibun  -- 愊暘帪娫(惍悢昩扨埵乯
                       僨僼僅儖僩偼 侾乮昩乯
               pmode -- 僾儘僢僩昞帵僨僶僀僗儌乕僪
                       0: XTERM 偍傛傃 PostScript僼傽僀儖弌椡乮僨僼僅儖僩乯
                       1: PostScrip 弌椡偺傒 
                       2: XTERM 弌椡偺傒 
               comment  -- 僐儊儞僩乮僌儔僼忋晹偵昞帵乯丅僗儁乕僗傪娷傑側偄偙偲丅
               soffset -- 僗僞乕僩帪崗僆僼僙僢僩乮昩乯
*/

#include <cpgplot.h>
#include <stdio.h>
#include <errno.h>
#include <fcntl.h>
#include <unistd.h>
#include <sys/stat.h>
//#include <sys/tdsio.h>
#include <math.h>
#ifdef linux
#include <sys/time.h>
#else
#include <time.h>
#endif
#include "libfft.h"
#include "libipvlbi.h"  // IP-VLBI梡娭悢孮(by T.Kondo)
//char dev[] = "/dev/tds0";

#ifndef PIRA
#define PIRA  1.745329252E-2    // PI/180  degree to radian
#endif
#ifndef PIDE
#define PIDE  57.29577951       // 180/PI radian to degree
#endif
#ifndef TWOPI
#define TWOPI 6.283185307       // 2PI
#endif
#ifndef PI
#define PI    3.141592653
#endif

/*    Offset from dB to dBm conversion from A/D multibit sampling data */ 
#define OFFSET1  16.0     //1 bit sampling full power convert to 0 dB (theoretical)
      /* following constants are provisional  (measured on 2002/5/18) */
#define OFFSET2   9.0     //2 bit sampling data offset measured at 200kHz for PC(vlbi11)
#define OFFSET4  -1.0     //4 bit sampling data offset measured at 200kHz for PC(vlbi11)
#define OFFSET8 -26.0     //8 bit sampling data offset measured at 200kHz for PC(vlbi11)

typedef struct ip_vlbi_file
{
	unsigned char bit32[4];
	unsigned char bit64[8];
	//char buf[65536];
	int sfreq;
	int icnt;   //  current position number
       char* fname; //僆乕僾儞偟偰偄傞僼傽僀儖柤
       float fsample;  //僒儞僾儕儞僌廃攇悢(Hz)
	int   adbit;  // AD價僢僩挿
	long int seconds; // 00h00m00s偐傜偺宱夁昩悢
       int hh, mm, ss;
	int spatn; // sync pattern
	int syncok; // =1 for 32bit all 1 detection
	int format; // =0 for old format , =1 for new format
	FILE *stream;
       fpos_t pos;
	int fh;
	int ch;
	int numch;
	int ierr;
       double dtime;
       unsigned int bytesread;
} ip_vlbi_file_t;

int checkheader(ip_vlbi_file_t *S)
{
       /* Read 64bit Data header check to get sampling parameters */
       int iok, i;

       S->bytesread = fread( S->bit64, sizeof( char ), 8, S->stream );
       //printf("bytes read 1 %d\n",bytesread); //debug

       if( ferror( S->stream ) || S->bytesread==0 )
       {
           printf("checkheader: File EOF or Read Error!\n");
	       i=fclose( S->stream ) ;
		   return -1;
	   }


       headerchkn(S->bit64,&iok,&S->spatn,&S->adbit,&S->sfreq,&S->numch,
                        &S->seconds,&S->hh,&S->mm,&S->ss);
       if (iok == 0)   // sync detection
	   {
          printf("New Format (header 64bits) Sync Detected!!\n");
		  printf("File : %s\n",S->fname);
          printf("  A/D(bits) = %d  # of ch = %d  Sampling Freq(kHz) = %d\n",
                    S->adbit,S->numch, S->sfreq);
          printf("  Time %02d:%02d:%02d   Sec in Day = %d\n",
                    S->hh,S->mm,S->ss,S->seconds);
          S->icnt=0;          
          return 0;
	   }
       else
       {
           printf("First 64bit is not header block!\n");
	       i=fclose( S->stream ) ;
           return -1;
       }
}

int fwdNsec(unsigned char* ibuf, int nsec, int numb, int imax,ip_vlbi_file_t *dat)
{
    int n, i;
    if(nsec <= 0) return 0;
	//printf("fwdNsec: nsec, numb, imax %d %d %d\n",nsec,numb,imax); //debug
    for(n=0; n<nsec; n++)
    {
        for(i=0; i< imax; i++)
        {
            dat->bytesread = fread( ibuf, sizeof( char ), numb, dat->stream );
            //printf("i, bytesread %d %d %d\n",i,dat->bytesread,ferror( dat->stream ));
            if( ferror( dat->stream ) || dat->bytesread==0 ){
		        return -1;
            }
		}
        if(checkheader(&(*dat)) < 0){
             return -1;
        }
     }
     return 0;
}

void plotspe_a(float *spow, int npos, int sfreq, int adbit, char *lab, char *lab2, int mode)
/*
        mode = 0; log frequency scale   1; linear frequency scale
*/
{
      float cofst,fmaxkhz, fminkhz, fsampl, df, ymax, ymin, pmax, pmin, maxdbm, mindbm;
      float *fdat, *powdat;
      float x,y;
      char htitl[60], xtit[60], ytit[60];
      int n, npmax;

      strcpy(htitl,"POWER SPECTRUM");

      /* frequency axis generation */

      npmax=npos/2;

      fdat   = (float *)malloc(npmax*sizeof(float));
      powdat = (float *)malloc(npmax*sizeof(float));

      fsampl=1000.0*(float)sfreq;
      df=fsampl/(float)npos;
      if(mode == 0) {    // log scale
         fmaxkhz=log10(fsampl/2000.0);
         fminkhz=log10(df/1000.0);
         for(n=0; n<npmax; n++){fdat[n]=log10((float)n*df/1000.0); }
      } else {          //  linear scale
         fmaxkhz=fsampl/2000.0;
         fminkhz=0.0;
         for(n=0; n<npmax; n++){fdat[n]=(float)n*df/1000.0; }
      }

      strcpy(xtit,"FREQUENCY (kHz)");
      strcpy(ytit,"POWER (dBm)");
      if( adbit==1 ){strcpy(ytit,"POWER (dB)"); }
      ymax=10.0;
      ymin=-80.0;
      if (adbit == 8 ){ ymin=-100; }    // extend dynamic range 

      /* data max, min check */
      for(n=0; n<npmax; n++)
      {
         powdat[n]=spow[n];
         if(n==0)
         {
            pmax=powdat[n];
            pmin=powdat[n];
         }
         else
         {
            if(powdat[n] > pmax){ pmax=powdat[n]; }
            if(powdat[n] < pmin){ pmin=powdat[n]; }
         }
      }

      maxdbm=pmax;
      mindbm=pmin;

      //goto exitend; //debug

      /* open graphics device */
      /*   move to main */
       //if (cpgopen("?") < 1)
      // if (cpgopen("/GF") < 1)
      // {
      //    goto exitend;
      // }

      cpgsch(2.0);  // set text height 2.0  (default is 1)
      cpgslw(2);    // set line width

      /* define coordinate range of graph */
      if(mode == 0){   
         cpgenv(fminkhz, fmaxkhz, ymin, ymax, 0, 10);   // log scale
      } else {         
         cpgenv(fminkhz, fmaxkhz, ymin, ymax, 0, 0);   // linear scale
      }
      
      /* label the axes */
       cpglab(xtit,ytit,htitl);

      /* sub label write*/
       cpgsch(1.5);
       cpgmtxt("T",1.5,0.5,0.5,lab);
       cpgmtxt("T",0.5,0.5,0.5,lab2);
       cpgsch(2.0);
      
      /* check data */
      if (pmax == pmin)
      {
          cpgmtxt("B",-6.0,0.5,0.5,"NO TIME VARIABLE");
          cpgmtxt("B",-4.0,0.5,0.5,"COMPONENT FOUND");
          /* label the axes */
          
      }
      else
      {


          /* plot the line graph */
          fdat[0]=fdat[1];
          powdat[0]=powdat[1];  // to avoid ugry plot
          //cpgslw(1);    // set line width
          cpgline(npmax,fdat,powdat); 

      }
      /* close the graphics device */
      /* move to main */
      //cpgclos();
  exitend: 
     free(powdat);
     free(fdat);
}

void plotautocor(float *spow, int npos, int sfreq, int adbit, char *lab, char *lab2, int amode)
/*
        mode = 0; log frequency scale   1; linear frequency scale
*/
{
      float cofst,tmax, tmin, fsampl, dt, ymax, ymin, pmax, pmin, maxdbm, mindbm;
      float *tdat, *pdat;
      float x,y;
      char htitl[60], xtit[60], ytit[60];
      int n, k, npmax, nhalf;

      strcpy(htitl,"AUTO CORRELATION (ABS)");

      /* time axis generation */

      if(amode <= 1) {
         npmax=npos;
      } else {
         npmax=amode;
      }
      tdat = (float *)malloc(npmax*sizeof(float));
      pdat = (float *)malloc(npmax*sizeof(float));
      nhalf=npmax/2;

      fsampl=1000.0*(float)sfreq;
      dt=1.0/fsampl;
      for(n=0; n<npmax; n++){
          tdat[n]=(float)(n-nhalf)*dt*1.0e6;   // usec unit
      }
      tmax=tdat[npmax-1];
      tmin=tdat[0];

      strcpy(xtit,"DELAY TIME (us)");
      strcpy(ytit,"CORRELATION AMP");
      ymax=1.1;
      ymin=0.0;

      /* data range limited copy and max, min check */
     
      for(n=0; n<npmax; n++)
      {
         k=n-nhalf;
         if(k < 0) k=npos+k;
         pdat[n]=spow[k];
         if(pdat[n]<0) pdat[n]=-pdat[n];  // to take absolute value
         if(n==0)
         {
            pmax=pdat[n];
            pmin=pdat[n];
         }
         else
         {
            if(pdat[n] > pmax){ pmax=pdat[n]; }
            if(pdat[n] < pmin){ pmin=pdat[n]; }
         }
      }
      

      cpgsch(2.0);  // set text height 2.0  (default is 1)
      cpgslw(2);    // set line width

      /* define coordinate range of graph */
      cpgenv(tmin, tmax, ymin, ymax, 0, 0);   // linear scale
      
      /* label the axes */
       cpglab(xtit,ytit,htitl);

      /* sub label write*/
       cpgsch(1.5);
       cpgmtxt("T",1.5,0.5,0.5,lab);
       cpgmtxt("T",0.5,0.5,0.5,lab2);
       cpgsch(2.0);
      /* check data */
      if (pmax == pmin)
      {
          cpgmtxt("B",-6.0,0.5,0.5,"NO DATA FOR PLOT");
          
      }
      else
      {
      
        /* plot the line graph */
         cpgline(npmax,tdat,pdat); 
      }

      /* move to main */
  exitend: 
     free(pdat);
     free(tdat);
}


int main(int argc, char *argv[])
{
       int spatn, sfreq, adbit, numch, kpos;
       int sfreqt,adbitt,numcht;
       int amode;    // >0 auto correlation plot mode 
       long int seconds;
       long int startsec,startsecdata;
       long int startoffset; // start time offset in sec
       int hh, mm, ss;
       int nloop,sekibun;
       int i,iok, imax, m, k, kaisu;