/************************************************************************
* Routines for decoding INT16, INT32, FLOAT32, FLOAT64, STEIM1,
* STEIM2, GEOSCOPE (24bit and gain ranged), CDSN, SRO and DWWSSN
* encoded data.
*
* modified: 2017.283
************************************************************************/
#include <memory.h>
#include <stdio.h>
#include <stdlib.h>
#include "libmseed.h"
#include "unpackdata.h"
/* Control for printing debugging information */
int decodedebug = 0;
/* Extract bit range. Byte order agnostic & defined when used with unsigned values */
#define EXTRACTBITRANGE(VALUE, STARTBIT, LENGTH) ((VALUE >> STARTBIT) & ((1U << LENGTH) - 1))
#define MAX12 0x7FFul /* maximum 12 bit positive # */
#define MAX14 0x1FFFul /* maximum 14 bit positive # */
#define MAX16 0x7FFFul /* maximum 16 bit positive # */
#define MAX24 0x7FFFFFul /* maximum 24 bit positive # */
/************************************************************************
* msr_decode_int16:
*
* Decode 16-bit integer data and place in supplied buffer as 32-bit
* integers.
*
* Return number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_int16 (int16_t *input, int samplecount, int32_t *output,
int outputlength, int swapflag)
{
int16_t sample;
int idx;
if (samplecount <= 0)
return 0;
if (!input || !output || outputlength <= 0)
return -1;
for (idx = 0; idx < samplecount && outputlength >= (int)sizeof (int32_t); idx++)
{
sample = input[idx];
if (swapflag)
ms_gswap2a (&sample);
output[idx] = (int32_t)sample;
outputlength -= sizeof (int32_t);
}
return idx;
} /* End of msr_decode_int16() */
/************************************************************************
* msr_decode_int32:
*
* Decode 32-bit integer data and place in supplied buffer as 32-bit
* integers.
*
* Return number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_int32 (int32_t *input, int samplecount, int32_t *output,
int outputlength, int swapflag)
{
int32_t sample;
int idx;
if (samplecount <= 0)
return 0;
if (!input || !output || outputlength <= 0)
return -1;
for (idx = 0; idx < samplecount && outputlength >= (int)sizeof (int32_t); idx++)
{
sample = input[idx];
if (swapflag)
ms_gswap4a (&sample);
output[idx] = sample;
outputlength -= sizeof (int32_t);
}
return idx;
} /* End of msr_decode_int32() */
/************************************************************************
* msr_decode_float32:
*
* Decode 32-bit float data and place in supplied buffer as 32-bit
* floats.
*
* Return number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_float32 (float *input, int samplecount, float *output,
int outputlength, int swapflag)
{
float sample;
int idx;
if (samplecount <= 0)
return 0;
if (!input || !output || outputlength <= 0)
return -1;
for (idx = 0; idx < samplecount && outputlength >= (int)sizeof (float); idx++)
{
memcpy (&sample, &input[idx], sizeof (float));
if (swapflag)
ms_gswap4a (&sample);
output[idx] = sample;
outputlength -= sizeof (float);
}
return idx;
} /* End of msr_decode_float32() */
/************************************************************************
* msr_decode_float64:
*
* Decode 64-bit float data and place in supplied buffer as 64-bit
* floats, aka doubles.
*
* Return number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_float64 (double *input, int samplecount, double *output,
int outputlength, int swapflag)
{
double sample;
int idx;
if (samplecount <= 0)
return 0;
if (!input || !output || outputlength <= 0)
return -1;
for (idx = 0; idx < samplecount && outputlength >= (int)sizeof (double); idx++)
{
memcpy (&sample, &input[idx], sizeof (double));
if (swapflag)
ms_gswap8a (&sample);
output[idx] = sample;
outputlength -= sizeof (double);
}
return idx;
} /* End of msr_decode_float64() */
/************************************************************************
* msr_decode_steim1:
*
* Decode Steim1 encoded miniSEED data and place in supplied buffer
* as 32-bit integers.
*
* Return number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_steim1 (int32_t *input, int inputlength, int samplecount,
int32_t *output, int outputlength, char *srcname,
int swapflag)
{
int32_t *outputptr = output; /* Pointer to next output sample location */
uint32_t frame[16]; /* Frame, 16 x 32-bit quantities = 64 bytes */
int32_t X0 = 0; /* Forward integration constant, aka first sample */
int32_t Xn = 0; /* Reverse integration constant, aka last sample */
int maxframes = inputlength / 64;
int frameidx;
int startnibble;
int nibble;
int widx;
int diffcount;
int idx;
union dword {
int8_t d8[4];
int16_t d16[2];
int32_t d32;
} * word;
if (inputlength <= 0)
return 0;
if (!input || !output || outputlength <= 0 || maxframes <= 0)
return -1;
if (decodedebug)
ms_log (1, "Decoding %d Steim1 frames, swapflag: %d, srcname: %s\n",
maxframes, swapflag, (srcname) ? srcname : "");
for (frameidx = 0; frameidx < maxframes && samplecount > 0; frameidx++)
{
/* Copy frame, each is 16x32-bit quantities = 64 bytes */
memcpy (frame, input + (16 * frameidx), 64);
/* Save forward integration constant (X0) and reverse integration constant (Xn)
and set the starting nibble index depending on frame. */
if (frameidx == 0)
{
if (swapflag)
{
ms_gswap4a (&frame[1]);
ms_gswap4a (&frame[2]);
}
X0 = frame[1];
Xn = frame[2];
startnibble = 3; /* First frame: skip nibbles, X0, and Xn */
if (decodedebug)
ms_log (1, "Frame %d: X0=%d Xn=%d\n", frameidx, X0, Xn);
}
else
{
startnibble = 1; /* Subsequent frames: skip nibbles */
if (decodedebug)
ms_log (1, "Frame %d\n", frameidx);
}
/* Swap 32-bit word containing the nibbles */
if (swapflag)
ms_gswap4a (&frame[0]);
/* Decode each 32-bit word according to nibble */
for (widx = startnibble; widx < 16 && samplecount > 0; widx++)
{
/* W0: the first 32-bit contains 16 x 2-bit nibbles for each word */
nibble = EXTRACTBITRANGE (frame[0], (30 - (2 * widx)), 2);
word = (union dword *)&frame[widx];
diffcount = 0;
switch (nibble)
{
case 0: /* 00: Special flag, no differences */
if (decodedebug)
ms_log (1, " W%02d: 00=special\n", widx);
break;
case 1: /* 01: Four 1-byte differences */
diffcount = 4;
if (decodedebug)
ms_log (1, " W%02d: 01=4x8b %d %d %d %d\n",
widx, word->d8[0], word->d8[1], word->d8[2], word->d8[3]);
break;
case 2: /* 10: Two 2-byte differences */
diffcount = 2;
if (swapflag)
{
ms_gswap2a (&word->d16[0]);
ms_gswap2a (&word->d16[1]);
}
if (decodedebug)
ms_log (1, " W%02d: 10=2x16b %d %d\n", widx, word->d16[0], word->d16[1]);
break;
case 3: /* 11: One 4-byte difference */
diffcount = 1;
if (swapflag)
ms_gswap4a (&word->d32);
if (decodedebug)
ms_log (1, " W%02d: 11=1x32b %d\n", widx, word->d32);
break;
} /* Done with decoding 32-bit word based on nibble */
/* Apply accumulated differences to calculate output samples */
if (diffcount > 0)
{
for (idx = 0; idx < diffcount && samplecount > 0; idx++, outputptr++)
{
if (outputptr == output) /* Ignore first difference, instead store X0 */
*outputptr = X0;
else if (diffcount == 4) /* Otherwise store difference from previous sample */
*outputptr = *(outputptr - 1) + word->d8[idx];
else if (diffcount == 2)
*outputptr = *(outputptr - 1) + word->d16[idx];
else if (diffcount == 1)
*outputptr = *(outputptr - 1) + word->d32;
samplecount--;
}
}
} /* Done looping over nibbles and 32-bit words */
} /* Done looping over frames */
/* Check data integrity by comparing last sample to Xn (reverse integration constant) */
if (outputptr != output && *(outputptr - 1) != Xn)
{
ms_log (1, "%s: Warning: Data integrity check for Steim1 failed, Last sample=%d, Xn=%d\n",
srcname, *(outputptr - 1), Xn);
}
return (outputptr - output);
} /* End of msr_decode_steim1() */
/************************************************************************
* msr_decode_steim2:
*
* Decode Steim2 encoded miniSEED data and place in supplied buffer
* as 32-bit integers.
*
* Return number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_steim2 (int32_t *input, int inputlength, int samplecount,
int32_t *output, int outputlength, char *srcname,
int swapflag)
{
int32_t *outputptr = output; /* Pointer to next output sample location */
uint32_t frame[16]; /* Frame, 16 x 32-bit quantities = 64 bytes */
int32_t X0 = 0; /* Forward integration constant, aka first sample */
int32_t Xn = 0; /* Reverse integration constant, aka last sample */
int32_t diff[7];
int32_t semask;
int maxframes = inputlength / 64;
int frameidx;
int startnibble;
int nibble;
int widx;
int diffcount;
int dnib;
int idx;
union dword {
int8_t d8[4];
int32_t d32;
} * word;
if (inputlength <= 0)
return 0;
if (!input || !output || outputlength <= 0 || maxframes <= 0)
return -1;
if (decodedebug)
ms_log (1, "Decoding %d Steim2 frames, swapflag: %d, srcname: %s\n",
maxframes, swapflag, (srcname) ? srcname : "");
for (frameidx = 0; frameidx < maxframes && samplecount > 0; frameidx++)
{
/* Copy frame, each is 16x32-bit quantities = 64 bytes */
memcpy (frame, input + (16 * frameidx), 64);
/* Save forward integration constant (X0) and reverse integration constant (Xn)
and set the starting nibble index depending on frame. */
if (frameidx == 0)
{
if (swapflag)
{
ms_gswap4a (&frame[1]);
ms_gswap4a (&frame[2]);
}
X0 = frame[1];
Xn = frame[2];
startnibble = 3; /* First frame: skip nibbles, X0, and Xn */
if (decodedebug)
ms_log (1, "Frame %d: X0=%d Xn=%d\n", frameidx, X0, Xn);
}
else
{
startnibble = 1; /* Subsequent frames: skip nibbles */
if (decodedebug)
ms_log (1, "Frame %d\n", frameidx);
}
/* Swap 32-bit word containing the nibbles */
if (swapflag)
ms_gswap4a (&frame[0]);
/* Decode each 32-bit word according to nibble */
for (widx = startnibble; widx < 16 && samplecount > 0; widx++)
{
/* W0: the first 32-bit quantity contains 16 x 2-bit nibbles */
nibble = EXTRACTBITRANGE (frame[0], (30 - (2 * widx)), 2);
diffcount = 0;
switch (nibble)
{
case 0: /* nibble=00: Special flag, no differences */
if (decodedebug)
ms_log (1, " W%02d: 00=special\n", widx);
break;
case 1: /* nibble=01: Four 1-byte differences */
diffcount = 4;
word = (union dword *)&frame[widx];
for (idx = 0; idx < diffcount; idx++)
{
diff[idx] = word->d8[idx];
}
if (decodedebug)
ms_log (1, " W%02d: 01=4x8b %d %d %d %d\n", widx, diff[0], diff[1], diff[2], diff[3]);
break;
case 2: /* nibble=10: Must consult dnib, the high order two bits */
if (swapflag)
ms_gswap4a (&frame[widx]);
dnib = EXTRACTBITRANGE (frame[widx], 30, 2);
switch (dnib)
{
case 0: /* nibble=10, dnib=00: Error, undefined value */
ms_log (2, "%s: Impossible Steim2 dnib=00 for nibble=10\n", srcname);
return -1;
break;
case 1: /* nibble=10, dnib=01: One 30-bit difference */
diffcount = 1;
semask = 1ul << (30 - 1); /* Sign extension from bit 30 */
diff[0] = EXTRACTBITRANGE (frame[widx], 0, 30);
diff[0] = (diff[0] ^ semask) - semask;
if (decodedebug)
ms_log (1, " W%02d: 10,01=1x30b %d\n", widx, diff[0]);
break;
case 2: /* nibble=10, dnib=10: Two 15-bit differences */
diffcount = 2;
semask = 1ul << (15 - 1); /* Sign extension from bit 15 */
for (idx = 0; idx < diffcount; idx++)
{
diff[idx] = EXTRACTBITRANGE (frame[widx], (15 - idx * 15), 15);
diff[idx] = (diff[idx] ^ semask) - semask;
}
if (decodedebug)
ms_log (1, " W%02d: 10,10=2x15b %d %d\n", widx, diff[0], diff[1]);
break;
case 3: /* nibble=10, dnib=11: Three 10-bit differences */
diffcount = 3;
semask = 1ul << (10 - 1); /* Sign extension from bit 10 */
for (idx = 0; idx < diffcount; idx++)
{
diff[idx] = EXTRACTBITRANGE (frame[widx], (20 - idx * 10), 10);
diff[idx] = (diff[idx] ^ semask) - semask;
}
if (decodedebug)
ms_log (1, " W%02d: 10,11=3x10b %d %d %d\n", widx, diff[0], diff[1], diff[2]);
break;
}
break;
case 3: /* nibble=11: Must consult dnib, the high order two bits */
if (swapflag)
ms_gswap4a (&frame[widx]);
dnib = EXTRACTBITRANGE (frame[widx], 30, 2);
switch (dnib)
{
case 0: /* nibble=11, dnib=00: Five 6-bit differences */
diffcount = 5;
semask = 1ul << (6 - 1); /* Sign extension from bit 6 */
for (idx = 0; idx < diffcount; idx++)
{
diff[idx] = EXTRACTBITRANGE (frame[widx], (24 - idx * 6), 6);
diff[idx] = (diff[idx] ^ semask) - semask;
}
if (decodedebug)
ms_log (1, " W%02d: 11,00=5x6b %d %d %d %d %d\n",
widx, diff[0], diff[1], diff[2], diff[3], diff[4]);
break;
case 1: /* nibble=11, dnib=01: Six 5-bit differences */
diffcount = 6;
semask = 1ul << (5 - 1); /* Sign extension from bit 5 */
for (idx = 0; idx < diffcount; idx++)
{
diff[idx] = EXTRACTBITRANGE (frame[widx], (25 - idx * 5), 5);
diff[idx] = (diff[idx] ^ semask) - semask;
}
if (decodedebug)
ms_log (1, " W%02d: 11,01=6x5b %d %d %d %d %d %d\n",
widx, diff[0], diff[1], diff[2], diff[3], diff[4], diff[5]);
break;
case 2: /* nibble=11, dnib=10: Seven 4-bit differences */
diffcount = 7;
semask = 1ul << (4 - 1); /* Sign extension from bit 4 */
for (idx = 0; idx < diffcount; idx++)
{
diff[idx] = EXTRACTBITRANGE (frame[widx], (24 - idx * 4), 4);
diff[idx] = (diff[idx] ^ semask) - semask;
}
if (decodedebug)
ms_log (1, " W%02d: 11,10=7x4b %d %d %d %d %d %d %d\n",
widx, diff[0], diff[1], diff[2], diff[3], diff[4], diff[5], diff[6]);
break;
case 3: /* nibble=11, dnib=11: Error, undefined value */
ms_log (2, "%s: Impossible Steim2 dnib=11 for nibble=11\n", srcname);
return -1;
break;
}
break;
} /* Done with decoding 32-bit word based on nibble */
/* Apply differences to calculate output samples */
if (diffcount > 0)
{
for (idx = 0; idx < diffcount && samplecount > 0; idx++, outputptr++)
{
if (outputptr == output) /* Ignore first difference, instead store X0 */
*outputptr = X0;
else /* Otherwise store difference from previous sample */
*outputptr = *(outputptr - 1) + diff[idx];
samplecount--;
}
}
} /* Done looping over nibbles and 32-bit words */
} /* Done looping over frames */
/* Check data integrity by comparing last sample to Xn (reverse integration constant) */
if (outputptr != output && *(outputptr - 1) != Xn)
{
ms_log (1, "%s: Warning: Data integrity check for Steim2 failed, Last sample=%d, Xn=%d\n",
srcname, *(outputptr - 1), Xn);
}
return (outputptr - output);
} /* End of msr_decode_steim2() */
/* Defines for GEOSCOPE encoding */
#define GEOSCOPE_MANTISSA_MASK 0x0FFFul /* mask for mantissa */
#define GEOSCOPE_GAIN3_MASK 0x7000ul /* mask for gainrange factor */
#define GEOSCOPE_GAIN4_MASK 0xf000ul /* mask for gainrange factor */
#define GEOSCOPE_SHIFT 12 /* # bits in mantissa */
/************************************************************************
* msr_decode_geoscope:
*
* Decode GEOSCOPE gain ranged data (demultiplexed only) encoded
* miniSEED data and place in supplied buffer as 32-bit floats.
*
* Return number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_geoscope (char *input, int samplecount, float *output,
int outputlength, int encoding,
char *srcname, int swapflag)
{
int idx = 0;
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 (!input || !output)
return -1;
if (samplecount <= 0 || outputlength <= 0)
return -1;
/* Make sure we recognize this as a GEOSCOPE encoding format */
if (encoding != DE_GEOSCOPE24 &&
encoding != DE_GEOSCOPE163 &&
encoding != DE_GEOSCOPE164)
{
ms_log (2, "msr_decode_geoscope(%s): unrecognized GEOSCOPE encoding: %d\n",
srcname, encoding);
return -1;
}
for (idx = 0; idx < samplecount && outputlength >= (int)sizeof (float); idx++)
{
switch (encoding)
{
case DE_GEOSCOPE24:
sample32.i = 0;
if (swapflag)
for (k = 0; k < 3; k++)
sample32.b[2 - k] = input[k];
else
for (k = 0; k < 3; k++)
sample32.b[1 + k] = input[k];
mantissa = sample32.i;
/* Take 2's complement for mantissa for overflow */
if ((unsigned long)mantissa > MAX24)
mantissa -= 2 * (MAX24 + 1);
/* Store */
dsample = (double)mantissa;
break;
case DE_GEOSCOPE163:
memcpy (&sint, input, 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, input, 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 */
output[idx] = (float)dsample;
outputlength -= sizeof (float);
/* Increment edata pointer depending on size */
switch (encoding)
{
case DE_GEOSCOPE24:
input += 3;
break;
case DE_GEOSCOPE163:
case DE_GEOSCOPE164:
input += 2;
break;
}
}
return idx;
} /* End of msr_decode_geoscope() */
/* Defines for CDSN encoding */
#define CDSN_MANTISSA_MASK 0x3FFFul /* mask for mantissa */
#define CDSN_GAINRANGE_MASK 0xC000ul /* mask for gainrange factor */
#define CDSN_SHIFT 14 /* # bits in mantissa */
/************************************************************************
* msr_decode_cdsn:
*
* Decode CDSN gain ranged data encoded miniSEED data and place in
* supplied buffer as 32-bit integers.
*
* Notes from original rdseed routine:
* CDSN data are compressed according to the formula
*
* sample = M * (2 exp G)
*
* where
* sample = seismic data sample
* M = mantissa; biased mantissa B is written to tape
* G = exponent of multiplier (i.e. gain range factor);
* key K is written to tape
* exp = exponentiation operation
* B = M + 8191, biased mantissa, written to tape
* K = key to multiplier exponent, written to tape
* K may have any of the values 0 - 3, as follows:
* 0 => G = 0, multiplier = 2 exp 0 = 1
* 1 => G = 2, multiplier = 2 exp 2 = 4
* 2 => G = 4, multiplier = 2 exp 4 = 16
* 3 => G = 7, multiplier = 2 exp 7 = 128
* Data are stored on tape in two bytes as follows:
* fedc ba98 7654 3210 = bit number, power of two
* KKBB BBBB BBBB BBBB = form of SEED data
* where K = key to multiplier exponent and B = biased mantissa
*
* Masks to recover key to multiplier exponent and biased mantissa
* from tape are:
* fedc ba98 7654 3210 = bit number = power of two
* 0011 1111 1111 1111 = 0x3fff = mask for biased mantissa
* 1100 0000 0000 0000 = 0xc000 = mask for gain range key
*
* Return number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_cdsn (int16_t *input, int samplecount, int32_t *output,
int outputlength, int swapflag)
{
int32_t idx = 0;
int32_t mantissa; /* mantissa */
int32_t gainrange; /* gain range factor */
int32_t mult = -1; /* multiplier for gain range */
uint16_t sint;
int32_t sample;
if (samplecount <= 0)
return 0;
for (idx = 0; idx < samplecount && outputlength >= (int)sizeof (int32_t); idx++)
{
memcpy (&sint, &input[idx], sizeof (int16_t));
if (swapflag)
ms_gswap2a (&sint);
/* Recover mantissa and gain range factor */
mantissa = (sint & CDSN_MANTISSA_MASK);
gainrange = (sint & CDSN_GAINRANGE_MASK) >> CDSN_SHIFT;
/* Determine multiplier from the gain range factor and format definition
* because shift operator is used later, these are powers of two */
if (gainrange == 0)
mult = 0;
else if (gainrange == 1)
mult = 2;
else if (gainrange == 2)
mult = 4;
else if (gainrange == 3)
mult = 7;
/* Unbias the mantissa */
mantissa -= MAX14;
/* Calculate sample from mantissa and multiplier using left shift
* mantissa << mult is equivalent to mantissa * (2 exp (mult)) */
sample = ((uint32_t)mantissa << mult);
/* Save sample in output array */
output[idx] = sample;
outputlength -= sizeof (int32_t);
}
return idx;
} /* End of msr_decode_cdsn() */
/* Defines for SRO encoding */
#define SRO_MANTISSA_MASK 0x0FFFul /* mask for mantissa */
#define SRO_GAINRANGE_MASK 0xF000ul /* mask for gainrange factor */
#define SRO_SHIFT 12 /* # bits in mantissa */
/************************************************************************
* msr_decode_sro:
*
* Decode SRO gain ranged data encoded miniSEED data and place in
* supplied buffer as 32-bit integers.
*
* 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 number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_sro (int16_t *input, int samplecount, int32_t *output,
int outputlength, char *srcname, int swapflag)
{
int32_t idx = 0;
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 (samplecount <= 0)
return 0;
add2gr = 0;
mult = -1;
add2result = 10;
for (idx = 0; idx < samplecount && outputlength >= (int)sizeof (int32_t); idx++)
{
memcpy (&sint, &input[idx], 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 ((unsigned long)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_decode_sro(%s): SRO gain ranging exponent out of range: %d\n",
srcname, exponent);
return MS_GENERROR;
}
/* Calculate sample as mantissa * 2^exponent */
sample = mantissa * ((uint64_t)1 << exponent);
/* Save sample in output array */
output[idx] = sample;
outputlength -= sizeof (int32_t);
}
return idx;
} /* End of msr_decode_sro() */
/************************************************************************
* msr_decode_dwwssn:
*
* Decode DWWSSN encoded miniSEED data and place in supplied buffer
* as 32-bit integers.
*
* Return number of samples in output buffer on success, -1 on error.
************************************************************************/
int
msr_decode_dwwssn (int16_t *input, int samplecount, int32_t *output,
int outputlength, int swapflag)
{
int32_t idx = 0;
int32_t sample;
uint16_t sint;
if (samplecount < 0)
return 0;
for (idx = 0; idx < samplecount && outputlength >= (int)sizeof (int32_t); idx++)
{
memcpy (&sint, &input[idx], sizeof (uint16_t));
if (swapflag)
ms_gswap2a (&sint);
sample = (int32_t)sint;
/* Take 2's complement for sample */
if ((unsigned long)sample > MAX16)
sample -= 2 * (MAX16 + 1);
/* Save sample in output array */
output[idx] = sample;
outputlength -= sizeof (int32_t);
}
return idx;
} /* End of msr_decode_dwwssn() */