speana2.c

实现ipvlbi数据记录

/**************************************************************/
/*                                                            */
/*   Program Name :  speana2                                  */
/*                                                            */
/*      Version:  1.00  2003-04-12                            */
/*                                                            */
/*                                                            */
/*      Copyright (c) 2002 T.Kondo/CRL All Right Reserved     */
/*                                                            */
/**************************************************************/
/*
 俬俹亅倁俴俛俬儃乕僪庢摼僨乕僞偺僗儁僋僩儖夝愅 偦偺俀
        乮僆僼儔僀儞張棟乯

          T.Kondo 2003.4.12

    憱傜偣曽 speana filename [smode [sekibun1 [skibun [soffset [comment [f1khz [f2khz]]]]]]]
         偙偙偱 filename -- 僨乕僞僼傽僀儖柤
                            0 偲偡傞偲僼傽僀儖偼
                           tds.data傪巊偆
               smode -- 廃攇悢暘夝擻儌乕僪
                       0: 係侽俋俇揰丂乮栺1倠俫倸@係俵俫倸僒儞僾儕儞僌乯乯
                       1: 4194304揰丂乮栺1俫倸@係俵俫倸僒儞僾儕儞僌乯
               sekibun1  -- 扨埵愊暘帪娫(惍悢昩扨埵乯
                       僨僼僅儖僩偼 侾乮昩乯
               sekibun  -- 慡扨埵愊暘帪娫(惍悢昩扨埵乯
                       僨僼僅儖僩偼 俁侽侽乮昩乯
               soffset -- 僗僞乕僩帪崗僆僼僙僢僩乮昩乯
                       僨僼僅儖僩偼 侽乮昩乯
               comment  -- 僐儊儞僩乮僌儔僼忋晹偵昞帵乯丅僗儁乕僗傪娷傑側偄偙偲丅
                           徣棯偟偨応崌偼夛榖儌乕僪擖椡偵側傞
                           拲丗僗儁乕僗傪娷傓僐儊儞僩偼夛榖儌乕僪偱擖椡偡傞偙偲丅
               f1khz  ---  僷儚乕嵟戝傪僒乕僠偡傞廃攇悢壓尷乮僨僼僅儖僩偼0)
               f2khz  ---  僷儚乕嵟戝傪僒乕僠偡傞廃攇悢忋尷乮僨僼僅儖僩偼價僨僆忋尷廃攇悢)
*/

#include <cpgplot.h>
#include <stdio.h>
#include <errno.h>
#include <fcntl.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.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, 
               double f1khz, double f2khz, float *pmaxdat, float *fdatmax, int nmax)
/*
        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;
      float fatmax;
      char htitl[60], xtit[60], ytit[60];
      int n, npmax,k, ncount,nmaxpos,i;
      //float pmaxdat[4],fdatmax[4];
	  //nmax=4;
      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;
         //fmaxkhz=990.0;
         //fminkhz=920.0;
      	 if(f1khz != 0.0) fminkhz=(float)f1khz;
      	 if(f2khz != 0.0) fmaxkhz=(float)f2khz;
         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;
      //ymax=-34.0;
      //ymin=-40.0;
      if (adbit == 8 ){ ymin=-100; }    // extend dynamic range 

      /* data max, min check */
	  //printf("max search start %d\n",npmax);
	  
	  for(n=0; n<npmax; n++) powdat[n]=spow[n];
	
      pmin= 1.0e30;
	  for(ncount=0; ncount<nmax; ncount++){
          pmax=-1.0e30;
		  k=0;
		  for(n=0; n<npmax; n++){
            if(fdat[n]>=fminkhz && fdat[n]<=fmaxkhz){
               if(powdat[n] > pmax){ 
                   pmax=powdat[n];
                   fatmax=fdat[n];
               	   pmaxdat[ncount]=pmax;
               	   fdatmax[ncount]=fatmax;
               	   nmaxpos=n;
               }
           	   if(ncount==0) {
           		    if(powdat[n] < pmin) pmin=powdat[n]; 
           	   }
               k++;
           	}
          }
		 powdat[nmaxpos]=-10.0e30;
	  	 //printf("pmax,fatmax, i %f %f %d\n",pmax,fatmax,i);
      }
	 /*
      printf("Number of points = %d\n",k);
      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]);
	 */
      maxdbm=pmaxdat[0];
      mindbm=pmin;
      ymax=pmaxdat[0];
      ymin=pmin;
	  //printf("ymax,ymin %f %f\n",ymax,ymin);
	
	  for(n=0; n<npmax; n++) powdat[n]=spow[n];
      //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,nloop1,sekibun,sekibun1,sekibun1st,npmax;
       int i,iok, imax, m, k, n, kaisu;
       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