（13，17）卷积码_咬尾卷积码_咬尾卷积_

%% PURPOSE: % Implement the Log-MAP algorithm for binary-convolutional-codes invented by Bahl etal. Computes the A-Posteriori-Probability % of each i/p bit being 1 or 0 given the received sequence. The channel must be Discrete-Memoryless. % %% INPUTS: % TREL = matlab trellis structure from poly2trellis() % RX = real-valued received bit LLRs. RX can be linearly-quantized. %% USAGE: %% HISTORY: % 10-Aug-11: initial version % 15-Aug-11: verified BER=0 on noiseless case for all codes. %% Author: Li Qiang % Copyright 2011. All Rights Reserved function [LLRu, LLRc] = LogMAPdecode(TREL, RX) if 0 return; end EncoderN = log2(TREL.numOutputSymbols); % cardinality of set of all possible o/p bitvecs EncoderK = log2(TREL.numInputSymbols); % cardinality of set of all possible i/p bitvecs RX = reshape(RX, 1, numel(RX)); % Make RX a row NumBran = length(RX)/EncoderN; % number of branches sent by TX to RX if rem(length(RX), EncoderN) ~= 0, disp('ERROR: Make length(RX) = multiple of EncoderN') end %% Pre-compute {contributing prior-states, corresponding i/p bitvec, branch %% o/p bitvec} for each current-state. Reqd for forward-recursion. This is %% independent of RX. prevStates = NaN(TREL.numStates, TREL.numInputSymbols); prevStateIn = NaN(TREL.numStates, TREL.numInputSymbols); prevStateOut = NaN(TREL.numStates, TREL.numInputSymbols); for dest_st = 0 : (TREL.numStates -1) % Find all origin-states from-which there could be a transition into dest_st % and store in row-(dest_st +1) [row, col] = find(TREL.nextStates == dest_st); prevStates(dest_st +1, :) = row' -1; % origin-state of transition prevStateIn(dest_st +1, :) = col' -1; % Corresponding i/p bit-vec causing transition prevStateOut(dest_st +1, :) = TREL.outputs( sub2ind(size(TREL.nextStates), row, col) )'; % Corresponding o/p branch bit-vec (in octal) end %% Compute branch-metrics gamma*(m', m) fev branch fev time using segments of RX. There %% are numOutputSymbols distinct branches. col-i contains all possible %% branch-metrics at time i, i= 1,...,NumBran. BranMetric = NaN(TREL.numOutputSymbols, NumBran); for time = 1 : NumBran, rxvec = RX((time-1)*EncoderN +1 : time*EncoderN); % [1, N], [N+1, 2*N] , ... for BranIdx = 1 : TREL.numOutputSymbols bitvec = bitget(BranIdx-1, EncoderN : -1 : 1); % All possible binary-vectors of length EncoderN BranMetric(BranIdx, time) = (1-bitvec)*rxvec'; % metric at "time" for branch = binary-representation(BranIdx-1) end end %% Alpha*(t, m) = sum_{m'=0,...,NS-1} [ Alpha*(t-1, m') * gamma*(m', m)] %% Initialize Alpha*(time=0) as below, since TX ensures encoder starts trellis at state=0 at time t=0. Alpha = NaN(NumBran +1, TREL.numStates); Alpha(1, 1) = 0; Alpha(1, 2:end) = -1e9; for time = 1 : NumBran, for dest_st = 0 : (TREL.numStates -1) sum1 = -1e9; for idx = 1 : length( prevStates(dest_st+1, :) ) orig_s = prevStates(dest_st+1, idx); enc_out_bitv = prevStateOut(dest_st+1, idx); sum1 = ln_of_sum_of_exps( sum1, Alpha(time -1 +1, orig_s +1) + BranMetric(oct2dec(enc_out_bitv) +1, time) ); end Alpha(time+1, dest_st+1) = sum1; end Alpha(time+1, :) = Alpha(time+1, :)-max(Alpha(time+1, :)); end %% Beta*(t-1, m) = sum_{m'=0,...,NS-1} [ Beta*(t, m') * gamma*(m', m)] %% Initialize Beta*(time= NumBran) as below, since TX ensures encoder terminates trellis at state=0 at time t=NumBran Beta = NaN(NumBran +1, TREL.numStates); Beta(NumBran +1, 1) = -log(TREL.numStates); Beta(NumBran +1, 2:end) = -log(TREL.numStates); for time = (NumBran -1) : -1 : 0 for orig_st = 0 : (TREL.numStates -1) sum2 = -1e9; for idx = 1 : length( TREL.nextStates(orig_st+1, :) ) dest_st = TREL.nextStates(orig_st+1, idx); enc_out_bitv = TREL.outputs(orig_st+1, idx); sum2 = ln_of_sum_of_exps( sum2, Beta(time +2, dest_st +1) + BranMetric(oct2dec(enc_out_bitv) +1, time+1) ); end Beta(time+1, orig_st+1) = sum2; end Beta(time+1,:) = Beta(time+1,:) - max(Beta(time+1,:)); end %% Final LLR computation for information bits for each time %% L(U(k)|Y) = ln( . / .) LLRu = zeros(1,NumBran); for time = NumBran : -1 : 1 %% prob of all transitions arising from i/p bit = 0. Works only for EncoderK = 1. total_prob_bit0 = -1e9; for orig_st = 0 : (TREL.numStates -1) dest_st = TREL.nextStates(orig_st+1, 1); branch = oct2dec(TREL.outputs(orig_st+1, 1)); incr_prob = Alpha(time, orig_st+1) + Beta(time+1, dest_st+1) + BranMetric(branch+1, time); total_prob_bit0 = ln_of_sum_of_exps( total_prob_bit0, incr_prob); end %% prob of all transitions arising from i/p bit = 1. Works only for EncoderK = 1. total_prob_bit1 = -1e9; for orig_st = 0 : (TREL.numStates -1) dest_st = TREL.nextStates(orig_st+1, 2); branch = oct2dec(TREL.outputs(orig_st+1, 2)); incr_prob = Alpha(time, orig_st+1) + Beta(time+1, dest_st+1) + BranMetric(branch+1, time); total_prob_bit1 = ln_of_sum_of_exps( total_prob_bit1, incr_prob); end LLRu(time) = total_prob_bit0 - total_prob_bit1; end %% Final LLR computation for code bits for each time LLRc = zeros(1, EncoderN*NumBran); for time = NumBran : -1 : 1 total_prob_bit = -1e9*ones(1,TREL.numOutputSymbols); for code_sym = 0 : TREL.numOutputSymbols - 1 [row,col]=find(oct2dec(TREL.outputs) == code_sym); for st_idx = 1 : length(row) dest_state = TREL.nextStates(row(st_idx), col(st_idx)); incr_prob = Alpha(time, row(st_idx)) + Beta(time+1, dest_state+1) + BranMetric(code_sym+1, time); total_prob_bit(code_sym+1) = ln_of_sum_of_exps( total_prob_bit(code_sym+1), incr_prob); end end % total_prob_bit = total_prob_bit - log(sum(exp(total_prob_bit))); % Normalization total_prob_bit = total_prob_bit - max(total_prob_bit); % Normalization %% symbol LLR->bit LLR, works only for EncoderN=2 fst_bit0 = ln_of_sum_of_exps(total_prob_bit(1), total_prob_bit(2)); fst_bit1 = ln_of_sum_of_exps(total_prob_bit(3), total_prob_bit(4)); LLRc(2*time-1) = fst_bit0 - fst_bit1; %% for the second bit in a code symbol sec_bit0 = ln_of_sum_of_exps(total_prob_bit(1), total_prob_bit(3)); sec_bit1 = ln_of_sum_of_exps(total_prob_bit(2), total_prob_bit(4)); LLRc(2*time) = sec_bit0 - sec_bit1; end %限幅 NUMMAX=20; for tt=1:length(LLRu) if abs(LLRu(tt))>NUMMAX LLRu(tt)=sign(LLRu(tt))*NUMMAX; end end for tt=1:length(LLRc) if abs(LLRc(tt))>NUMMAX LLRc(tt)=sign(LLRc(tt))*NUMMAX; end end return; %%========================================================================= %% Function ln( e^x + e^y) = max(x, y) + ln(1 + e^-|x-y|) function a = ln_of_sum_of_exps(x, y) a = max(x, y) + log( 1 + exp(-abs(x - y)) ); return;