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