octave_packages/communications-1.1.1/reedmullerdec.m

   1 ## Copyright (C) 2007, 2011 Muthiah Annamalai <muthiah.annamalai@uta.edu>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {}  reedmullerdec (@var{VV},@var{G},@var{R},@var{M})
  18 ##
  19 ## Decode the received code word @var{VV} using  the RM-generator matrix @var{G},
  20 ## of order @var{R}, @var{M}, returning the code-word C. We use the standard
  21 ## majority logic vote method due to Irving S. Reed. The received word has to be
  22 ## a matrix of column size equal to to code-word size (i.e @math{2^m}). Each row
  23 ## is treated as a separate received word.
  24 ##
  25 ## The second return value is the message @var{M} got from @var{C}.
  26 ##
  27 ## G is obtained from definition type construction of Reed Muller code,
  28 ## of order @var{R}, length @math{2^M}. Use the function reedmullergen,
  29 ## for the generator matrix for the (@var{R},@var{M}) order RM code.
  30 ##
  31 ## Faster code constructions (also easier) exist, but since
  32 ## finding permutation order of the basis vectors, is important, we
  33 ## stick with the standard definitions. To use decoder
  34 ## function reedmullerdec,  you need to use this specific
  35 ## generator function.
  36 ##
  37 ## see: Lin & Costello, Ch.4, "Error Control Coding", 2nd Ed, Pearson.
  38 ##
  39 ## @example
  40 ## @group
  41 ## G=reedmullergen(2,4);
  42 ## M=[rand(1,11)>0.5];
  43 ## C=mod(M*G,2);
  44 ## [dec_C,dec_M]=reedmullerdec(C,G,2,4)
  45 ##
  46 ## @end group
  47 ## @end example
  48 ##
  49 ## @end deftypefn
  50 ## @seealso{reedmullergen,reedmullerenc}
  51
  52 ## FIXME: make possible to use different generators, if permutation
  53 ## matrix (i.e polynomial vector elements of rows of G are given
  54 function [C,CMM]=reedmullerdec(VV,G,R,M)
  55     if ( nargin < 4 )
  56        print_usage();
  57     end
  58
  59    %
  60    % we do a R+1 level majority logic decoding.
  61    % at each order of polynomial modifying the code-word.
  62    %
  63    U=0:M-1; %allowed basis vectors in V2^M.
  64    C=-1*ones(size(VV)); %preset the output word.
  65    [Rows,Cols]=size(G);%rows shadows 'rows()'
  66
  67    %
  68    %first get the row index of G & its corresponding permutation
  69    %elements.
  70    %
  71    P{1}=[0];
  72    for idx=1:M
  73      P{idx+1}=idx;
  74    end
  75    idx=idx+1;
  76
  77    Ufull=1:M;
  78    r=2;
  79    while r <= R
  80      TMP=nchoosek(Ufull,r);
  81      for idy=1:nchoosek(M,r)
  82        P{idx+idy}=TMP(idy,:);
  83      end
  84      idx=idx+idy;
  85      r=r+1;
  86    end
  87
  88   %
  89   %enter majority logic decoding loop, R+1 order polynomial,
  90   %but we do it here for n-k times, both are equivalent.
  91   %
  92
  93   NCODES=size(VV);
  94   NCODES=NCODES(1);
  95   v_adjust=[];
  96
  97   for row_v=1:1:NCODES
  98    V=VV(row_v,:);
  99    CM=-1*ones(1,Rows);
 100
 101    %
 102    %Now start at bottom row, and get the index set,
 103    %for each until the 2nd most row.
 104    %
 105
 106    %special case, r=0, parity check, so just sum-up.
 107    if ( R == 0 )
 108       wt=__majority_logic_vote(V);
 109       CMM(row_v,:)=wt;
 110       C(row_v,:)=mod(wt*G,2);
 111       continue;
 112    end
 113
 114    order=R;
 115    Gadj=G;
 116    prev_len=length(P{Rows});
 117    for idx=Rows:-1:1
 118      %
 119      %adjust the 'V' received vector, at change of each order.
 120      %
 121      if ( prev_len ~= length(P{idx})  || idx == 1 ) %force for_ idx=1
 122          v_adjust=mod(CM(idx+1:end)*Gadj(idx+1:end,:),2);
 123          Gadj(idx+1:end,:)=0;
 124          V=mod(V+ v_adjust,2); % + = - in GF(2).
 125          order = order - 1;
 126          if ( order == 0 ) %special handling of the all-1's basis vector.
 127            CM(idx)=__majority_logic_vote(V);
 128            break
 129          end
 130       end
 131
 132       prev_len=length(P{idx});
 133       Si=P{idx};% index identifier
 134       Si=sort(Si,'descend');
 135
 136       %generate index elements
 137       B=__binvec(0:(2.^length(Si)-1));
 138       WTS=2.^[Si-1];
 139       %actual index set elements.
 140       S=sum(B.*repmat(WTS,[2^length(Si),1]),2);
 141
 142       %doing the operation set difference U \ S to get SCi
 143       SCi=U;
 144       Si_diff=Si-1;
 145       rmidx=[];
 146       for idy=1:M
 147         if( any( Si_diff==SCi(idy) ) )
 148           rmidx=[rmidx, idy];
 149         end
 150       end
 151       SCi(rmidx)=[];
 152       SCi=sort(SCi,'descend');
 153
 154       %corner case RM(r=m,m) case
 155       if (length(SCi) > 0 )
 156         %generate the set SC,
 157         B=__binvec(0:(2.^length(SCi)-1));
 158         WTS=2.^[SCi];
 159         %actual index set elements.
 160         SC=sum(B.*repmat(WTS,[2^length(SCi),1]),2);
 161       else
 162         SC=[0]; %default, has to be empty set mathematically;
 163       end
 164
 165       %
 166       %next compute the checksums & form the weights.
 167       %
 168       wts=[]; %clear prev history
 169       for id_el = 1:length(SC)
 170         sc_el=SC(id_el);
 171         elems=sc_el + S;
 172         elems=elems+1; %adjust indexing
 173         wt=mod(sum(V(elems)),2);%add elements of V, rx vector.
 174         wts(id_el)= wt;%this is checksum
 175       end
 176
 177       %
 178       %do the majority logic vote.
 179       %
 180       CM(idx)=__majority_logic_vote(wts);
 181   end
 182
 183   CMM(row_v,:)=CM;
 184   C(row_v,:)=mod(CM*G,2);
 185   end
 186   return;
 187 end
 188
 189 %
 190 % utility functions
 191 %
 192 function bvec=__binvec(dec_vec)
 193      maxlen=ceil(log2(max(dec_vec)+1));
 194      x=[]; bvec=zeros(length(dec_vec),maxlen);
 195      for idx=maxlen:-1:1
 196          tmp=mod(dec_vec,2);
 197          bvec(:,idx)=tmp.';
 198          dec_vec=(dec_vec-tmp)./2;
 199      end
 200      return
 201 end
 202
 203 %
 204 % crude majority logic decoding; force the = case to 0 by default.
 205 %
 206 function wt=__majority_logic_vote(wts)
 207       wt=sum(wts)-sum(1-wts);%count no of 1's - no of 0's.
 208       if ( wt ~= 0 )
 209          wt = (wt > 0);
 210          %else
 211          %wt = rand() > 0.5; %break the tie.
 212         %end
 213       end
 214 end
 215
 216 %
 217 % majority logic decoding, tie-break using random.
 218 %
 219 function wt=__majority_logic_vote_random(wts)
 220     wt=(1+sign( sum(wts)-sum(1-wts) ))/2;
 221     if ( wt == 0.5 )
 222         wt = (rand()>0.5);
 223     end
 224 end
 225
 226 % test cases
 227 %G=[1 1 1 1,1 1 1 1;
 228 %   0 1 0 1,0 1 0 1;
 229 %   0 0 1 1,0 0 1 1;
 230 %   0 0 0 0 1 1 1 1];
 231 %m=[1 0 0 1];
 232 %c=mod(m*G,2);
 233 %c(1)=1-c(1); %corrects errors!
 234 %[dc,dm]=reedmullerdec(c,G,1,3)
 235 %pause
 236 %
 237 %G=reedmullergen(1,4);
 238 %m=[1 0 0 0 1];
 239 %c=mod(m*G,2);
 240 %[dc,dm]=reedmullerdec(c,G,1,4)
 241 %pause
 242 %
 243 %G=reedmullergen(3,4);
 244 %m=[ones(1,15)];
 245 %c=mod(m*G,2);
 246 %[dc,dm]=reedmullerdec(c,G,3,4)
 247 %pause
 248 %
 249 %G=reedmullergen(2,3);
 250 %m=[0 0 0 1 1 1 1]
 251 %c=mod(m*G,2)
 252 %[dc,dm]=reedmullerdec(c,G,2,3)
 253 %pause
 254 %
 255 %G=reedmullergen(3,3);
 256 %c1=mod([ones(1,8)]*G,2);
 257 %c2=mod([ones(1,4),zeros(1,4)]*G,2);
 258 %[dC,dM]=reedmullerdec([c2;c2;c1;c2],G,3,3)
 259 %
 260 % %special case of repetition code.
 261 % G=reedmullergen(0,3);
 262 % G
 263 % c1=1*G;
 264 % c2=0*G; C=[c1; c2]
 265 % [dC,dM]=reedmullerdec(C,G,0,3)
 266 %