--- /dev/null
+## Copyright (C) 2007, 2011 Muthiah Annamalai <muthiah.annamalai@uta.edu>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} reedmullerdec (@var{VV},@var{G},@var{R},@var{M})
+##
+## Decode the received code word @var{VV} using the RM-generator matrix @var{G},
+## of order @var{R}, @var{M}, returning the code-word C. We use the standard
+## majority logic vote method due to Irving S. Reed. The received word has to be
+## a matrix of column size equal to to code-word size (i.e @math{2^m}). Each row
+## is treated as a separate received word.
+##
+## The second return value is the message @var{M} got from @var{C}.
+##
+## G is obtained from definition type construction of Reed Muller code,
+## of order @var{R}, length @math{2^M}. Use the function reedmullergen,
+## for the generator matrix for the (@var{R},@var{M}) order RM code.
+##
+## Faster code constructions (also easier) exist, but since
+## finding permutation order of the basis vectors, is important, we
+## stick with the standard definitions. To use decoder
+## function reedmullerdec, you need to use this specific
+## generator function.
+##
+## see: Lin & Costello, Ch.4, "Error Control Coding", 2nd Ed, Pearson.
+##
+## @example
+## @group
+## G=reedmullergen(2,4);
+## M=[rand(1,11)>0.5];
+## C=mod(M*G,2);
+## [dec_C,dec_M]=reedmullerdec(C,G,2,4)
+##
+## @end group
+## @end example
+##
+## @end deftypefn
+## @seealso{reedmullergen,reedmullerenc}
+
+## FIXME: make possible to use different generators, if permutation
+## matrix (i.e polynomial vector elements of rows of G are given
+function [C,CMM]=reedmullerdec(VV,G,R,M)
+ if ( nargin < 4 )
+ print_usage();
+ end
+
+ %
+ % we do a R+1 level majority logic decoding.
+ % at each order of polynomial modifying the code-word.
+ %
+ U=0:M-1; %allowed basis vectors in V2^M.
+ C=-1*ones(size(VV)); %preset the output word.
+ [Rows,Cols]=size(G);%rows shadows 'rows()'
+
+ %
+ %first get the row index of G & its corresponding permutation
+ %elements.
+ %
+ P{1}=[0];
+ for idx=1:M
+ P{idx+1}=idx;
+ end
+ idx=idx+1;
+
+ Ufull=1:M;
+ r=2;
+ while r <= R
+ TMP=nchoosek(Ufull,r);
+ for idy=1:nchoosek(M,r)
+ P{idx+idy}=TMP(idy,:);
+ end
+ idx=idx+idy;
+ r=r+1;
+ end
+
+ %
+ %enter majority logic decoding loop, R+1 order polynomial,
+ %but we do it here for n-k times, both are equivalent.
+ %
+
+ NCODES=size(VV);
+ NCODES=NCODES(1);
+ v_adjust=[];
+
+ for row_v=1:1:NCODES
+ V=VV(row_v,:);
+ CM=-1*ones(1,Rows);
+
+ %
+ %Now start at bottom row, and get the index set,
+ %for each until the 2nd most row.
+ %
+
+ %special case, r=0, parity check, so just sum-up.
+ if ( R == 0 )
+ wt=__majority_logic_vote(V);
+ CMM(row_v,:)=wt;
+ C(row_v,:)=mod(wt*G,2);
+ continue;
+ end
+
+ order=R;
+ Gadj=G;
+ prev_len=length(P{Rows});
+ for idx=Rows:-1:1
+ %
+ %adjust the 'V' received vector, at change of each order.
+ %
+ if ( prev_len ~= length(P{idx}) || idx == 1 ) %force for_ idx=1
+ v_adjust=mod(CM(idx+1:end)*Gadj(idx+1:end,:),2);
+ Gadj(idx+1:end,:)=0;
+ V=mod(V+ v_adjust,2); % + = - in GF(2).
+ order = order - 1;
+ if ( order == 0 ) %special handling of the all-1's basis vector.
+ CM(idx)=__majority_logic_vote(V);
+ break
+ end
+ end
+
+ prev_len=length(P{idx});
+ Si=P{idx};% index identifier
+ Si=sort(Si,'descend');
+
+ %generate index elements
+ B=__binvec(0:(2.^length(Si)-1));
+ WTS=2.^[Si-1];
+ %actual index set elements.
+ S=sum(B.*repmat(WTS,[2^length(Si),1]),2);
+
+ %doing the operation set difference U \ S to get SCi
+ SCi=U;
+ Si_diff=Si-1;
+ rmidx=[];
+ for idy=1:M
+ if( any( Si_diff==SCi(idy) ) )
+ rmidx=[rmidx, idy];
+ end
+ end
+ SCi(rmidx)=[];
+ SCi=sort(SCi,'descend');
+
+ %corner case RM(r=m,m) case
+ if (length(SCi) > 0 )
+ %generate the set SC,
+ B=__binvec(0:(2.^length(SCi)-1));
+ WTS=2.^[SCi];
+ %actual index set elements.
+ SC=sum(B.*repmat(WTS,[2^length(SCi),1]),2);
+ else
+ SC=[0]; %default, has to be empty set mathematically;
+ end
+
+ %
+ %next compute the checksums & form the weights.
+ %
+ wts=[]; %clear prev history
+ for id_el = 1:length(SC)
+ sc_el=SC(id_el);
+ elems=sc_el + S;
+ elems=elems+1; %adjust indexing
+ wt=mod(sum(V(elems)),2);%add elements of V, rx vector.
+ wts(id_el)= wt;%this is checksum
+ end
+
+ %
+ %do the majority logic vote.
+ %
+ CM(idx)=__majority_logic_vote(wts);
+ end
+
+ CMM(row_v,:)=CM;
+ C(row_v,:)=mod(CM*G,2);
+ end
+ return;
+end
+
+%
+% utility functions
+%
+function bvec=__binvec(dec_vec)
+ maxlen=ceil(log2(max(dec_vec)+1));
+ x=[]; bvec=zeros(length(dec_vec),maxlen);
+ for idx=maxlen:-1:1
+ tmp=mod(dec_vec,2);
+ bvec(:,idx)=tmp.';
+ dec_vec=(dec_vec-tmp)./2;
+ end
+ return
+end
+
+%
+% crude majority logic decoding; force the = case to 0 by default.
+%
+function wt=__majority_logic_vote(wts)
+ wt=sum(wts)-sum(1-wts);%count no of 1's - no of 0's.
+ if ( wt ~= 0 )
+ wt = (wt > 0);
+ %else
+ %wt = rand() > 0.5; %break the tie.
+ %end
+ end
+end
+
+%
+% majority logic decoding, tie-break using random.
+%
+function wt=__majority_logic_vote_random(wts)
+ wt=(1+sign( sum(wts)-sum(1-wts) ))/2;
+ if ( wt == 0.5 )
+ wt = (rand()>0.5);
+ end
+end
+
+% test cases
+%G=[1 1 1 1,1 1 1 1;
+% 0 1 0 1,0 1 0 1;
+% 0 0 1 1,0 0 1 1;
+% 0 0 0 0 1 1 1 1];
+%m=[1 0 0 1];
+%c=mod(m*G,2);
+%c(1)=1-c(1); %corrects errors!
+%[dc,dm]=reedmullerdec(c,G,1,3)
+%pause
+%
+%G=reedmullergen(1,4);
+%m=[1 0 0 0 1];
+%c=mod(m*G,2);
+%[dc,dm]=reedmullerdec(c,G,1,4)
+%pause
+%
+%G=reedmullergen(3,4);
+%m=[ones(1,15)];
+%c=mod(m*G,2);
+%[dc,dm]=reedmullerdec(c,G,3,4)
+%pause
+%
+%G=reedmullergen(2,3);
+%m=[0 0 0 1 1 1 1]
+%c=mod(m*G,2)
+%[dc,dm]=reedmullerdec(c,G,2,3)
+%pause
+%
+%G=reedmullergen(3,3);
+%c1=mod([ones(1,8)]*G,2);
+%c2=mod([ones(1,4),zeros(1,4)]*G,2);
+%[dC,dM]=reedmullerdec([c2;c2;c1;c2],G,3,3)
+%
+% %special case of repetition code.
+% G=reedmullergen(0,3);
+% G
+% c1=1*G;
+% c2=0*G; C=[c1; c2]
+% [dC,dM]=reedmullerdec(C,G,0,3)
+%