tayttt
/
StartServerWPF
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
master
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '4f686ff5d0'
${ noResults }
StartServerWPF/mseedC/libmseed/unpackdata.c

952 lines
28 KiB
C
Raw Blame History

/************************************************************************
* Routines for decoding INT16, INT32, FLOAT32, FLOAT64, STEIM1,
* STEIM2, GEOSCOPE (24bit and gain ranged), CDSN, SRO and DWWSSN
* encoded data.
*
* This file is part of the miniSEED Library.
*
* Copyright (c) 2023 Chad Trabant, EarthScope Data Services
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
************************************************************************/
#include <memory.h>
#include <stdio.h>
#include <stdlib.h>
#include "libmseed.h"
#include "unpackdata.h"
/* 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, int64_t samplecount, int32_t *output,
int64_t 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_gswap2 (&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, int64_t samplecount, int32_t *output,
int64_t 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_gswap4 (&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, int64_t samplecount, float *output,
int64_t 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_gswap4 (&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, int64_t samplecount, double *output,
int64_t 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_gswap8 (&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, int64_t samplecount,
int32_t *output, int64_t outputlength, const char *srcname,
int swapflag)
{
uint32_t frame[16]; /* Frame, 16 x 32-bit quantities = 64 bytes */
int32_t diff[60]; /* Difference values for a frame, max is 15 x 4 (8-bit samples) */
int32_t Xn = 0; /* Reverse integration constant, aka last sample */
int64_t outputidx;
int maxframes = inputlength / 64;
int diffidx;
int frameidx;
int startnibble;
int nibble;
int widx;
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;
/* Make sure output buffer is sufficient for all output samples */
if (outputlength < (samplecount * sizeof (int32_t)))
{
ms_log (2, "%s(%s) Output buffer not large enough for decoded samples\n",
__func__, srcname);
return -1;
}
#if DECODE_DEBUG
ms_log (0, "Decoding %d Steim1 frames, swapflag: %d, srcname: %s\n",
maxframes, swapflag, (srcname) ? srcname : "");
#endif
for (frameidx = 0, outputidx = 0;
frameidx < maxframes && outputidx < samplecount;
frameidx++)
{
/* Copy frame, each is 16x32-bit quantities = 64 bytes */
memcpy (frame, input + (16 * frameidx), 64);
diffidx = 0;
/* 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_gswap4 (&frame[1]);
ms_gswap4 (&frame[2]);
}
output[0] = frame[1];
outputidx++;
Xn = frame[2];
startnibble = 3; /* First frame: skip nibbles, X0, and Xn */
#if DECODE_DEBUG
ms_log (0, "Frame %d: X0=%d Xn=%d\n", frameidx, output[0], Xn);
#endif
}
else
{
startnibble = 1; /* Subsequent frames: skip nibbles */
#if DECODE_DEBUG
ms_log (0, "Frame %d\n", frameidx);
#endif
}
/* Swap 32-bit word containing the nibbles */
if (swapflag)
ms_gswap4 (&frame[0]);
/* Decode each 32-bit word according to nibble */
for (widx = startnibble; widx < 16; 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];
switch (nibble)
{
case 0: /* 00: Special flag, no differences */
#if DECODE_DEBUG
ms_log (0, " W%02d: 00=special\n", widx);
#endif
break;
case 1: /* 01: Four 1-byte differences */
for (idx = 0; idx < 4; idx++)
{
diff[diffidx++] = word->d8[idx];
}
#if DECODE_DEBUG
ms_log (0, " W%02d: 01=4x8b %d %d %d %d\n", widx,
diff[diffidx - 4], diff[diffidx - 3], diff[diffidx - 2], diff[diffidx - 1]);
#endif
break;
case 2: /* 10: Two 2-byte differences */
for (idx = 0; idx < 2; idx++)
{
if (swapflag)
{
ms_gswap2 (&word->d16[idx]);
}
diff[diffidx++] = word->d16[idx];
}
#if DECODE_DEBUG
ms_log (0, " W%02d: 10=2x16b %d %d\n", widx,
diff[diffidx - 2], diff[diffidx - 1]);
#endif
break;
case 3: /* 11: One 4-byte difference */
if (swapflag)
{
ms_gswap4 (&word->d32);
}
diff[diffidx++] = word->d32;
#if DECODE_DEBUG
ms_log (0, " W%02d: 11=1x32b %d\n", widx, diff[diffidx - 1]);
#endif
break;
} /* Done with decoding 32-bit word based on nibble */
} /* Done looping over nibbles and 32-bit words */
/* Apply differences in this frame to calculate output samples,
* ignoring first difference for first frame */
for (idx = (frameidx == 0) ? 1 : 0;
idx < diffidx && outputidx < samplecount;
idx++, outputidx++)
{
output[outputidx] = output[outputidx - 1] + diff[idx];
}
} /* Done looping over frames */
/* Check data integrity by comparing last sample to Xn (reverse integration constant) */
if (outputidx == samplecount && output[outputidx - 1] != Xn)
{
ms_log (1, "%s: Warning: Data integrity check for Steim1 failed, Last sample=%d, Xn=%d\n",
srcname, output[outputidx - 1], Xn);
}
return (int)outputidx;
} /* 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, int64_t samplecount,
int32_t *output, int64_t outputlength, const char *srcname,
int swapflag)
{
uint32_t frame[16]; /* Frame, 16 x 32-bit quantities = 64 bytes */
int32_t diff[105]; /* Difference values for a frame, max is 15 x 7 (4-bit samples) */
int32_t Xn = 0; /* Reverse integration constant, aka last sample */
int64_t outputidx;
int maxframes = inputlength / 64;
int diffidx;
int frameidx;
int startnibble;
int nibble;
int widx;
int dnib;
int idx;
union dword {
int8_t d8[4];
int32_t d32;
} *word;
/* Bitfield specifications for sign extension of various bit-width values */
struct {signed int x:4;} s4;
struct {signed int x:5;} s5;
struct {signed int x:6;} s6;
struct {signed int x:10;} s10;
struct {signed int x:15;} s15;
struct {signed int x:30;} s30;
if (inputlength <= 0)
return 0;
if (!input || !output || outputlength <= 0 || maxframes <= 0)
return -1;
/* Make sure output buffer is sufficient for all output samples */
if (outputlength < (samplecount * sizeof (int32_t)))
{
ms_log (2, "%s(%s) Output buffer not large enough for decoded samples\n",
__func__, srcname);
return -1;
}
#if DECODE_DEBUG
ms_log (0, "Decoding %d Steim2 frames, swapflag: %d, srcname: %s\n",
maxframes, swapflag, (srcname) ? srcname : "");
#endif
for (frameidx = 0, outputidx = 0;
frameidx < maxframes && outputidx < samplecount;
frameidx++)
{
/* Copy frame, each is 16x32-bit quantities = 64 bytes */
memcpy (frame, input + (16 * frameidx), 64);
diffidx = 0;
/* 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_gswap4 (&frame[1]);
ms_gswap4 (&frame[2]);
}
output[0] = frame[1];
outputidx++;
Xn = frame[2];
startnibble = 3; /* First frame: skip nibbles, X0, and Xn */
#if DECODE_DEBUG
ms_log (0, "Frame %d: X0=%d Xn=%d\n", frameidx, output[0], Xn);
#endif
}
else
{
startnibble = 1; /* Subsequent frames: skip nibbles */
#if DECODE_DEBUG
ms_log (0, "Frame %d\n", frameidx);
#endif
}
/* Swap 32-bit word containing the nibbles */
if (swapflag)
ms_gswap4 (&frame[0]);
/* Decode each 32-bit word according to nibble */
for (widx = startnibble; widx < 16; widx++)
{
/* W0: the first 32-bit quantity contains 16 x 2-bit nibbles (high order bits) */
nibble = EXTRACTBITRANGE (frame[0], (30 - (2 * widx)), 2);
switch (nibble)
{
case 0: /* nibble=00: Special flag, no differences */
#if DECODE_DEBUG
ms_log (0, " W%02d: 00=special\n", widx);
#endif
break;
case 1: /* nibble=01: Four 8-bit differences, starting at high order bits */
word = (union dword *)&frame[widx];
for (idx = 0; idx < 4; idx++)
{
diff[diffidx++] = word->d8[idx];
}
#if DECODE_DEBUG
ms_log (0, " W%02d: 01=4x8b %d %d %d %d\n", widx,
diff[diffidx - 4], diff[diffidx - 3], diff[diffidx - 2], diff[diffidx - 1]);
#endif
break;
case 2: /* nibble=10: Must consult dnib, the high order two bits */
if (swapflag)
ms_gswap4 (&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 */
diff[diffidx++] = (s30.x = EXTRACTBITRANGE (frame[widx], 0, 30));
#if DECODE_DEBUG
ms_log (0, " W%02d: 10,01=1x30b %d\n", widx, diff[diffidx - 1]);
#endif
break;
case 2: /* nibble=10, dnib=10: Two 15-bit differences, starting at high order bits */
for (idx = 0; idx < 2; idx++)
{
diff[diffidx++] = (s15.x = EXTRACTBITRANGE (frame[widx], (15 - idx * 15), 15));
}
#if DECODE_DEBUG
ms_log (0, " W%02d: 10,10=2x15b %d %d\n", widx,
diff[diffidx - 2], diff[diffidx - 1]);
#endif
break;
case 3: /* nibble=10, dnib=11: Three 10-bit differences, starting at high order bits */
for (idx = 0; idx < 3; idx++)
{
diff[diffidx++] = (s10.x = EXTRACTBITRANGE (frame[widx], (20 - idx * 10), 10));
}
#if DECODE_DEBUG
ms_log (0, " W%02d: 10,11=3x10b %d %d %d\n", widx,
diff[diffidx - 3], diff[diffidx - 2], diff[diffidx - 1]);
#endif
break;
}
break;
case 3: /* nibble=11: Must consult dnib, the high order two bits */
if (swapflag)
ms_gswap4 (&frame[widx]);
dnib = EXTRACTBITRANGE (frame[widx], 30, 2);
switch (dnib)
{
case 0: /* nibble=11, dnib=00: Five 6-bit differences, starting at high order bits */
for (idx = 0; idx < 5; idx++)
{
diff[diffidx++] = (s6.x = EXTRACTBITRANGE (frame[widx], (24 - idx * 6), 6));
}
#if DECODE_DEBUG
ms_log (0, " W%02d: 11,00=5x6b %d %d %d %d %d\n", widx,
diff[diffidx - 5], diff[diffidx - 4], diff[diffidx - 3], diff[diffidx - 2],
diff[diffidx - 1]);
#endif
break;
case 1: /* nibble=11, dnib=01: Six 5-bit differences, starting at high order bits */
for (idx = 0; idx < 6; idx++)
{
diff[diffidx++] = (s5.x = EXTRACTBITRANGE (frame[widx], (25 - idx * 5), 5));
}
#if DECODE_DEBUG
ms_log (0, " W%02d: 11,01=6x5b %d %d %d %d %d %d\n", widx,
diff[diffidx - 6], diff[diffidx - 5], diff[diffidx - 4], diff[diffidx - 3],
diff[diffidx - 2], diff[diffidx - 1]);
#endif
break;
case 2: /* nibble=11, dnib=10: Seven 4-bit differences, starting at high order bits */
for (idx = 0; idx < 7; idx++)
{
diff[diffidx++] = (s4.x = EXTRACTBITRANGE (frame[widx], (24 - idx * 4), 4));
}
#if DECODE_DEBUG
ms_log (0, " W%02d: 11,10=7x4b %d %d %d %d %d %d %d\n", widx,
diff[diffidx - 7], diff[diffidx - 6], diff[diffidx - 5], diff[diffidx - 4],
diff[diffidx - 3], diff[diffidx - 2], diff[diffidx - 1]);
#endif
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 */
} /* Done looping over nibbles and 32-bit words */
/* Apply differences in this frame to calculate output samples,
* ignoring first difference for first frame */
for (idx = (frameidx == 0) ? 1 : 0;
idx < diffidx && outputidx < samplecount;
idx++, outputidx++)
{
output[outputidx] = output[outputidx - 1] + diff[idx];
}
} /* Done looping over frames */
/* Check data integrity by comparing last sample to Xn (reverse integration constant) */
if (outputidx == samplecount && output[outputidx - 1] != Xn)
{
ms_log (1, "%s: Warning: Data integrity check for Steim2 failed, Last sample=%d, Xn=%d\n",
srcname, output[outputidx - 1], Xn);
}
return (int)outputidx;
} /* 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.
*
* \ref MessageOnError - this function logs a message on error
************************************************************************/
int
msr_decode_geoscope (char *input, int64_t samplecount, float *output,
int64_t outputlength, int encoding,
const 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, "%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_gswap2 (&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_gswap2 (&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, int64_t samplecount, int32_t *output,
int64_t 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_gswap2 (&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, int64_t samplecount, int32_t *output,
int64_t outputlength, const 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_gswap2 (&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, "%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, int64_t samplecount, int32_t *output,
int64_t 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_gswap2 (&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() */
Reference in New Issue View Git Blame Copy Permalink