X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fcommunications-1.1.1%2Freedmullerdec.m;fp=octave_packages%2Fcommunications-1.1.1%2Freedmullerdec.m;h=5feec06419374c7167bbcb2747c0f01d77af1782;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/communications-1.1.1/reedmullerdec.m b/octave_packages/communications-1.1.1/reedmullerdec.m new file mode 100644 index 0000000..5feec06 --- /dev/null +++ b/octave_packages/communications-1.1.1/reedmullerdec.m @@ -0,0 +1,266 @@ +## Copyright (C) 2007, 2011 Muthiah Annamalai +## +## 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 . + +## -*- 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) +%