A low-footprint compression technique for when efficiency really matters.
March 01, 2003
URL:http://www.drdobbs.com/easy-analog-data-compression/184401624
A low-footprint compression technique for when efficiency really matters. When I needed to implement a lossless data compressor for 16-bit ECGs (electrocardiograms), EEGs (electroencephalograms), and other medical signals in C on a 16-bit integer processor, I discovered how little value the majority of compression techniques were for my problem. Nearly all of the techniques I found on the Internet and in computer science literature rely on two fundamental techniques. The first technique consists of repetitive patterns (known as the Lempel-Ziv family of algorithms, such as gzip) that can be encoded in an efficient way. The second technique uses a family of approaches (known as entropy coding or Huffman coding) that exploit the fact that the individual numbers/letters found in the raw data are usually not equally probable. Both techniques make use of look-up tables -- the larger the tables, the better the chance for good compression ratios. That was another challenge in my project. My algorithm should use less than 2KB of data memory for up to four independent channels and should not need more than 1,000 instruction words. I needed something very simple but effective.
To the untrained human eye, an ECG appears to be extremely repetitive (see Figure 1). Each heartbeat yields a waveform that looks very much like the preceding heartbeat, but the data values differ slightly. The barely visible differences of the high-resolution signals make the compression with standard techniques nearly impossible. In an experiment, I assessed that a standard version of gzip reduces the size of my ECG by barely five percent although I consider myself quite healthy. An EEG curve (see Figure 2) is so irregular that zip, pack, and others do an even worse job.
However, I know something that Lempel, Ziv, and Huffman do not. My data consists of digitized signals originating from an analog source. These signals are known to exhibit relatively few big jumps in their signal amplitude. For signals that change slightly over time, "delta encoding" (i.e., the coding of the signed differences between two sample values instead of the sample values themselves), saves a lot of memory. Suppose, for example, you know that the absolute difference between two adjacent samples of the 16-bit signal never exceeds 1,023; you can code the signal as one absolute value of 16 bits followed by the delta values of 11 bits each.
If a new delta value cannot be represented by the current bit width, then iBitWidthOfCode must be increased. This time the decoder cannot anticipate the change and must be informed. The encoder does this by storing a reserved letter piOVFLCODE[] that exists for every possible delta bit width. When the decoder encounters such a reserved letter, it knows it must increase its iBitWidthOfCode by one. To reduce excessive coding of piOVFLCODE[] when iBitWidthOfCode gets so small that merely signal noise is coded, I have introduced the lower bound MIN_BIT_WIDTH.
I also have reserved two other letters: piALIGNCODE[] indicates that the next code is aligned with the next word boundary, an indication I need to implement flushing; piABSCODE[] indicates that an absolute code (a raw sample value) follows, aligned with word boundaries.
These reserved letters are chosen to flank the largest and smallest possible two's complement numbers of the delta code. If a delta value matches a reserved letter, the bit width must be increased before the value can be encoded without ambiguities. Hence, I use the reserved letters in storeDelta16 to check if iBitWidthOfCode must be increased. If this is the case, then the appropriate piOVFLCODE[] is stored, and storeDelta16() is called recursively with an increased bit width.
The function storeX, called in storeDelta16, concatenates the codes of variable bit lengths that are passed as parameters. As soon as a WORD (16 bits) is complete, it passes it to writeCode. Any code fraction that is cut off when completing a WORD is kept in the buffer iXcode. The function writeCode represents the interface of a device driver, which is supposed to accept 16-bit-entities only. In my application, data needs to be transmitted via a serial interface and written to flash memory. An additional requirement is to be able to flush the buffer iXcode at any time. This is possible by calling flushX, which checks if there is anything buffered in iXcode. If this is the case, it will add the aforementioned piALIGNCODE[] to the output and pass the content of iXcode to writeCode, although not all of its 16 bits are meaningful. When the decoding algorithm reads the reserved letter piALIGNCODE[], it knows that any undecoded remaining bits of the WORD at hand can be dismissed, and the next delta code will be found in the next WORD.
However, if you can extrapolate the signal curve and predict the new sample value with reasonable accuracy, then you may encode the difference between the prediction and the actual sample value to further reduce the size of the delta code. The function encoder_storeADEPT shown in Listing 1 performs these compression steps.
I call the compression technique ADEPT (adaptive delta encoding with prediction). The prediction step itself may be chosen according to a model of the signal, heuristics, similarities to other signals recorded simultaneously, etc. Signal prediction is a complicated matter that could fill entire books or at least book chapters [1]. I decided to implement a simple linear extrapolator instead of spending the rest of my days trying to understand those books. The predictor given here records the history of the last three samples, tries to fit a line through them, and extrapolates linearly.
As you can see in encoder_storeADEPT, I have arranged the computations of the predictor such that first I do the division and then the multiplication. This makes an integer overflow less likely in case you rewrite the function to deal with 32-bit data. If the prediction fails completely, it may happen that the generated delta code cannot be represented as a 16-bit number anymore. Then the algorithm resets the predictor, writes the reserved letter piABSCODE[], and finally writes the sample directly, word aligned, without delta encoding. If no absolute code is needed, encoder_storeADEPT passes the delta value to storeDelta16.
Note that the simple predictor presented above assumes a rather continuous signal. It improves the compression for EEG signals or breathing curves, but performs very poorly for ECGs. Having no prediction step at all is sometimes better than having an inappropriate one: heartbeats of a very high amplitude may irritate the presented predictor so much that pure adaptive delta encoding outperforms adaptive delta encoding with prediction.
Decode does this division into delta values in a while loop. As a bitmask for this job, it uses mpiABSCODE[]. The obtained delta code iSubword might represent a negative number, but it might be a positive 32-bit two's complement number and therefore must be sign extended. Think of the number 0xFF. In an 8-bit two's complement representation, it is interpreted as -1. It must be negative because the MSB (most significant bit) is 1. In a 32-bit system, this is interpreted as the positive number 255. Its MSB is 0. To interpret the number correctly in the 32-bit system, the bits to the left of the sign bit of the small-scale representation must hold a copy of this sign bit. In our example, 0xFF will be converted to 0xFFFFFFFF, which is interpreted as -1. This technique is called sign extension.
After the sign extension of iSubword, Decode checks if the obtained code is a reserved letter. If so, it changes the object's state to be able to react appropriately; else it reconstructs the original data by calling the member function DeltaDecode and stores the result in the vector named by the constructor.
Listing 3 shows an example program that uses ADEPT compression and decompression.
ADEPT performs well if the difference between two subsequent samples is small and/or the signal form is fairly linear. In the ideal case (i.e., when a strictly linear signal is compressed), ADEPT exhibits a compression ratio of 16:MIN_BIT_WIDTH. With a reasonable value for the minimal bit width, it is obvious that the best compression ratio that can be achieved with ADEPT is between 16:3 and 16:4, or 18 and 25percent, respectively.
Obviously a noise-free, repetitive signal can be more compressed by zip and pack than by ADEPT, which cannot go beyond the mentioned threshold. However, signals coming from natural sources are rarely exactly repetitive, therefore the aforementioned bad performance of gzip on an ECG. However, a comparison of ADEPT to gzip is unfair. Each of the algorithms has a very different application domain. Hence the attempt to compress an EEG with gzip must fail, and you should compare ADEPT to other techniques.
The strengths of ADEPT are most of all simplicity, no need for floating-point arithmetic, small memory requirements, and zero delay. The latter needs explanation: compression techniques based on so-called FIR filters inherently introduce a time delay. This may be very annoying when real-time performance is desired.
An obvious weakness of the technique is its dependency on scale (i.e., on the gain factor of the analog/digital converter); a higher gain value will lead to larger delta values for the digitized data and thus needs a wider code. Simplicity pays its price when you compare ADEPT to specialized algorithms, as shown in Table 1. Note that such comparisons are very difficult, because they are hardly ever objective. This is especially true when different input data is used to assess the performance of the respective methods, which is the case here.
Changes of the signal scale may increase or decrease the compression ratio of ADEPT given in Table 1 by about five percent. Most likely, the same is true for the two reference ECG compressors, but no figures were given in [2]. The EEG compressor in [3] is probably the most powerful known, but it is so specialized that it cannot encode arbitrary signals without losses.
A nice feature of the algorithms ALZ77 and ADPCM Subband Coding is a smooth transition from lossless to lossy compression, a step that is difficult with DPCM, ADEPT, or entropy coding. When tuning ADEPT, you should note the comparatively little compression gain when improving the prediction because of the logarithmic dependency between the delta code width and the prediction error.
[2] R. Nigel Horspool and Warren J. Windels. "An LZ Approach to ECG Compression," Proceedings of the 1994 IEEE 7th Symposium on Computer-Based Medical Systems. See also <www.csr.uvic.ca/~nigelh/Publications/ECG-compression.pdf>.
[3] Zlatko Sijercic, Gyan Agarwal, and Charles Anderson. "EEG Signal Compression with ADPCM Subband Coding," Proceedings of the 39th Midwest Symposium on Circuits and Systems, Aug 1996. See also <www.cs.colostate.edu/~anderson/res/eeg/>.
[4] John G. Proakis. Digital Communications, 3rd Edition (McGraw-Hill, 1996), Chapter 3.
Figure 1: An ECG signal
Figure 2: An EEG signal
Listing 1: Partial listing of the encoder/compressor encoder.c
static iBitWidthOfCode = 6; static INT16 iXcode = 0; static int iXbitsTaken = 0; /* that many bits of iXcode are valid */ static const INT16 piABSCODE[17] = { +0x0000, +0x0001, +0x0003, +0x0007, +0x000F, +0x001F, +0x003F, +0x007F, +0x00FF, +0x01FF, +0x03FF, +0x07FF, +0x0FFF, +0x1FFF, +0x3FFF, +0x7FFF, (INT16) +0xFFFF }; static const INT16 piOVFLCODE[16] = { +0x0000, +0x0000, +0x0002, +0x0006, +0x000E, +0x001E, +0x003E, +0x007E, +0x00FE, +0x01FE, +0x03FE, +0x07FE, +0x0FFE, +0x1FFE, +0x3FFE, +0x7FFE }; static const INT16 piALIGNCODE[16] = { -0x0001, -0x0002, -0x0004, -0x0008, -0x0010, -0x0020, -0x0040, -0x0080, -0x0100, -0x0200, -0x0400, -0x0800, -0x1000, -0x2000, -0x4000, -0x8000 }; static INT32 piLastSamplesDiv3[3] = {0, 0, 0}; void writeCode(int iSample) { flash_write((WORD) iSample); } void storeX(INT16 iCode) { int iEmpty = 16 - iXbitsTaken; /* available bits in the word */ /* storeX() expects a iBitWidthOfCode-bit value; * however, negative numbers are sign extended up to the 16-bit MSB; * this sign extension must thus be removed */ iCode &= piABSCODE[iBitWidthOfCode]; if (iEmpty > iBitWidthOfCode) /* more bits available than needed? */ { /* yes: write entire bit pattern */ iCode <<= iXbitsTaken; /* position new sub-word correctly */ iXcode |= iCode; /* keep it to -or- w. next pattern */ iXbitsTaken += iBitWidthOfCode; } else if (iEmpty == iBitWidthOfCode) /* if exactly as much bits are */ { /* available as needed, then */ iCode <<= iXbitsTaken; /* position new sub-word, */ iCode |= iXcode; /* -or- new sample to kept val */ iXcode = 0x0000; /* reset the kept values */ iXbitsTaken = 0; writeCode(iCode); /* store new sample + kept val.*/ } else { /* we have to split the bit pattern into two words */ register WORD wTemp = iCode; /* make a copy */ iCode <<= iXbitsTaken; /* position new sub-word */ iCode |= iXcode; /* create the 1st part */ wTemp >>= iEmpty; /* position part two */ iXcode = wTemp; /* and keep it for next write */ iXbitsTaken = iBitWidthOfCode - iEmpty; writeCode(iCode); /* store 1st part */ } } void flushX() { /* Check if the buffer of storeX() is not empty. */ if (iXbitsTaken != 0) { /* record that the next sample that follows will be aligned at * WORD boundaries */ storeX(piALIGNCODE[iBitWidthOfCode - 1]); /* check if the ALIGNCODE already managed the alignment */ if (iXbitsTaken != 0) { /* if not, then flush the buffer and enforce the alignment */ writeCode(iXcode); } iXcode = 0; iXbitsTaken = 0; } } /* flushX() */ void encoder_flush() { flushX(); } void storeDelta16(INT16 iDeltaCode) { static iSamples2MSBidentical = 0; if ((iDeltaCode >= piOVFLCODE[iBitWidthOfCode - 1]) || (iDeltaCode <= piALIGNCODE[iBitWidthOfCode - 1])) { /* overflow/underflow: store OVFLCODE and increase bit width */ storeX(piOVFLCODE[iBitWidthOfCode - 1]); ++iBitWidthOfCode; assert(iBitWidthOfCode < 17); iSamples2MSBidentical = 0; storeDelta16(iDeltaCode); /* retry with incr. width */ } else { int iMaxBitMask = 3 << (iBitWidthOfCode - 2); int iDeltaBitWidthAfterStoring = 0; if ((iDeltaCode & iMaxBitMask) == 0x0000 || (iDeltaCode & iMaxBitMask) == iMaxBitMask) { /* the two MSBits of iSample are identical */ iSamples2MSBidentical++; if ((iSamples2MSBidentical == BIT_REDUCTION_LENGTH) && (iBitWidthOfCode > MIN_BIT_WIDTH)) { /* already the BIT_REDUCTION_LENGTH-th sample where * MSBit is not in use --> reduce the code width */ iDeltaBitWidthAfterStoring = -1; iSamples2MSBidentical = 0; } } else { iSamples2MSBidentical = 0; } storeX(iDeltaCode); iBitWidthOfCode += iDeltaBitWidthAfterStoring; } } /* storeDelta16() */ void encoder_storeADEPT(INT16 iSample) { static BOOL bFirstValue = TRUE; /* start-up needed */ INT32 iCodedSample; /* the delta code */ /* estimate the value x'(t) for iSample out of the last * iSamples x(t-1), x(t-2), x(t-3) with the formula * x'(t) = -2/3*x(t-3) + 1/3*x(t-2) + 4/3*x(t-1). * This value x'(t) is taken as basis for the delta encoding */ INT32 iPrediction; /* linear extrapolation */ iPrediction = 4 * piLastSamplesDiv3[2] + piLastSamplesDiv3[1] - 2 * piLastSamplesDiv3[0]; iCodedSample = iSample - iPrediction; /* delta encoding */ piLastSamplesDiv3[0] = piLastSamplesDiv3[1]; piLastSamplesDiv3[1] = piLastSamplesDiv3[2]; piLastSamplesDiv3[2] = iSample/3; /* div here to prevent overflow */ if ((iCodedSample >= piOVFLCODE[15]) || (iCodedSample <= piALIGNCODE[15]) || (bFirstValue)) { /* signal that a 32-bit absolute value follows */ storeX(piABSCODE[iBitWidthOfCode - 1]); /* check alignment; if not aligned, * flush the remaining bits of the buffer */ if (iXbitsTaken != 0) { writeCode(iXcode); iXbitsTaken = 0; iXcode = 0x0000; } /* reset delta encoding and store the absolute value */ piLastSamplesDiv3[0] = iSample/3; piLastSamplesDiv3[1] = iSample/3; piLastSamplesDiv3[2] = iSample/3; writeCode(iSample); bFirstValue = FALSE; } else { storeDelta16((INT16) iCodedSample); } }
Listing 2: Partial listing of decoder.cpp
class TDecoder { private: TDecoder(); // no default constructor wanted public: TDecoder(vector<int> &tiDestination); // init the TDecoder and pass a vector whereto decoded data shall // be written; public: void Decode(WORD wCode); // decode the word wCode and push_back the result into tiDestination private: int DeltaDecode(int iDelta); int miXsample; // the untouched fraction of the last 16-bit code passed to Decode() int miXbitsTaken; // denotes what number of bits are valid in miXsample int miBitWidthOfCode; // holds the bit width of the code int miAbsValueCoded; // this holds the number of WORDs left that code an absolute value int mpiLastSamplesDiv3[3]; // usually holds the last 3 decoded samples int miSamples2MSBidentical; // this counter is used for the reduction of the bit width (in // the adaption) vector<int> *mptiDestination; // pointer to the result vector static const int mpiABSCODE[17]; static const int mpiOVFLCODE[16]; static const int mpiALIGNCODE[16]; static const int mpiSIGNEXTENSION[16]; }; const int TDecoder::mpiABSCODE[17] = { ... }; // as in listing 1 const int TDecoder::mpiOVFLCODE[16] = { ... }; // as in listing 1 const int TDecoder::mpiALIGNCODE[16] = { ... }; // as in listing 1 const int TDecoder::mpiSIGNEXTENSION[16] = { // the value at index i denotes a bitmask for 32-bits sign extension // for an i-bits value 0xFFFFFFFE, 0xFFFFFFFC, 0xFFFFFFF8, 0xFFFFFFF0, 0xFFFFFFE0, 0xFFFFFFC0, 0xFFFFFF80, 0xFFFFFF00, 0xFFFFFE00, 0xFFFFFC00, 0xFFFFF800, 0xFFFFF000, 0xFFFFE000, 0xFFFFC000, 0xFFFF8000, 0xFFFF0000 }; TDecoder::TDecoder(vector<int> &tiDestination): miXsample(0), miXbitsTaken(0), miBitWidthOfCode(6), miAbsValueCoded(0), miSamples2MSBidentical(0), mptiDestination(&tiDestination) { assert( sizeof(int) >= 4 ); // TDecoder does not run on 16-bit // platforms for(int i = 0; i < 3; i++) mpiLastSamplesDiv3[i] = 0; } int TDecoder::DeltaDecode(int iDelta) { int iPrediction = 4 * mpiLastSamplesDiv3[2] + mpiLastSamplesDiv3[1] - 2 * mpiLastSamplesDiv3[0]; int iDecodedSample = iPrediction + iDelta; mpiLastSamplesDiv3[0] = mpiLastSamplesDiv3[1]; mpiLastSamplesDiv3[1] = mpiLastSamplesDiv3[2]; mpiLastSamplesDiv3[2] = iDecodedSample / 3; return iDecodedSample; } void TDecoder::Decode(WORD wCode) { unsigned uiCode = wCode; // use full 32-bit register width if( miAbsValueCoded == 1 ) { miXsample = uiCode; if (uiCode & 0x8000) miXsample -=0x10000; // unsigned -> int16 // conversion (*mptiDestination).push_back(miXsample); mpiLastSamplesDiv3[0] = miXsample/3; // update for // delta encoding mpiLastSamplesDiv3[1] = miXsample/3; mpiLastSamplesDiv3[2] = miXsample/3; miAbsValueCoded--; // switch back to delta // encoding miXsample = 0x0000; // reset buffered sub-word } else { // a delta encoded value, not necessarily aligned with a word // boundary has been received uiCode <<= miXbitsTaken; // shift it left to and combine it with miXsample |= uiCode; // a buffered rest of preceding iCode:s miXbitsTaken += 16; // 16 new bits to decode were passed while( miXbitsTaken >= miBitWidthOfCode ) { // cut out a complete code; reduce contents of the buffer by // the code cut out; reduce number of valid bits accordingly int iSubword = miXsample & mpiABSCODE[miBitWidthOfCode]; int iSignBitMask = 1 << (miBitWidthOfCode - 1); miXsample >>= miBitWidthOfCode; miXbitsTaken -= miBitWidthOfCode; int iMaxBitMask = 3 << (miBitWidthOfCode - 2); if( iSignBitMask & iSubword) { // a negative number was coded; sign extension is needed // as well for the comparison to mpiALIGNCODE[] as for // delta decoding iSubword |= mpiSIGNEXTENSION[miBitWidthOfCode - 1]; } if( iSubword == mpiOVFLCODE[miBitWidthOfCode - 1] ) { miBitWidthOfCode++; // increase bit width of code miSamples2MSBidentical = 0; } else if( iSubword == mpiABSCODE[miBitWidthOfCode - 1] ) { miAbsValueCoded = 1; miXbitsTaken = 0; miXsample = 0; } else if( iSubword == mpiALIGNCODE[miBitWidthOfCode - 1] ) { miXbitsTaken = 0; miXsample = 0; } else { if( (iSubword & iMaxBitMask) == 0x0000 || (iSubword & iMaxBitMask) == iMaxBitMask) { miSamples2MSBidentical++; if ((miSamples2MSBidentical == BIT_REDUCTION_LENGTH) && (miBitWidthOfCode > MIN_BIT_WIDTH)) { miBitWidthOfCode--; miSamples2MSBidentical = 0; } } else { miSamples2MSBidentical = 0; } int iNewSample = DeltaDecode( iSubword ); (*mptiDestination).push_back( iNewSample ); } } // while } } // TDecoder::Decode()
Listing 3: An example program using ADEPT compression and decompression
void test() { ifstream nIntFile("IntegerSignal.txt"); int iValue; vector<int> tiSignal; // load the test signal while(nIntFile >> iValue) { tiSignal.push_back(iValue); } for (vector<int>::iterator tiIt = tiSignal.begin(); tiIt != tiSignal.end(); tiIt++) { encoder_storeADEPT(*tiIt); } encoder_flush(); WORD wCoded; vector<int> tiReconstruction; TDecoder tChn0Decoder(tiReconstruction); // decode what is read from Flash and store it in tiReconstruction while(flash_read(wCoded)) { tChn0Decoder.Decode(wCoded); } assert(tiReconstruction == tiSignal); }
Table 1: Comparison of ADEPT to lossless ECG and EEG compressors
Terms of Service | Privacy Statement | Copyright © 2024 UBM Tech, All rights reserved.