Roman Khassraf | aa8fa99 | 2015-06-02 09:20:03 +0200 | [diff] [blame] | 1 | /* |
| 2 | * Copyright 2008, 2009, 2014 Free Software Foundation, Inc. |
| 3 | * Copyright 2014 Range Networks, Inc. |
| 4 | * |
| 5 | * |
| 6 | * This software is distributed under the terms of the GNU Affero Public License. |
| 7 | * See the COPYING file in the main directory for details. |
| 8 | * |
| 9 | * This use of this software may be subject to additional restrictions. |
| 10 | * See the LEGAL file in the main directory for details. |
| 11 | |
| 12 | This program is free software: you can redistribute it and/or modify |
| 13 | it under the terms of the GNU Affero General Public License as published by |
| 14 | the Free Software Foundation, either version 3 of the License, or |
| 15 | (at your option) any later version. |
| 16 | |
| 17 | This program is distributed in the hope that it will be useful, |
| 18 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 19 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 20 | GNU Affero General Public License for more details. |
| 21 | |
| 22 | You should have received a copy of the GNU Affero General Public License |
| 23 | along with this program. If not, see <http://www.gnu.org/licenses/>. |
| 24 | |
| 25 | */ |
| 26 | |
| 27 | |
| 28 | |
| 29 | |
| 30 | #include "BitVector.h" |
| 31 | #include "ViterbiR204.h" |
| 32 | #include <iostream> |
| 33 | #include <stdio.h> |
| 34 | #include <sstream> |
| 35 | #include <string.h> |
| 36 | |
| 37 | using namespace std; |
| 38 | |
| 39 | |
| 40 | /** |
| 41 | Apply a Galois polymonial to a binary seqeunce. |
| 42 | @param val The input sequence. |
| 43 | @param poly The polynomial. |
| 44 | @param order The order of the polynomial. |
| 45 | @return Single-bit result. |
| 46 | */ |
| 47 | unsigned ViterbiBase::applyPoly(uint64_t val, uint64_t poly, unsigned order) |
| 48 | { |
| 49 | uint64_t prod = val & poly; |
| 50 | unsigned sum = prod; |
| 51 | for (unsigned i=1; i<order; i++) sum ^= prod>>i; |
| 52 | return sum & 0x01; |
| 53 | } |
| 54 | |
| 55 | unsigned ViterbiBase::applyPoly(uint64_t val, uint64_t poly) |
| 56 | { |
| 57 | uint64_t prod = val & poly; |
| 58 | prod = (prod ^ (prod >> 32)); |
| 59 | prod = (prod ^ (prod >> 16)); |
| 60 | prod = (prod ^ (prod >> 8)); |
| 61 | prod = (prod ^ (prod >> 4)); |
| 62 | prod = (prod ^ (prod >> 2)); |
| 63 | prod = (prod ^ (prod >> 1)); |
| 64 | return prod & 0x01; |
| 65 | } |
| 66 | |
| 67 | |
| 68 | |
| 69 | //void BitVector::encode(const ViterbiR2O4& coder, BitVector& target) |
| 70 | void ViterbiR2O4::encode(const BitVector& in, BitVector& target) const |
| 71 | { |
| 72 | const ViterbiR2O4& coder = *this; |
| 73 | size_t sz = in.size(); |
| 74 | |
| 75 | assert(sz*coder.iRate() == target.size()); |
| 76 | |
| 77 | // Build a "history" array where each element contains the full history. |
| 78 | uint32_t history[sz]; |
| 79 | uint32_t accum = 0; |
| 80 | for (size_t i=0; i<sz; i++) { |
| 81 | accum = (accum<<1) | in.bit(i); |
| 82 | history[i] = accum; |
| 83 | } |
| 84 | |
| 85 | // Look up histories in the pre-generated state table. |
| 86 | char *op = target.begin(); |
| 87 | for (size_t i=0; i<sz; i++) { |
| 88 | unsigned index = coder.cMask() & history[i]; |
| 89 | for (unsigned g=0; g<coder.iRate(); g++) { |
| 90 | *op++ = coder.stateTable(g,index); |
| 91 | } |
| 92 | } |
| 93 | } |
| 94 | |
| 95 | |
| 96 | ViterbiR2O4::ViterbiR2O4() |
| 97 | { |
| 98 | assert(mDeferral < 32); |
| 99 | // (pat) The generator polynomials are: G0 = 1 + D**3 + D**4; and G1 = 1 + D + D**3 + D**4 |
| 100 | mCoeffs[0] = 0x019; // G0 = D**4 + D**3 + 1; represented as binary 11001, |
| 101 | mCoeffs[1] = 0x01b; // G1 = + D**4 + D**3 + D + 1; represented as binary 11011 |
| 102 | computeStateTables(0); |
| 103 | computeStateTables(1); |
| 104 | computeGeneratorTable(); |
| 105 | } |
| 106 | |
| 107 | |
| 108 | void ViterbiR2O4::initializeStates() |
| 109 | { |
| 110 | for (unsigned i=0; i<mIStates; i++) vitClear(mSurvivors[i]); |
| 111 | for (unsigned i=0; i<mNumCands; i++) vitClear(mCandidates[i]); |
| 112 | } |
| 113 | |
| 114 | |
| 115 | |
| 116 | // (pat) The state machine has 16 states. |
| 117 | // Each state has two possible next states corresponding to 0 or 1 inputs to original encoder. |
| 118 | // which are saved in mStateTable in consecutive locations. |
| 119 | // In other words the mStateTable second index is ((current_state <<1) + encoder_bit) |
| 120 | // g is 0 or 1 for the first or second bit of the encoded stream, ie, the one we are decoding. |
| 121 | void ViterbiR2O4::computeStateTables(unsigned g) |
| 122 | { |
| 123 | assert(g<mIRate); |
| 124 | for (unsigned state=0; state<mIStates; state++) { |
| 125 | // 0 input |
| 126 | uint32_t inputVal = state<<1; |
| 127 | mStateTable[g][inputVal] = applyPoly(inputVal, mCoeffs[g], mOrder+1); |
| 128 | // 1 input |
| 129 | inputVal |= 1; |
| 130 | mStateTable[g][inputVal] = applyPoly(inputVal, mCoeffs[g], mOrder+1); |
| 131 | } |
| 132 | } |
| 133 | |
| 134 | void ViterbiR2O4::computeGeneratorTable() |
| 135 | { |
| 136 | for (unsigned index=0; index<mIStates*2; index++) { |
| 137 | mGeneratorTable[index] = (mStateTable[0][index]<<1) | mStateTable[1][index]; |
| 138 | } |
| 139 | } |
| 140 | |
| 141 | |
| 142 | void ViterbiR2O4::branchCandidates() |
| 143 | { |
| 144 | // Branch to generate new input states. |
| 145 | const vCand *sp = mSurvivors; |
| 146 | for (unsigned i=0; i<mNumCands; i+=2) { |
| 147 | // extend and suffix |
| 148 | const uint32_t iState0 = (sp->iState) << 1; // input state for 0 |
| 149 | const uint32_t iState1 = iState0 | 0x01; // input state for 1 |
| 150 | const uint32_t oStateShifted = (sp->oState) << mIRate; // shifted output (by 2) |
| 151 | const float cost = sp->cost; |
| 152 | int bec = sp->bitErrorCnt; |
| 153 | sp++; |
| 154 | // 0 input extension |
| 155 | mCandidates[i].cost = cost; |
| 156 | // mCMask is the low 5 bits, ie, full width of mGeneratorTable. |
| 157 | mCandidates[i].oState = oStateShifted | mGeneratorTable[iState0 & mCMask]; |
| 158 | mCandidates[i].iState = iState0; |
| 159 | mCandidates[i].bitErrorCnt = bec; |
| 160 | // 1 input extension |
| 161 | mCandidates[i+1].cost = cost; |
| 162 | mCandidates[i+1].oState = oStateShifted | mGeneratorTable[iState1 & mCMask]; |
| 163 | mCandidates[i+1].iState = iState1; |
| 164 | mCandidates[i+1].bitErrorCnt = bec; |
| 165 | } |
| 166 | } |
| 167 | |
| 168 | |
| 169 | void ViterbiR2O4::getSoftCostMetrics(const uint32_t inSample, const float *matchCost, const float *mismatchCost) |
| 170 | { |
| 171 | const float *cTab[2] = {matchCost,mismatchCost}; |
| 172 | for (unsigned i=0; i<mNumCands; i++) { |
| 173 | vCand& thisCand = mCandidates[i]; |
| 174 | // We examine input bits 2 at a time for a rate 1/2 coder. |
| 175 | // (pat) mismatched will end up with bits in it for previous transitions, |
| 176 | // but we only use the bottom two bits of mismatched so it is ok. |
| 177 | const unsigned mismatched = inSample ^ (thisCand.oState); |
| 178 | // (pat) TODO: Are these two tests swapped? |
| 179 | thisCand.cost += cTab[mismatched&0x01][1] + cTab[(mismatched>>1)&0x01][0]; |
| 180 | if (mismatched & 1) { thisCand.bitErrorCnt++; } |
| 181 | if (mismatched & 2) { thisCand.bitErrorCnt++; } |
| 182 | } |
| 183 | } |
| 184 | |
| 185 | |
| 186 | void ViterbiR2O4::pruneCandidates() |
| 187 | { |
| 188 | const vCand* c1 = mCandidates; // 0-prefix |
| 189 | const vCand* c2 = mCandidates + mIStates; // 1-prefix |
| 190 | for (unsigned i=0; i<mIStates; i++) { |
| 191 | if (c1[i].cost < c2[i].cost) mSurvivors[i] = c1[i]; |
| 192 | else mSurvivors[i] = c2[i]; |
| 193 | } |
| 194 | } |
| 195 | |
| 196 | |
| 197 | const ViterbiR2O4::vCand& ViterbiR2O4::minCost() const |
| 198 | { |
| 199 | int minIndex = 0; |
| 200 | float minCost = mSurvivors[0].cost; |
| 201 | for (unsigned i=1; i<mIStates; i++) { |
| 202 | const float thisCost = mSurvivors[i].cost; |
| 203 | if (thisCost>=minCost) continue; |
| 204 | minCost = thisCost; |
| 205 | minIndex=i; |
| 206 | } |
| 207 | return mSurvivors[minIndex]; |
| 208 | } |
| 209 | |
| 210 | |
| 211 | const ViterbiR2O4::vCand* ViterbiR2O4::vstep(uint32_t inSample, const float *probs, const float *iprobs, bool isNotTailBits) |
| 212 | { |
| 213 | branchCandidates(); |
| 214 | // (pat) tail bits do not affect cost or error bit count of any branch. |
| 215 | if (isNotTailBits) getSoftCostMetrics(inSample,probs,iprobs); |
| 216 | pruneCandidates(); |
| 217 | return &minCost(); |
| 218 | } |
| 219 | |
| 220 | |
| 221 | void ViterbiR2O4::decode(const SoftVector &in, BitVector& target) |
| 222 | { |
| 223 | ViterbiR2O4& decoder = *this; |
| 224 | const size_t sz = in.size(); |
| 225 | const unsigned oSize = in.size() / decoder.iRate(); |
| 226 | const unsigned deferral = decoder.deferral(); |
| 227 | const size_t ctsz = sz + deferral*decoder.iRate(); |
| 228 | assert(sz <= decoder.iRate()*target.size()); |
| 229 | |
| 230 | // Build a "history" array where each element contains the full history. |
| 231 | // (pat) We only use every other history element, so why are we setting them? |
| 232 | uint32_t history[ctsz]; |
| 233 | { |
| 234 | BitVector bits = in.sliced(); |
| 235 | uint32_t accum = 0; |
| 236 | for (size_t i=0; i<sz; i++) { |
| 237 | accum = (accum<<1) | bits.bit(i); |
| 238 | history[i] = accum; |
| 239 | } |
| 240 | // Repeat last bit at the end. |
| 241 | // (pat) TODO: really? Does this matter? |
| 242 | for (size_t i=sz; i<ctsz; i++) { |
| 243 | accum = (accum<<1) | (accum & 0x01); |
| 244 | history[i] = accum; |
| 245 | } |
| 246 | } |
| 247 | |
| 248 | // Precompute metric tables. |
| 249 | float matchCostTable[ctsz]; |
| 250 | float mismatchCostTable[ctsz]; |
| 251 | { |
| 252 | const float *dp = in.begin(); |
| 253 | for (size_t i=0; i<sz; i++) { |
| 254 | // pVal is the probability that a bit is correct. |
| 255 | // ipVal is the probability that a bit is incorrect. |
| 256 | float pVal = dp[i]; |
| 257 | if (pVal>0.5F) pVal = 1.0F-pVal; |
| 258 | float ipVal = 1.0F-pVal; |
| 259 | // This is a cheap approximation to an ideal cost function. |
| 260 | if (pVal<0.01F) pVal = 0.01; |
| 261 | if (ipVal<0.01F) ipVal = 0.01; |
| 262 | matchCostTable[i] = 0.25F/ipVal; |
| 263 | mismatchCostTable[i] = 0.25F/pVal; |
| 264 | } |
| 265 | |
| 266 | // pad end of table with unknowns |
| 267 | // Note that these bits should not contribute to Bit Error Count. |
| 268 | for (size_t i=sz; i<ctsz; i++) { |
| 269 | matchCostTable[i] = 0.5F; |
| 270 | mismatchCostTable[i] = 0.5F; |
| 271 | } |
| 272 | } |
| 273 | |
| 274 | { |
| 275 | decoder.initializeStates(); |
| 276 | // Each sample of history[] carries its history. |
| 277 | // So we only have to process every iRate-th sample. |
| 278 | const unsigned step = decoder.iRate(); |
| 279 | // input pointer |
| 280 | const uint32_t *ip = history + step - 1; |
| 281 | // output pointers |
| 282 | char *op = target.begin(); |
| 283 | const char *const opt = target.end(); // (pat) Not right if target is larger than needed; should be: op + sz/2; |
| 284 | // table pointers |
| 285 | const float* match = matchCostTable; |
| 286 | const float* mismatch = mismatchCostTable; |
| 287 | size_t oCount = 0; |
| 288 | const ViterbiR2O4::vCand *minCost = NULL; |
| 289 | while (op<opt) { |
| 290 | // Viterbi algorithm |
| 291 | assert(match-matchCostTable<(float)sizeof(matchCostTable)/sizeof(matchCostTable[0])-1); |
| 292 | assert(mismatch-mismatchCostTable<(float)sizeof(mismatchCostTable)/sizeof(mismatchCostTable[0])-1); |
| 293 | minCost = decoder.vstep(*ip, match, mismatch, oCount < oSize); |
| 294 | ip += step; |
| 295 | match += step; |
| 296 | mismatch += step; |
| 297 | // output |
| 298 | if (oCount>=deferral) *op++ = (minCost->iState >> deferral)&0x01; |
| 299 | oCount++; |
| 300 | } |
| 301 | // Dont think minCost == NULL can happen. |
| 302 | mBitErrorCnt = minCost ? minCost->bitErrorCnt : 0; |
| 303 | } |
| 304 | } |
| 305 | |
| 306 | // vim: ts=4 sw=4 |