unpackdata.c

C编写的格式转换程序

    ms_log (1, "%s: forward/reverse integration constants:  X0: %d  XN: %d\n",
	    UNPACK_SRCNAME, *px0, *pxn);
  
  /* Decode compressed data in each frame */
  for (fn = 0; fn < num_data_frames; fn++)
    {
      
      ctrl = pf->ctrl;
      if ( swapflag ) ms_gswap4a (&ctrl);
      
      for (wn = 0; wn < VALS_PER_FRAME; wn++)
	{
	  if (nd >= num_samples) break;
	  
	  compflag = (ctrl >> ((VALS_PER_FRAME-wn-1)*2)) & 0x3;
	  
	  switch (compflag)
	    {
	    case STEIM2_SPECIAL_MASK:
	      /* Headers info -- skip it */
	      break;
	      
	    case STEIM2_BYTE_MASK:
	      /* Next 4 bytes are 4 1-byte differences */
	      for (i=0; i < 4 && nd < num_samples; i++, nd++)
		*diff++ = pf->w[wn].byte[i];
	      break;
	      
	    case STEIM2_123_MASK:
	      val = pf->w[wn].fw;
	      if ( swapflag ) ms_gswap4a (&val);
	      dnib =  val >> 30 & 0x3;
	      switch (dnib)
		{
		case 1:	/* 1 30-bit difference */
		  bits = 30; n = 1; m1 = 0x3fffffff; m2 = 0x20000000; break;
		case 2:	/* 2 15-bit differences */
		  bits = 15; n = 2; m1 = 0x00007fff; m2 = 0x00004000; break;
		case 3:	/* 3 10-bit differences */
		  bits = 10; n = 3; m1 = 0x000003ff; m2 = 0x00000200; break;
		default:	/*  should NEVER get here  */
		  ms_log (2, "msr_unpack_steim2(%s): invalid compflag, dnib, fn, wn = %d, %d, %d, %d\n", 
			  UNPACK_SRCNAME, compflag, dnib, fn, wn);
		  return MS_STBADCOMPFLAG;
		}
	      /*  Uncompress the differences */
	      for (i=(n-1)*bits; i >= 0 && nd < num_samples; i-=bits, nd++)
		{
		  *diff = (val >> i) & m1;
		  *diff = (*diff & m2) ? *diff | ~m1 : *diff;
		  diff++;
		}
	      break;
	      
	    case STEIM2_567_MASK:
	      val = pf->w[wn].fw;
	      if ( swapflag ) ms_gswap4a (&val);
	      dnib =  val >> 30 & 0x3;
	      switch (dnib)
		{
		case 0:	/*  5 6-bit differences  */
		  bits = 6; n = 5; m1 = 0x0000003f; m2 = 0x00000020; break;
		case 1:	/*  6 5-bit differences  */
		  bits = 5; n = 6; m1 = 0x0000001f; m2 = 0x00000010; break;
		case 2:	/*  7 4-bit differences  */
		  bits = 4; n = 7; m1 = 0x0000000f; m2 = 0x00000008; break;
		default:
		  ms_log (2, "msr_unpack_steim2(%s): invalid compflag, dnib, fn, wn = %d, %d, %d, %d\n", 
			  UNPACK_SRCNAME, compflag, dnib, fn, wn);
		  return MS_STBADCOMPFLAG;
		}
	      /* Uncompress the differences */
	      for (i=(n-1)*bits; i >= 0 && nd < num_samples; i-=bits, nd++)
		{
		  *diff = (val >> i) & m1;
		  *diff = (*diff & m2) ? *diff | ~m1 : *diff;
		  diff++;
		}
	      break;
	      
	    default:
	      /* Should NEVER get here */
	      ms_log (2, "msr_unpack_steim2(%s): invalid compflag, fn, wn = %d, %d, %d - nsamp: %d\n",
		      UNPACK_SRCNAME, compflag, fn, wn, nd);
	      return MS_STBADCOMPFLAG;
	    }
	}
      ++pf;
    }
    
  /* Test if the number of samples implied by the data frames is the
   * same number indicated in the header.
   */
  if ( nd != num_samples )
    {
      ms_log (2, "msr_unpack_steim2(%s): number of samples indicated in header (%d) does not equal data (%d)\n",
	      UNPACK_SRCNAME, num_samples, nd);
    }

  /*	For now, assume sample count in header to be correct.		*/
  /*	One way of "trimming" data from a block is simply to reduce	*/
  /*	the sample count.  It is not clear from the documentation	*/
  /*	whether this is a valid or not, but it appears to be done	*/
  /*	by other program, so we should not complain about its effect.	*/
  
  nr = req_samples;
  
  /* Compute first value based on last_value from previous buffer.	*/
  /* The two should correspond in all cases EXCEPT for the first	*/
  /* record for each component (because we don't have a valid xn from	*/
  /* a previous record).  Although the Steim compression algorithm	*/
  /* defines x(-1) as 0 for the first record, this only works for the	*/
  /* first record created since coldstart of the datalogger, NOT the	*/
  /* first record of an arbitrary starting record.	                */
  
  /* In all cases, assume x0 is correct, since we don't have x(-1).	*/
  data = databuff;
  diff = diffbuff;
  last_data = *px0;
  if (nr > 0)
    *data = *px0;

  /* Compute all but first values based on previous value               */
  prev = data - 1;
  while (--nr > 0 && --nd > 0)
    last_data = *++data = *++diff + *++prev;
  
  /* If a short count was requested compute the last sample in order    */
  /* to perform the integrity check comparison                          */
  while (--nd > 0)
    last_data = *++diff + last_data;
  
  /* Verify that the last value is identical to xn = rev. int. constant */
  if (last_data != *pxn)
    {
      ms_log (2, "%s: Data integrity check for Steim-2 failed, last_data=%d, xn=%d\n",
	      UNPACK_SRCNAME, last_data, *pxn);
    }
  
  return ((req_samples < num_samples) ? req_samples : num_samples);
}  /* End of msr_unpack_steim2() */


/* Defines for GEOSCOPE encoding */
#define GEOSCOPE_MANTISSA_MASK 0x0fff   /* mask for mantissa */
#define GEOSCOPE_GAIN3_MASK 0x7000      /* mask for gainrange factor */
#define GEOSCOPE_GAIN4_MASK 0xf000      /* mask for gainrange factor */
#define GEOSCOPE_SHIFT 12               /* # bits in mantissa */

/************************************************************************
 *  msr_unpack_geoscope:                                                *
 *                                                                      *
 *  Unpack GEOSCOPE gain ranged data (demultiplexed only) encoded       *
 *  miniSEED data and place in supplied buffer.                         *
 *                                                                      *
 *  Return: # of samples returned.                                      *
 ************************************************************************/
int msr_unpack_geoscope
 (const char   *edata,		/* ptr to encoded data.			*/
  int		num_samples,	/* number of data samples in total.     */
  int		req_samples,	/* number of data desired by caller.	*/
  float	       *databuff,	/* ptr to unpacked data array.		*/
  int           encoding,       /* specific GEOSCOPE encoding type      */
  int		swapflag)	/* if data should be swapped.	        */
{
  int nd = 0;		/* # of data points in packet.		*/
  int mantissa;		/* mantissa from SEED data */
  int gainrange;	/* gain range factor */
  int exponent;		/* total exponent */
  int k;
  uint64_t exp2val;
  int16_t sint;
  double dsample = 0.0;
  
  union {
    uint8_t b[4];
    uint32_t i;
  } sample32;
  
  if (num_samples < 0) return 0;
  if (req_samples < 0) return 0;

  /* Make sure we recognize this as a GEOSCOPE encoding format */
  if ( encoding != DE_GEOSCOPE24 &&
       encoding != DE_GEOSCOPE163 &&
       encoding != DE_GEOSCOPE164 )
    {
      ms_log (2, "msr_unpack_geoscope(%s): unrecognized GEOSCOPE encoding: %d\n",
	      UNPACK_SRCNAME, encoding);
      return -1;
    }
  
  for (nd=0; nd<req_samples && nd<num_samples; nd++)
    {
      switch (encoding)
	{
	case DE_GEOSCOPE24:
	  sample32.i = 0;
	  if ( swapflag )
	    for (k=0; k < 3; k++)
	      sample32.b[2-k] = edata[k];
	  else
	    for (k=0; k < 3; k++)
	      sample32.b[1+k] = edata[k];
	  
	  mantissa = sample32.i;

	  /* Take 2's complement for mantissa for overflow */
	  if (mantissa > MAX24) 
	    mantissa -= 2 * (MAX24 + 1);
	  
	  /* Store */
	  dsample = (double) mantissa;
	  
	  break;
	case DE_GEOSCOPE163:
	  memcpy (&sint, edata, sizeof(int16_t));
	  if ( swapflag ) ms_gswap2a(&sint);
	  
	  /* Recover mantissa and gain range factor */
	  mantissa = (sint & GEOSCOPE_MANTISSA_MASK);
	  gainrange = (sint & GEOSCOPE_GAIN3_MASK) >> GEOSCOPE_SHIFT;
	  
	  /* Exponent is just gainrange for GEOSCOPE */
	  exponent = gainrange;
	  
	  /* Calculate sample as mantissa / 2^exponent */
	  exp2val = (uint64_t) 1 << exponent;
	  dsample = ((double) (mantissa-2048)) / exp2val;
	  
	  break;
	case DE_GEOSCOPE164:
	  memcpy (&sint, edata, sizeof(int16_t));
	  if ( swapflag ) ms_gswap2a(&sint);
	  
	  /* Recover mantissa and gain range factor */
	  mantissa = (sint & GEOSCOPE_MANTISSA_MASK);
	  gainrange = (sint & GEOSCOPE_GAIN4_MASK) >> GEOSCOPE_SHIFT;
	  
	  /* Exponent is just gainrange for GEOSCOPE */
	  exponent = gainrange;
	  
	  /* Calculate sample as mantissa / 2^exponent */
	  exp2val = (uint64_t) 1 << exponent;
	  dsample = ((double) (mantissa-2048)) / exp2val;
	  
	  break;
	}
      
      /* Save sample in output array */
      databuff[nd] = (float) dsample;
      
      /* Increment edata pointer depending on size */
      switch (encoding)
	{
	case DE_GEOSCOPE24:
	  edata += 3;
	  break;
	case DE_GEOSCOPE163:
	case DE_GEOSCOPE164:
	  edata += 2;
	  break;
	}
    }
  
  return nd;
}  /* End of msr_unpack_geoscope() */


/* Defines for SRO encoding */
#define SRO_MANTISSA_MASK 0x0fff   /* mask for mantissa */
#define SRO_GAINRANGE_MASK 0xf000  /* mask for gainrange factor */
#define SRO_SHIFT 12               /* # bits in mantissa */

/************************************************************************
 *  msr_unpack_sro:                                                     *
 *                                                                      *
 *  Unpack SRO gain ranged data encoded miniSEED data and place in      *
 *  supplied buffer.                                                    *
 *                                                                      *
 *  Notes from original rdseed routine:                                 *
 *  SRO data are represented according to the formula                   *
 *                                                                      *
 *  sample = M * (b exp {[m * (G + agr)] + ar})                         *
 *                                                                      *
 *  where                                                               *
 *	sample = seismic data sample                                    *
 *	M      = mantissa                                               *
 *	G      = gain range factor                                      *
 *	b      = base to be exponentiated = 2 for SRO                   *
 *	m      = multiplier  = -1 for SRO                               *
 *	agr    = term to be added to gain range factor = 0 for SRO      *
 *	ar     = term to be added to [m * (gr + agr)]  = 10 for SRO     *
 *	exp    = exponentiation operation                               *
 *	Data are stored in two bytes as follows:                        *
 *		fedc ba98 7654 3210 = bit number, power of two          *
 *		GGGG MMMM MMMM MMMM = form of SEED data                 *
 *		where G = gain range factor and M = mantissa            *
 *	Masks to recover gain range and mantissa:                       *
 *		fedc ba98 7654 3210 = bit number = power of two         *
 *		0000 1111 1111 1111 = 0x0fff     = mask for mantissa    *
 *		1111 0000 0000 0000 = 0xf000     = mask for gain range  *
 *                                                                      *
 *  Return: # of samples returned.                                      *
 ************************************************************************/
int msr_unpack_sro
 (int16_t      *edata,		/* ptr to encoded data.			*/
  int		num_samples,	/* number of data samples in total.     */
  int		req_samples,	/* number of data desired by caller.	*/
  int32_t      *databuff,	/* ptr to unpacked data array.		*/
  int		swapflag)	/* if data should be swapped.	        */
{
  int32_t nd = 0;	/* sample count */
  int32_t mantissa;	/* mantissa */
  int32_t gainrange;	/* gain range factor */
  int32_t add2gr;       /* added to gainrage factor */
  int32_t mult;         /* multiplier for gain range */
  int32_t add2result;   /* added to multiplied gain rage */
  int32_t exponent;	/* total exponent */
  uint16_t sint;
  int32_t sample;
  
  if (num_samples < 0) return 0;
  if (req_samples < 0) return 0;
  
  add2gr = 0;
  mult = -1;
  add2result = 10;
  
  for (nd=0; nd<req_samples && nd<num_samples; nd++)
    {
      memcpy (&sint, &edata[nd], sizeof(int16_t));
      if ( swapflag ) ms_gswap2a(&sint);
      
      /* Recover mantissa and gain range factor */
      mantissa = (sint & SRO_MANTISSA_MASK);
      gainrange = (sint & SRO_GAINRANGE_MASK) >> SRO_SHIFT;
      
      /* Take 2's complement for mantissa */
      if ( mantissa > MAX12 )
	mantissa -= 2 * (MAX12 + 1);
      
      /* Calculate exponent, SRO exponent = 0..10 */
      exponent = (mult * (gainrange + add2gr)) + add2result;
      
      if ( exponent < 0 || exponent > 10 )
	{
	  ms_log (2, "msr_unpack_sro(%s): SRO gain ranging exponent out of range: %d\n",
		  UNPACK_SRCNAME, exponent);
	  return MS_GENERROR;
	}
      
      /* Calculate sample as mantissa * 2^exponent */
      sample = mantissa * ( (uint64_t) 1 << exponent );
      
      /* Save sample in output array */
      databuff[nd] = sample;
    }
  
  return nd;
}  /* End of msr_unpack_sro() */


/************************************************************************
 *  msr_unpack_dwwssn:                                                  *
 *                                                                      *
 *  Unpack DWWSSN encoded miniSEED data and place in supplied buffer.   *
 *                                                                      *
 *  Return: # of samples returned.                                      *
 ************************************************************************/
int msr_unpack_dwwssn
 (int16_t      *edata,		/* ptr to encoded data.			*/
  int		num_samples,	/* number of data samples in total.     */
  int		req_samples,	/* number of data desired by caller.	*/
  int32_t      *databuff,	/* ptr to unpacked data array.		*/
  int		swapflag)	/* if data should be swapped.	        */
{
  int32_t nd = 0;	/* sample count */
  int32_t sample;
  uint16_t sint;
  
  if (num_samples < 0) return 0;
  if (req_samples < 0) return 0;
  
  for (nd=0; nd<req_samples && nd<num_samples; nd++)
    {
      memcpy (&sint, &edata[nd], sizeof(uint16_t));
      if ( swapflag ) ms_gswap2a(&sint);
      sample = (int32_t) sint;
      
      /* Take 2's complement for sample */
      if ( sample > MAX16 )
	sample -= 2 * (MAX16 + 1);
      
      /* Save sample in output array */
      databuff[nd] = sample;
    }
  
  return nd;
}  /* End of msr_unpack_dwwssn() */