Coding Theory Programs Built in Octave. Part II.

#
1;

#================================
function Badd = add(x1,x2)
# description: this function takes two binary vectors x1,x2, and returns x1+x2
# mode 2.
# Note: the representation is binary, not in powers.
# usgae: Badd = add(x1,x2)   
      n=rows(x1);
      x3=zeros(n,1);
      for i=1:n
            m=x1(i)+x2(i);
            if (m==1)  (x3(i)=1); endif;
            if (m==2)  (x3(i)=0); endif;
      endfor;
      Badd=x3;
endfunction
#=====================================
function TwoFact=TF(n)
Prod=1;
for i=1:n
      Prod=Prod*2;
endfor;
TF=Prod;
endfunction
#=====================================

function BtoD=BtD(x)
[n,m]=size(x);
N=max(n,m);
s=x(1);
for i=2:N
      s=s+x(i)*TF(i-1);
endfor;
BtD=s;
endfunction

#=====================================
function w=PolyWt(v)

# description: this function returns the weight of the binary polynomial
# v(x) which is represented as a vector v. The vector v contains the
# powers of v(x).
# usage: w=PolyWt(v)

[n,m]=size(v);
if (m>n) (v=v');  endif;
n=rows(v);
w=n;
endfunction

#=====================================
function BtoD=BtD(x)

# Description: This function takes a number x in the binary form
# and returns the decimal form.
# Usage: BtoD=BtD(x)

[n,m]=size(x);
N=max(n,m);
s=x(1);
for i=2:N
      s=s+x(i)*TF(i-1);
endfor;
BtD=s;
endfunction

#=====================================
function DV=div(a,b)

# Description: This function finds a div b
# Usage: DV=div(a,b)
 
x=a/b;
y=round(x);
if (y==x) (DV=x); endif;
if (y<>x) (DV=y-1); endif;
endfunction
#=====================================
function MD=mod(a,b)

# Description: This function returns a mod b
# Usage: MD=mod(a,b)

x=a/b;
y=round(x);
if (y==x) (MD=0);endif;
if (y<>x) (MD=(x-(y-1))*b);endif;
endfunction
#=====================================
function DecToBin=DtoB(x,g)

# Description: Changes any positive integer from the decimal form to the
# required system (from 1 to 16, eg: 2 corresponds to binary). Note: x is
# the number to be changed, and g is the system to change to.
# Usage: DecToBin=DtoB(x,g)

dd=0;
k=[];
y=x;
for i=1:1000
      z=div(y,g);
      e=mod(y,g);
      dd=dd+1;
      k=[k,e];
      y=z;
      if (y==0)
            k=k';
            tt=zeros(dd,1);
            for i=1:dd
                  tt(i)=k(dd-i+1);
            endfor;
            DecToBin=tt;
            return;    
      endif;
endfor;
     
endfunction

#=====================================
function Code=CodeGen(M)

# description: This function takes a generator of a certain binary code and
# returns the code itself. Each row of the matrix M is one generator.
# usage: Code=CodeGen(M)

[m,n]=size(M);
Prod=1;
for k=1:m
      Prod=Prod*2;
endfor;
C=zeros(Prod,n);
DD=zeros(Prod,m);
for i=1:Prod-1
      TEMP=DtoB(i,2)';
      S=TEMP;
      if (columns(TEMP)<m)
            dif=m-columns(TEMP);
            q=zeros(1,dif);
            pp=[q,TEMP];
            S=pp;
      endif;
      DD(i,:)=S;
endfor;
for i=1:Prod-1
      MM=M;
      for j=1:m
            MM(j,:)=MM(j,:)*DD(i,j);
      endfor;
      summ=MM(1,:);
      for j=2:m
            summ=add(summ',MM(j,:)')';
      endfor;
      C(i,:)=summ;
endfor;

Code=C;
endfunction

#=====================================
function weight=wt(v)

# Description: This function takes a certain vector v and returns
# its weight; i.e. the number of ones.
# NOTE: the representation is binary, not in powers.
# usage: weight=wt(v)

[n,m]=size(v);
if (m>n) (v=v');  endif;
n=rows(v);
w=0;
for i=1:n
      if (v(i)==1) (w=w+1); endif;
endfor;
weight=w;
endfunction

#=====================================

function MatVecMul=Bmult(M,v)

# Description: This function takes a certain binary matrix M and a binary
# vector v, and returns M*v mod 2.
# NOTE: the representation is binary, not in powers.
# Usage: MatVecMul=Bmult(M,v)

#[n,m]=size(v);
#if (m>n) (v=v')  endif;
[n,m]=size(M);
w=zeros(n,1);
for i=1:n
      qq=zeros(m,1);
      for j=1:m
            qq(j)=M(i,j)*v(j);
      endfor;
      count=0;
      for k=1:m
            if (qq(k)==1) (count=count+1); endif;
      endfor;
      zz=count/2;
      yy=round(zz);
      if (yy==zz) (w(i)=0); endif;
      if (yy<>zz) (w(i)=1); endif;
endfor;
w=w; #new
MatVecMul=w;
           
endfunction
#=====================================
function [MVadd,flag1] = ad(M,v,check)
# description: This function takes a binary vector v,a binary matrix M, and
# returns M+v mod 2.
# usgae: [MVadd,flag1] = ad(M,v,check) 
      [nn,n]=size(M);flag1=0;M1=M;
      x=zeros(n,1);y=zeros(n,1);
      for i=1:nn
            for j=1:n
                  x(j)=M1(i,j);
            endfor;
            y=add(x,v);
            for k=1:8
                  if (check(k,:)==y')
                        (flag1=1);
                  endif;
            endfor;
            for j=1:n
                  M1(i,j)=y(j);
            endfor;

      endfor;
      MVadd=M1;
endfunction
#================
function BinElemAdd=BEA(x1,x2)

# description: This function takes the two binary vectors x1 and x2, and
# returns x1+x2 mod 2.
# NOTE: the representation is binary, not in powers.
# usage: BinElemAdd=BEA(x1,x2)

x=x1+x2;
BinElemAdd=0;
if (x==1) (BinElemAdd=1) endif;
endfunction
#================
function LinFSR1=lfsr1(x,g,n)

# description: This function returns does LFSR.
# NOTE: The polynomials g and x have binary representation (not in powers).
# This function performs one cycle only.
# usgae: LinFSR1=lfsr1(x,g,n)
[n1,m1]=size(x);
if (m1>n1) (x=x'); endif;
[n1,m1]=size(x);
[n2,m2]=size(g);
if (m2>n2) (g=g'); endif;
[n2,m2]=size(g);
b=x;
temp=b(n2-1);
for j=(n2-1):-1:2
      if (g(j)==1) (b(j)=BEA(b(j-1),temp)); endif;
      if (g(j)<>1) (b(j)=b(j-1));endif;
endfor;
b(1)=temp;
LinFSR1=b';
endfunction
#================
function LinFSR=lfsr(x,g,n)

# description: This function returns does LFSR.
# NOTE: The polynomials g and x have binary representation (not in powers).
# usgae: LinFSR=lfsr(x,g,n) 

[n1,m1]=size(x);
if (m1>n1) (x=x'); endif;
[n1,m1]=size(x);
[n2,m2]=size(g);
if (m2>n2) (g=g'); endif;
[n2,m2]=size(g);
b=zeros(n2-1,1);
count=0;
for i=n2-1:-1:1
      count=count+1;
      b(i)=x(n1-count+1);
endfor;
for i=(n1-(n2-1)):-1:1
      temp=b(n2-1);
      #b(n2-1)=b(n2-2);
      for j=(n2-1):-1:2
            if (g(j)==1) (b(j)=BEA(b(j-1),temp)); endif;
            if (g(j)<>1) (b(j)=b(j-1));endif;
      endfor;
      if (g(1)==1) (b(1)=BEA(x(i),temp));endif;
      if (g(1)<>1) (b(1)=x(i));endif;
endfor;
LinFSR=b';
endfunction
#====================================
function GenMatLfsr=GenMatLFSR(g,n)
# description:
# usage:
[n1,m1]=size(g);
if (m1>n1) (g=g');endif;
m=max(n1,m1)-1
k=n-m
M=zeros(k,m);
u=zeros(m+1,1)';
u(m+1)=1;
u=lfsr(u,g,n);
[n2,m2]=size(u);
if (m2<n2) (u=u');endif;
M(1,:)=u;
count=1;
for i=n-k+1:n-1
      count=count+1;
      u=lfsr1(u,g,n);
      [n2,m2]=size(u);
      if (m2<n2) (u=u');endif;
      M(count,:)=u;
endfor;
GenMatLfsr=[M,eye(k)];
endfunction




#====================================

function [pos,synd,cycDec1]=CycDec1(u,g,n,t)

# description: this function decodes a received vector u, where the code
# used is binary cyclic, g is the generating function and t is the number
# of errors which can be corrected."n" is the dimension of the code.
# NOTE: all the polynomials are represented in the binary form (not powers).
# This funtion returns the vector cycDec which is the binary representation
# (NOT in powers) of the decoded polynomial.
# NOTE: This function uses lfsr function and not long division.
# usage: [pos,synd,cycDec1]=CycDec1(u,g,n,t)

[n1,m1]=size(u);
if (m1>n1) (u=u');endif;
s=lfsr(u,g,n);
n1=max(n1,m1);
for i=1:n
      if (wt(s)<=t)
            q=n-i+1;
            [n2,m2]=size(s);
            if (m2>n2) (s=s');endif;
            m2=max(n2,m2);
            d1=n-n1;
            temp1=zeros(d1,1);
            u=u';
            d2=n-m2;
            temp2=zeros(d2,1)';
            s=s';
            SS=s;
            s=[s,temp2]
            S=s;
            s=zeros(n,1);
            S
            for j=1:n
                  if (S(j)==1);
                        temp=j+q;
                        if (temp>n);
                              s(temp-n)=1;  #modified
                        endif;
                        if (temp<=n) (s(temp)=1);endif;
                  endif;
            endfor;
            s'
            pos=q;
            synd=SS;
            cycDec1=add(u',s')';
            return;
      endif;
s=lfsr1(s,g,n);
endfor;
printf("The error is uncorrectable\n");
endfunction
#=============================

#========================
function [R,LD]=LongDiv(v,g)

# description: this function divides the two polynomials v and g, where v
# and g are elemnts in Z2[x]; i.e. the coefficients are either zeros or
# ones. The polynomial v, as well as g, is represented in its powers;i.e.
# if v(x)=x^6+x^3+X+1, then v is represented as the vector:[6,3,1,0].
# This function takes as an input the two vectors corresponding to the
# two polynomials to be divided, and it returns the vector LD whose elements
# are actually the coefficients v(x)/g(x); and the vector R, whose
# coefficients are actually the powers of the remainder; i.e. R(x)=
# v(x) mod g(x).
# e.g. if v(x)=x^6+x^5+X^3+x^2+1 and g(x)=x^3+x+1,
# then v=[6,5,3,2,0] and g=[3,1,0]. So LD=[3,2,1,0],
# which corresponds to the polynomial LD(x)=x^3+x^2+x+1, and
# R=[2], which corresponds to the polynomial R(x)=x^2.
# NOTE: the function with the highest power has to be written first.
# usage:  [R,LD]=LongDiv(v,g)

v=sort(v);g=sort(g);
[m1,n1]=size(v);
if (n1>m1) (v=v'); endif;
[m2,n2]=size(g);
if (n2>m2) (g=g'); endif;
n=rows(v);m=rows(g);temp=zeros(m,1);w=[];
u=100000*(n+1);
for i=1:u
      if (v==[])
            LD=w';R=0;return;
      endif;
      N=rows(v);

      if (v(N)<g(m))
            LD=w';R=v';return;
      endif;
      q=v(N)-g(m);
      w=[w,q];
      for j=1:m
            temp(j,1)=g(j)+q;
      endfor;
      count=0;
      vNEW=[];
      for j=1:m
            flag=0;
            for k=1:N
                  if (v(k)==temp(j)) (flag=1); endif;
            endfor;
            if (flag==0) (vNEW=[vNEW,temp(j)]); endif;
      endfor;
           
      for k=1:N
            for j=1:m
                  if (v(k)==temp(j))
                        v(k)=-1; count=count+1;
                  endif;
            endfor;
      endfor;
      z=zeros(N-count,1);pos=0;
      for k=1:N
           
            if (v(k)<>-1)
                  pos=pos+1;z(pos,1)=v(k);
            endif;
      endfor;
      z=[z',vNEW];
      z=sort(z);
      z=z';
      v=z;
endfor;
endfunction
#====================================
function BPA=BinPolAdd(p,q)

# description: this function adds the two polynomials p and q, where p
# and q are elemnts in Z2[x]; i.e. the coefficients are either zeros or
# ones. The polynomial p, as well as q, is represented in its powers;i.e.
# if p(x)=x^6+x^3+X+1, then p is represented as the vector:[6,3,1,0].
# This function takes as an input the two vectors corresponding to the
# two polynomials to be added, and it returns the vector BPA whose elements
# are actually the powers of the sum; e.g. if p(x)=x^6+X^3+x+1 and
# q(x)=X^7+x^3+x^2+1, then p=[6,3,1,0] and q=[7,3,2,0]. So BPA=[7,6,2,1],
# which corresponds to the polynomial BPA(x)=x^7+x^6+x^2+x.
# usage:  BPA=BinPolAdd(p,q)

p=sort(p);q=sort(q);
[m1,n1]=size(p);
if (n1>m1) (p=p'); endif;
[m2,n2]=size(q);
if (n2>m2) (q=q'); endif;
n=rows(p);m=rows(q);
count=0;
      vNEW=[];
      for j=1:m
            flag=0;
            for k=1:n
                  if (p(k)==q(j)) (flag=1); endif;
            endfor;
            if (flag==0) (vNEW=[vNEW,q(j)]); endif;
      endfor;    
      for k=1:n
            for j=1:m
                  if (p(k)==q(j))
                        p(k)=-1; count=count+1;
                  endif;
            endfor;
      endfor;
      z=zeros(n-count,1);pos=0;
      for k=1:n
           
            if (p(k)<>-1)
                  pos=pos+1;z(pos,1)=p(k);
            endif;
      endfor;
      z=[z',vNEW];
      z=sort(z);
      BPA=z;
endfunction

#====================================
function BMPM=BinMonPolyMult(p,q,n)

# description: this function multiplies the two polynomials p and q, where p
# and q are elemnts in Z2[x]; i.e. the coefficients are either zeros or
# ones. The polynomial p, as well as q, is represented in its powers;i.e.
# if p(x)=x^6+x^3+X+1, then p is represented as the vector:(6,3,1,0).
# NOTE: q must be a monomial represented in a form of a single vector;
# i.e. if q(x)=x^7, then q=[7], if q(x)=1, then q=[0].
# This function takes as an input the two vectors corresponding to the
# two polynomials to be multiplied, and it returns the vector BPA whose
# elements are actually the coefficients of the product mod (x^n - 1); e.g.
# if n>=10 and
# p(x)=x^6+X^3+x+1 and q(x)=X^4, then p=[6,3,1,0] and q=[4]. So BPA=[10,7,5,4],
# which corresponds to the polynomial BPA(x)=x^10+x^7+x^7+x^4.
# As another example, if n=7 and p(x),q(x) as above, then BPA=[0,3,4,5].
# NOTE: n here is the code lenght.
# usage:  BMPM=BinMonPolMult(p,q,n)

p=sort(p);
[m1,n1]=size(p);
if (n1>m1) (p=p'); endif;
n1=rows(p);
s=q(1);
for i=1:n1
      tt=p(i,1)+s;
      if (tt>=n)
            (p(i,1)=tt-n);   
      endif;
      if (tt<n)
            p(i,1)=tt;
      endif;
endfor;
BMPM=p';
endfunction
#====================================

function gcd=GCD(p,q)

# description: this function finds the greatest common divisor (gcd) of the two # polynomials p and q, where p and q are elemnts in Z2[x];
# i.e. the coefficients are either zeros or
# ones. The polynomial p, as well as q, is represented in its powers;i.e.
# if p(x)=x^6+x^3+X+1, then p is represented as the vector:[6,3,1,0].
# This function takes as an input the two vectors corresponding to the
# powers of the two polynomials which we want to find their greatest common
# divisor (gcd), and it returns the vector gcd whose elements are actually the
# coefficients of the gcd (not the powers).
# e.g. if p(x)= x^8+1, and q(x)=x^4+x^3+x^2+1 , then p=[0,8] and
# q=[0,2,3,4]. So gcd=[0,1], which corresponds to the polynomial
# gcd(x)=1+x.
# NOTE: the function with the highest power has to be written first.
# usage:  gcd=GCD(p,q)

for i=1:100
      [R,LD]=LongDiv(p,q);
      if (R==0)
            gcd=q;
            [m1,n1]=size(gcd);
            if (m1<n1) (gcd=gcd');endif;
            m=rows(gcd);
            n=gcd(m);
            s=zeros(n,1);
            for i=1:m
                  temp=gcd(i);
                  s(temp+1)=1;
            endfor;
            gcd=s';    
            return;
      endif;
      p=q;q=R;
endfor;
endfunction

#====================================

function gcdP=GCDP(p,q)

# description: this function finds the greatest common divisor (gcd) of the two # polynomials p and q, where p and q are elemnts in Z2[x];
# i.e. the coefficients are either zeros or
# ones. The polynomial p, as well as q, is represented in its powers;i.e.
# if p(x)=x^6+x^3+X+1, then p is represented as the vector:[6,3,1,0].
# This function takes as an input the two vectors corresponding to the
# powers of the two polynomials which we want to find their greatest common
# divisor (gcd), and it returns the vector gcd whose elements are actually the
# powers of the gcd (not the coefficients).
# e.g. if p(x)= x^8+1, and q(x)=x^4+x^3+x^2+1 , then p=[0,8] and
# q=[0,2,3,4]. So gcd=[0,1], which corresponds to the polynomial
# gcd(x)=1+x.
# NOTE: the function with the highest power has to be written first.
# usage:  gcdP=GCDP(p,q)

for i=1:100
      [R,LD]=LongDiv(p,q);
      if (R==0)
            gcd=q;
            return;
      endif;
      p=q;q=R;
endfor;
endfunction

#====================================
function BVVM=BinVecVecMult(v,w,n)

# description: This function multiplies the two binary vectors v and w,
# which are represented in powers; i.e. if v(x)=1+x+x^3+x^6, then v(x)
# is represented as v=(0,1,3,6).
# It returns the product represented in powers; i.e. v(x)*w(x). If
# v(x)*w(x)=1+x^3+x^8, then the function returns the vector (0,3,8). 
# usage: BVVM=BinVecVecMult(v,w,n)

[n1,m1]=size(v);
if (n1<m1) (v=v'); endif;
[n1,m1]=size(v);
[n2,m2]=size(w);
if (n2<m2) (w=w'); endif;
[n2,m2]=size(w);
p=w;
z=v(1);
q=[z];
s=BinMonPolyMult(p,q,n);
for i=2:n1
      z=v(i);
      q=[z];
      temp=BinMonPolyMult(p,q,n);
      s=BinPolAdd(temp,s);
endfor;
BVVM=s;
printf("The representation is in powers");printf("\n");
endfunction
#====================================
function MBVVM=MulBinVecVecMult(a,n)

# description: this function multiplies all the binary polynomials which are
# represented in vectors(not in powers) and placed as rows in the matrix a with # each other.
# It returns the product mod (x^n - 1), which is represented in powers.
# usage: MBVVM=MulBinVecVecMult(a,n)

[m,c]=size(a);
v=[];
for j=1:c
      if (a(1,j)==1)
            v=[v,j-1];
      endif;
endfor;
for i=2:m
      w=[];
      for j=1:c
      if (a(i,j)==1)
            w=[w,j-1];
      endif;
endfor;
      v=BinVecVecMult(v,w,n);
endfor;
MBVVM=v;
printf("The representation is in powers");printf("\n");
endfunction
#====================================

function [cycDEC,shift,pos]=cyclicDecoding(u,g,n,t)

# description: this function decodes a received vector u, where the code
# used is binary cyclic, g is the generating function and t is the number
# of errors which can be corrected."n" is the dimension of the code.
# NOTE: all the polynomials are represented in a similar way to what we
# did in previous functions; i.e. g(x) and u(x) are represented in their
# powers.
# This funtion returns the vector cycDec which is the binary representation
# (NOT in powers) of the decoded polynomial.
# usage: [cycDEC,shift,ii]=cyclicDecoding(u,g,n,t)

[n1,m1]=size(u);
if (m1>n1) (u=u');endif;
#uRows=rows(u);
[R,LD]=LongDiv(u,g);
s=R;
n1=max(n1,m1);

for i=1:n
      if (PolyWt(s)<=t)
            q=[n-i+1];
            p=s;
            BMPM=BinMonPolyMult(s,q,n);
            u=BinPolAdd(u,BMPM);
            temp=zeros(n,1);
            u=u';
            uRows=rows(u);
            for j=1:uRows
                  c=u(j);
                  temp(c+1)=1;
            endfor;
            pos=i;
            shift=q;
            cycDEC=temp';
            return;

      endif;
      BMPM=BinMonPolyMult(s,[1],n);
      [R,LD]=LongDiv(BMPM,g);
      s=R;
endfor;
printf("The error is uncorrectable\n");
endfunction
#=============================
function G=CycCodeGenMat(g,n)

# description: this function finds the generator matric G for the code C
# generated by the binary polynomial g(x) represented in the vector form g
# which contains the powers available in g(x);i.e. if g(x)=x^6+x^3+X+1, then p # is represented as the vector: g=(6,3,1,0).
# n here represents the length of the code.
# This function takes as an input the vector g corresponding to the generating
# polynomial, and it returns the generating matrix G, written in
# the standard form; i.e. the last columns of G are Ik, where k is the
# dimension of the code; i.e. k=n-deg(g).
# usage:  G=CycCodeGenMat(g,n)

g=sort(g);
[n1,m1]=size(g);
if (m1>n1) (g=g');endif;
m=max(g)
k=n-m;
B=zeros(k,m);
count=0;
for i=m:n-1
      count=count+1;
      v=[i];
      [R,LD]=LongDiv(v,g);   
      temp=zeros(m,1);
      R=sort(R);
      R=R';
      Rrows=rows(R);
      for j=1:Rrows
            e=R(j);
            temp(e+1)=1;
      endfor;
      temp=temp';
      B(count,:)=temp;
endfor;
G=[B,eye(k)];
endfunction
 

Back to Iyad Abu-Jeib's Homepage