added macros, travis and gitignore

17 files changed, 1698 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..b25c15b
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+*~
diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 0000000..1b45475
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,5 @@
+language: scilab
+
+before_install:
+ - sudo apt-get install scilab
+
diff --git a/macros/alignsignals.sci b/macros/alignsignals.sci
new file mode 100644
index 0000000..2013e3b
--- /dev/null
+++ b/macros/alignsignals.sci
@@ -0,0 +1,166 @@
+function varargout = alignsignals(x,y,varargin)
+//This function aligns the two input signals.
+//
+//Calling Sequence
+//[Xa Ya] = ALIGNSIGNALS(X,Y)
+//[Xa Ya] = ALIGNSIGNALS(X,Y,MAXLAG)
+//[Xa Ya] = ALIGNSIGNALS(X,Y,MAXLAG,1)
+//[Xa Ya D] = ALIGNSIGNALS(...)
+//
+//Description
+//[Xa Ya] = ALIGNSIGNALS(X,Y) aligns the two vectors X and Y by estimating
+//the delay D between the two. If Y is delayed with respect to X, D is
+//   positive , and X is delayed by D samples. If Y is advanced with respect
+//   to X, D is negative, and Y is delayed by -D samples.
+//
+//   [Xa Ya] = ALIGNSIGNALS(X,Y,MAXLAG) considers MAXLAG be the maximum correlation
+//   window size which is used to calculate the estimated delay D between X and Y. 
+//   MAXLAG is an integer-valued scalar. By default, MAXLAG is equal to MAX(LX,LY)-1.
+//   If MAXLAG is empty ([]),then default value is considered. If MAXLAG
+//   is negative, it is replaced by its absolute value.
+//
+//   [Xa Ya] = ALIGNSIGNALS(X,Y,MAXLAG,1) keeps the lengths of Xa
+//   and Ya the same as those of X and Y, respectively. 
+//   Here, 1 implies truncation of the intermediate vectors.
+//   Input argument 4 is 0 implies truncation_off (no truncation).
+//   D is positive implies D zeros are pre-pended to X, and the last D samples of X are truncated.
+//   D is negative implies -D zeros are pre-pended to Y, and the last -D samples
+//   of Y are truncated. That means, when D>=Length(X), all samples of X are lost.
+//   Similarly, when -D>=Length(Y), all samples of Y are lost. 
+//   Avoid assigning a specific value to MAXLAG when using the truncate=1 option, set MAXLAG to [].
+//
+//   [Xa Ya D] = ALIGNSIGNALS(...) returns the estimated delay D.
+//
+//   Examples
+//   X = [0 0 0 1 2 3 ];
+//   Y = [1 2 3 ];
+//   [Xa,Ya] = alignsignals(X,Y,[],1)
+//
+//   See also 
+//   finddelay
+//
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+// Check number of input arguments
+[out_a,inp_a]=argn(0)
+
+if inp_a<=1 | inp_a>4 then
+    error('comm:alignsignals: Invalid number of input arguments')
+end
+
+
+if out_a>3 then
+    error('comm:alignsignals: Invalid number of output arguments')
+end
+
+//Check input arguments
+if (~or(type(x)==[1 5 8]) | ~isvector(x) )
+    error('comm:alignsignals:Input argument 1 should be a vector of numbers');
+end
+
+[row_x,col_x] = size(x);
+len_x = length(x);
+
+
+if (~or(type(y)==[1 5 8]) | ~isvector(y) )
+    error('comm:alignsignals:Input argument 2 should be a vector of numbers');
+end
+
+[row_y,col_y] = size(y);
+len_y = length(y);
+
+//Check for MaxLag
+if inp_a==3 then
+    maxlag = varargin(1)
+else
+    maxlag = max(len_x,len_y)-1; //Default value
+end
+
+if  ~isempty(maxlag) then
+    if ( ~or(type(maxlag)==[1 5 8]) | ~isreal(maxlag) | length(maxlag)~=1 | ceil(maxlag)~=maxlag),
+        error('comm:alignsignals:Input argument 3 should be a scalar integer');   
+    elseif (( isnan(maxlag)) | isinf(maxlag)),
+        error('comm:alignsignals:Input argument 3 can not be Inf or NAN');
+    end
+else
+    maxlag = maxlag_default;
+end
+
+maxlag = double(abs(maxlag));
+
+//Check for truncate
+trunc_on=0;
+if inp_a==4
+    if (varargin(2))
+        trunc_on=1;
+    end
+end    
+
+
+// Estimate delay between X and Y.
+if inp_a==2
+    d = finddelay(x,y);
+else
+    d = finddelay(x,y,maxlag);
+end
+
+
+if d == 0   // X and Y are already aligned.
+    varargout(1) = x;
+    varargout(2) = y;    
+elseif d > 0    // Y is estimated to be delayed with respect to X.
+    if row_x>1 // X is entered as a column vector.
+        if trunc_on==0
+            varargout(1) = [zeros(d,1) ; x]; 
+        else
+            if d>=row_x
+                warning('comm:alignsignals:firstInputTruncated');
+                varargout(1) = zeros(row_x,1);
+            else
+                varargout(1) = [zeros(d,1) ; x(1:len_x-d)];
+            end                
+        end
+    else    // X is entered as a row vector.
+        if trunc_on==0
+            varargout(1) = [zeros(1,d) x];
+        else
+            if d>=col_x
+                warning('comm:alignsignals:firstInputTruncated');
+                varargout(1) = zeros(1,col_x);
+            else
+                varargout(1) = [zeros(1,d) x(1:len_x-d)];
+            end                
+        end    
+    end            
+    varargout(2) = y;
+else    // X is estimated to be delayed with respect to Y.
+    varargout(1) = x;
+    if row_y>1 // Y is entered as a column vector.
+        if trunc_on==0
+            varargout(2) = [zeros(-d,1) ; y];
+        else
+            if (-d)>=row_y
+                warning('comm:alignsignals:secondInputTruncated');
+                varargout(2) = zeros(row_y,1);
+            else
+                varargout(2) = [zeros(-d,1) ; y(1:len_y-(-d))];
+            end                
+        end
+    else    // Y is entered as a row vector.
+        if trunc_on==0
+            varargout(2)= [zeros(1,-d) y];
+        else
+            if (-d)>=col_y
+                warning('comm:alignsignals:secondInputTruncated');
+                varargout(2) = zeros(1,col_y);
+            else
+                varargout(2) = [zeros(1,-d) y(1:len_y-(-d))];
+            end   
+        end    
+    end    
+end    
+    
+varargout(3) = d;
+
+endfunction
diff --git a/macros/arithdeco.sci b/macros/arithdeco.sci
new file mode 100644
index 0000000..95a06b9
--- /dev/null
+++ b/macros/arithdeco.sci
@@ -0,0 +1,170 @@
+function [seq] = arithdeco(code, count, len)
+//   This function decodes the given code using arithmetic coding    
+    
+//   Calling sequence
+//   SEQ = ARITHDECO(CODE, COUNT, LEN) 
+//
+//   Description
+//   SEQ = ARITHDECO(CODE, COUNT, LEN) decodes the given received seq (CODE) to message using arithmetic coding.
+//   COUNT is vector whihc gives information about the source statistics (i.e. frequency of each symbol in the source alphabet)
+//   CODE is the binary arithmetic code
+    
+//   Source Alphabet is assumed to be {1,2,....N} where N is a positive integer
+//   Therefore, sequence should be finite and positive 
+//   Length of the COUNT should match the length of the source alphabet
+
+//   Examples 
+//   counts = [40 1 9];
+//   len = 11;
+//   seq = [1 3 2 1 1 1 3 3 1 1 2 ]
+//   code = arithenco(seq,counts);
+//   disp(code)
+//   dseq=arithdeco(code,counts,len)
+//   disp(dseq)
+//   disp(seq)
+    
+//   Bibliography 
+//   Sayood, K., Introduction to Data Compression, Morgan Kaufmann, 2000, Chapter 4, Section 4.4.3.
+    
+//   See also 
+//   arithenco
+    
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+
+//*************************************************************************************************************************************//
+    
+    //Input argument check
+    
+    [outa,inpa]=argn(0);
+    if(~inpa==3)
+    error("comm:arithenco:Wrong number of Input Arguments");
+    end
+    
+    [row_code,col_code]=size(code);
+    [row_sta,col_sta]=size(count);
+    
+    // Check to make sure that sequence is 1D
+    if(~(row_code==1|col_code==1))
+        error("comm:arithenco: Invalid dimensions: Input Arithmetic Encoded Sequence should be 1D ");
+    end
+    
+    // Check for source statistics matrix 
+    if(~(row_sta==1|col_sta==1))
+        error("comm:arithenco: Invalid dimensions: Argument 2 should be 1D ");
+    end
+    
+    if(~isreal(code) | or(code<0) )
+        error("comm:arithenco: Input sequence should be finite positive integer");
+    end
+    
+    if(~isreal(count) | or(count<0) )
+        error("comm:arithenco: Source statistics should be finite positive integer");
+    end
+    
+    if(~isreal(len) | or(len<=0) | ~isscalar(len))
+        error("comm:arithenco: length should be finite positive integer and scalar");
+    end
+    
+   //Check the incoming orientation and adjust if necessary
+    if (row_code > 1),
+        code = code.';
+    end
+    
+    if (row_sta > 1),
+        count = count.';
+    end
+    
+    // Check if the given code is binary
+    for i=1:length(code)
+    if (~(code(i)==1 | code(i)==0 ))
+        error("comm:arithenco:Input Arithmetic Encoded Sequence is not binary");
+    end
+    end
+    
+    [row_s,col_s]=size(code);
+    [row_c,col_c]=size(count);
+    
+    //Calculate the cumulative count
+    cum_count=[0,cumsum(count)];
+    total_count=cum_count(length(cum_count));
+    
+    //Initialization
+    m=ceil(log2(total_count)) + 2;
+    low=zeros(1,m);
+    up=ones(1,m);
+    dec_low=0;
+    dec_up=2^m-1;
+    seq= zeros(1,len);
+    seq_index=1;
+    k=m;
+    tag=code(1:m);
+    dec_tag=0;
+    value=0;
+    for i=1:length(tag)
+        dec_tag=dec_tag+tag(i)*2^(length(tag)-i);
+    end
+    
+    //loop till you decode entire seq
+    while (seq_index <= len)
+
+        // Compute value
+        value =floor( ((dec_tag-dec_low+1)*total_count-1)/(dec_up-dec_low+1) );
+        
+        //Decode the symbol and update it
+        c=find(cum_count <= value)
+        ptr=c(length(c))
+        
+        seq(seq_index)=ptr;
+        seq_index=seq_index+1;
+        
+        //Compute lower and upper bounds
+        dec_low_new = dec_low + floor( (dec_up-dec_low+1)*cum_count(ptr)/total_count );
+        dec_up = dec_low + floor( (dec_up-dec_low+1)*cum_count(ptr+1)/total_count )-1;
+        dec_low = dec_low_new;
+        
+        for i=1:m
+        low(i)=strtod(part(dec2bin(dec_low,m),i))
+        end
+        for i=1:m
+        up(i)=strtod(part(dec2bin(dec_up,m),i))
+        end
+
+        //Loop while E1, E2 or E3 condition
+        while(isequal(low(1),up(1))) | (isequal(low(2),1) & isequal(up(2),0))
+        
+        if (k==length(code)) then
+            break;
+        end
+        
+        k=k+1;
+            if isequal(low(1),up(1)) then //E1 or E2 holds
+                
+                //Left shift
+                low=[low(2:m) 0];
+                up=[up(2:m) 1];
+                tag=[tag(2:m) code(k)];
+                
+            elseif (isequal(low(2),1) & isequal(up(2),0)) then //for E3 
+                
+                //left shift
+                low=[low(2:m),0];
+                up=[up(2:m),1];
+                tag=[tag(2:m) code(k)];
+                
+                low(1)=bitxor(low(1),1);
+                up(1)=bitxor(up(1),1);
+                tag(1)=bitxor(tag(1),1);
+                
+            end
+        end
+        dec_low=0;dec_up=0;dec_tag=0;
+        for i=1:length(low)
+        dec_low=dec_low+low(i)*2^(length(low)-i);
+        dec_up=dec_up+up(i)*2^(length(up)-i);
+        dec_tag=dec_tag+tag(i)*2^(length(tag)-i);
+        end
+    end
+    
+endfunction
diff --git a/macros/arithenco.sci b/macros/arithenco.sci
new file mode 100644
index 0000000..a463afa
--- /dev/null
+++ b/macros/arithenco.sci
@@ -0,0 +1,166 @@
+function [code] = arithenco(seq, count)
+//   This function encodes the given sequence using aritmetic coding
+
+//   Calling sequence
+//   CODE = ARITHENCO(SEQ, COUNT)
+
+//   Description
+//   CODE = ARITHENCO(SEQ, COUNT) encodes the given sequence (SEQ) using arithmetic coding.
+//   COUNT is vector whihc gives information about the source statistics (i.e. frequency of each symbol in the source alphabet)
+//   CODE is the binary arithmetic code
+//   Source Alphabet is assumed to be {1,2,....N} where N is a positive integer
+//   Therefore, sequence should be finite and positive 
+//   Length of the COUNT should match the length of the source alphabet
+
+//   Examples    
+//   counts = [40 1 9];
+//   len = 4;
+//   seq = [1 3 2 1]
+//   code = arithenco(seq,counts);
+//   disp(code)
+    
+//   Bibliography
+//   Sayood, K., Introduction to Data Compression, Morgan Kaufmann, 2000, Chapter 4, Section 4.4.3.
+
+//   See also 
+//   arithdeco
+    
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+
+//*************************************************************************************************************************************//
+    
+    //Input argument check
+    [outa,inpa]=argn(0);
+    if(~inpa==2)
+    error("comm:arithenco:Wrong number of Input Arguments");
+    end
+    
+    [row_seq,col_seq]=size(seq);
+    [row_sta,col_sta]=size(count);
+    
+    // Check to make sure that sequence is 1D
+    if(~(row_seq==1|col_seq==1))
+        error("comm:arithenco: Invalid dimensions: Input Sequence should be 1D ");
+    end
+    
+    // Check for source statistics matrix 
+    if(~(row_sta==1|col_sta==1))
+        error("comm:arithenco: Invalid dimensions: Argument 2 should be 1D ");
+    end
+    
+    if(~isreal(seq) | or(seq<0) )
+        error("comm:arithenco: Input sequence should be finite positive integer");
+    end
+    
+    if(~isreal(count) | or(count<0) )
+        error("comm:arithenco: Source statistics should be finite positive integer");
+    end
+    
+    //Check if number of elements in source alphabet is equal to dimensions of source statistics matrix
+    if max(seq) > length(count)
+        error("comm:arithenco:Source alphabet size does not match with source statistics size");
+    end
+    
+    //Check the incoming orientation and adjust if necessary
+    
+    if (row_seq > 1),
+        seq = seq.';
+    end
+    
+    if (row_sta > 1),
+        count = count.';
+    end
+    
+    [row_s,col_s]=size(seq);
+    [row_c,col_c]=size(count);
+    
+    //Calculate the cumulative count
+    cum_count=[0,cumsum(count)];
+    scale3=0;
+    total_count=cum_count(length(cum_count));
+    
+    //Initialization
+    m=ceil(log2(total_count)) + 2;
+    low=zeros(1,m);
+    up=ones(1,m);
+    dec_low=0;
+    dec_up=2^m-1;
+    code_len = length(seq) * ( ceil(log2(length(count))) + 2 ) + m;
+    code = zeros(1, code_len);
+    code_index = 1;
+    
+    
+    // For each bit in the seq
+    
+    for k=1:length(seq)
+        symbol = seq(k);
+        // Compute the lower and upper bounds
+        dec_low_new = dec_low + floor( (dec_up-dec_low+1)*cum_count(symbol+1-1)/total_count );
+        dec_up = dec_low + floor( (dec_up-dec_low+1)*cum_count(symbol+1)/total_count )-1;
+        dec_low = dec_low_new;
+        
+        for i=1:m
+        low(i)=strtod(part(dec2bin(dec_low,m),i))
+        end
+        for i=1:m
+        up(i)=strtod(part(dec2bin(dec_up,m),i))
+        end
+
+        //Loop while E1, E2 or E3 condition
+        while(isequal(low(1),up(1))) | (isequal(low(2),1) & isequal(up(2),0))
+            // Check for E1 or E2 condition
+            if isequal(low(1),up(1)) then //E1 or E2 holds
+                // Transmit MSB
+                b=low(1);
+                code(code_index) = b;
+                code_index = code_index + 1;
+                
+                //Left shift
+                low=[low(2:m) 0];
+                up=[up(2:m) 1];
+                while (scale3>0)
+                code(code_index) = bitxor(b,1);
+                code_index = code_index + 1;
+                scale3 = scale3-1;
+                end
+            end
+            if (isequal(low(2),1) & isequal(up(2),0)) then //for E3 
+                //left shift
+                low=[low(2:m),0];
+                up=[up(2:m),1];
+                //disp(up,low,"up,low");
+
+                low(1)=bitxor(low(1),1);
+                up(1)=bitxor(up(1),1);
+
+                scale3=scale3+1;
+            end
+        end
+        dec_low=0;dec_up=0;
+        for i=1:length(low)
+        dec_low=dec_low+low(i)*2^(length(low)-i);
+        dec_up=dec_up+up(i)*2^(length(up)-i);
+        end
+    end
+    if scale3==0 then
+    code(code_index:code_index + m - 1) = low;
+    code_index = code_index + m;
+    end
+    if ~scale3==0 then
+    b=low(1);
+    code(code_index)=b;
+    code_index = code_index + 1;
+    
+    while (scale3>0)
+        code(code_index) = bitxor(b,1);
+        code_index = code_index + 1;
+        scale3 = scale3-1;
+    end
+    
+    code(code_index:code_index+m-2) = low(2:m);
+    code_index = code_index + m - 1;
+    end
+code = code(1:code_index-1);
+endfunction
diff --git a/macros/build_help.sce b/macros/build_help.sce
new file mode 100644
index 0000000..7e477d7
--- /dev/null
+++ b/macros/build_help.sce
@@ -0,0 +1,5 @@
+help_lang_dir = get_absolute_file_path("build_help.sce");
+TOOLBOX_TITLE="Communications Toolbox"
+tbx_build_help(TOOLBOX_TITLE, help_lang_dir);
+ok=add_help_chapter("Demo",get_absolute_file_path("build_help.sce"));
+clear help_lang_dir;
diff --git a/macros/finddelay.sci b/macros/finddelay.sci
new file mode 100644
index 0000000..eb8e96a
--- /dev/null
+++ b/macros/finddelay.sci
@@ -0,0 +1,187 @@
+function d = finddelay(x,y,varargin)
+//   This function returns the estimated delay between two input signals using crosscorrelation.
+//   If signals are periodic, delay with least absolute value is returned.
+
+//   Calling Sequence
+//   D = FINDDELAY(X,Y)
+//   D = FINDDELAY(...,MAXLAG)
+
+//   Description
+//   D = FINDDELAY(X,Y), returns estimated Delay D between X
+//   and Y. D is positive implies Y is delayed with respect to X and vice versa. 
+//   If X, Y are matrices, then D is a row vector corresponding to delay between columns of X and Y
+
+//   D = FINDDELAY(...,MAXLAG), uses MAXLAG as the maximum correlation
+//   window size used to find the estimated delay(s) between X and Y:
+//  
+//  > If MAXLAG is an integer-valued scalar, and X and Y are row or column
+//      vectors or matrices, the vector D of estimated delays is found by
+//      cross-correlating (the columns of) X and Y over a range of lags
+//      -MAXLAG:MAXLAG. 
+//  > If MAXLAG is an integer-valued row or column vector, and one input is vector 
+//      and another be matirx (let X is a row or column vector ,
+//      and Y is a matrix) then the vector D of estimated delays is found by
+//      cross-correlating X and column J of Y over a range of lags
+//      -MAXLAG(J):MAXLAG(J), for J=1:Number of columns of Y.
+//  > If MAXLAG is an integer-valued row or column vector, and X and Y are 
+//      both matrices. then vector D of estimated delays is found by
+//      cross-correlating corresponding columns of X and Y over a range of lags
+//      -MAXLAG(J):MAXLAG(J).
+//
+//   By default, MAXLAG is equal to MAX(LX,LY)-1 for vectors,
+
+//   Examples
+//   X = [ 0 0 1 2 3 ];
+//   Y = [ 0 0 0 1 2 3]; 
+//   D = finddelay(X,Y,2)
+//   disp(D)
+//   X = [ 0 1 0 0 ; 1 0 2 1 ;0 0 0 2 ];
+//   Y = [ 0 0 1 0 ;1 0 0 2 ; 0 0 0 0 ];
+//   D = finddelay(X,Y)
+//   disp(D)
+
+//   See also 
+//   alignsignals
+
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+
+//*************************************************************************************************************************************//
+
+
+// Check number of input arguments
+[out_a,inp_a]=argn(0)
+
+if inp_a<=1 | inp_a>3 then
+    error('comm:finddelay: Invalid number of input arguments')
+end
+
+if out_a>1 then
+    error('comm:finddelay: Invalid number of output arguments')
+end
+
+//Error Checking of input arguments 
+if (~or(type(x)==[1 5 8]) | (isempty(x) ) | (ndims(x)>2) | ~or(type(y)==[1 5 8]) | (isempty(y) ) | (ndims(y)>2))
+    error('comm:finddelay:Input arguments must be numeric');
+end
+
+if isvector(x) 
+    x = x';
+end
+
+if isvector(y) 
+    y = y';
+end
+  
+[row_x,col_x] = size(x);
+[row_y,col_y] = size(y);
+
+x = double(x);
+y = double(y);
+
+
+// Check if matrices are of compatible
+if ~isvector(x) & ~isvector(y)
+    if col_x~=col_y
+        error('comm:finddelay:When both inputs are matrices, they must have the same number of columns.')
+    end
+end
+
+// Check for maxlag
+if inp_a==3  
+    
+    if ( ndims(varargin(1))>2 | ~isreal(varargin(1)) | isempty(varargin(1)) | or(isnan(varargin(1))) | or(isinf(varargin(1))) | varargin(1) ~= ceil(varargin(1))),
+        error('comm:finddelay:Input argument 3 should be a finite integer vector')
+    end
+
+    if ( (isvector(x)) & (isvector(y)) & (length(varargin(1))>1) )
+        error('comm:finddelay: If x and y are both vectors, maxlag should be a scalar')
+    elseif ( (isvector(y)) & (length(varargin(1))>1) & (length(varargin(1))~=col_x) ),
+        error('comm:finddelay: If maxlag is a row/column vector, it should be of same length as the number of columns of X or Y');
+    elseif ( (isvector(x)) & (length(varargin(1))>1) & (length(varargin(1))~=col_y) ),
+        error('comm:finddelay: If maxlag is a row/column vector, it should be of same length as the number of columns of X or Y');
+    elseif ( (length(varargin(1))>1) & (length(varargin(1))~=col_x) & (length(varargin(1))~=col_y) ),
+        error('comm:finddelay: If X and Y are matrices, MAXLAG should be the same length as the number of columns of X and Y.');
+    else
+        if isempty(varargin(1))
+            maxlag = max(row_x,row_y)-1; //default value
+        else
+            maxlag = double(abs(varargin(1)));
+        end
+    end
+else
+    maxlag = max(row_x,row_y)-1;
+end
+
+
+max_col=max(col_x,col_y);
+
+if (length(maxlag)==1)
+    maxlag = repmat(maxlag,1,max_col);
+end
+
+if col_x<max_col
+    x = repmat(x,1,max_col);
+elseif col_y<max_col
+    y = repmat(y,1,max_col);
+end    
+    
+
+// Initialise cross-correlation matrix .
+maxlag_max = max(maxlag);
+c_normalized = zeros(2*maxlag_max+1,max_col);
+index_max = zeros(1,max_col);
+max_c = zeros(1,max_col);
+
+
+// Compute cross-correlation matrix:
+sq_x = abs(x).^2
+sq_y = abs(y).^2
+cxx0 = sum(sq_x,"r");
+cyy0 = sum(sq_y,"r");
+
+for i = 1:max_col
+    if ( (cxx0(i)==0) | (cyy0(i)==0) )
+        c_normalized(:,i) = zeros(2*maxlag_max+1,1);
+    else
+        c_normalized(maxlag_max-maxlag(i)+1:maxlag_max-maxlag(i)+2*maxlag(i)+1,i) ...
+            = abs(xcorr(x(:,i),y(:,i),maxlag(i)))/sqrt(cxx0(i)*cyy0(i));
+    end
+end
+
+// Find lowest positive or zero indices of lags (negative delays) giving the
+// largest absolute values of normalized cross-correlations. 
+
+[max_pos,index_max_pos] = max(c_normalized(maxlag_max+1:$,:),"r");    
+// Find lowest negative indices of lags (positive delays) giving the largest
+// absolute values of normalized cross-correlations. 
+A=c_normalized(1:maxlag_max,:)
+[max_neg,index_max_neg] = max(A($:-1:1,:),"r");
+
+if isempty(max_neg)
+    index_max = maxlag_max + index_max_pos;
+else
+    for i=1:max_col
+        if max_pos(i)>max_neg(i)
+            index_max(i) = maxlag_max + index_max_pos(i);
+            max_c(i) = max_pos(i);
+        elseif max_pos(i)<max_neg(i)
+            index_max(i) = maxlag_max + 1 - index_max_neg(i);
+            max_c(i) = max_neg(i);
+        elseif max_pos(i)==max_neg(i)
+            if index_max_pos(i)<=index_max_neg(i)
+                index_max(i) = maxlag_max + index_max_pos(i);
+                max_c(i) = max_pos(i);
+            else
+                index_max(i) = maxlag_max + 1 - index_max_neg(i);
+                max_c(i) = max_neg(i);
+            end 
+        end   
+    end
+end
+
+// Subtract delays.
+d =  (maxlag_max + 1) - index_max;
+
+endfunction
diff --git a/macros/gfcosets.sci b/macros/gfcosets.sci
new file mode 100644
index 0000000..ed69f86
--- /dev/null
+++ b/macros/gfcosets.sci
@@ -0,0 +1,93 @@
+function gfcs = gfcosets(m, p)
+//   This function produces cyclotomic cosets for a Galois field GF(P)
+//
+//   Calling Sequence
+//   GFCS = GFCOSETS(M)
+//   GFCS = GFCOSETS(M, P)
+// 
+//   Description
+//   GFCS = GFCOSETS(M) produces cyclotomic cosets mod(2^M - 1). Each row of the
+//   output GFCS contains one cyclotomic coset.
+//
+//   GFCS = GFCOSETS(M, P) produces cyclotomic cosets mod(P^M - 1), where 
+//   P is a prime number.
+//       
+//   Because the length of the cosets varies in the complete set, %nan is used to
+//   fill out the extra space in order to make all variables have the same
+//   length in the output matrix GFCS.
+
+
+//   Examples
+//   c = gfcosets(2,3)
+//   disp(c)
+
+//   See also
+//   gfminpol, gfprimdf,  gfroots
+
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+
+//*************************************************************************************************************************************//
+
+//Input argument check
+[out_a,inp_a]=argn(0)
+
+// Error Checking
+if inp_a < 2
+    p = 2;
+elseif ( isempty(p) | ~isscalar(p) | abs(p)~=p | floor(p)~=p | length(factor(p))~=1 | p==1)
+    error('comm:gfcosets: P should be a positive prime number');
+end
+
+if ( isempty(m) | ~isscalar(m) | ~isreal(m) | floor(m)~=m | m<1 )
+    error('comm:gfcosets: M should be a positive integer');
+end
+
+//Initialization
+if ( ( p == 2 ) & ( m == 1 ) )
+    i = [];
+else
+    i = 1;
+end
+
+n = p^m - 1;
+
+gfcs = [];
+
+mk = ones(1, n - 1);
+
+while ~isempty(i)
+
+   mk(i) = 0;
+   s = i;
+   j = s;
+   pk = modulo(p*s, n);
+   //compute cyclotomic coset for s=i
+   while (pk > s)
+          mk(pk) = 0;
+          j = [j pk];
+          pk = modulo(pk * p, n);
+   end;
+
+   // append the coset to gfcs
+   [row_cs, col_cs] = size(gfcs);
+   [row_j, col_j ]  = size(j);
+   if (col_cs == col_j) | (row_cs == 0)
+          gfcs = [gfcs; j];
+   elseif (col_cs > col_j)
+          gfcs = [gfcs; [j, ones(1, col_cs - col_j) * %nan]];
+   else
+
+          gfcs = [[gfcs, ones(row_cs, col_j - col_cs) * %nan]; j];
+   end;
+   i = find(mk == 1,1);        // find the index of next number.
+
+end;
+
+// adding "0" to the first coset
+[row_cs, col_cs] = size(gfcs);
+gfcs = [[0, ones(1, col_cs - 1) * %nan]; gfcs];
+
+endfunction
+
diff --git a/macros/gflineq.sci b/macros/gflineq.sci
new file mode 100644
index 0000000..3b34529
--- /dev/null
+++ b/macros/gflineq.sci
@@ -0,0 +1,170 @@
+function [x, sflag] = gflineq(a, b, p)
+//   This function finds a solution for linear equation Ax = b over a prime Galois field.
+//
+//   Calling Sequence
+//   [X, SFLAG] = GFLINEQ(A, B) 
+//   [X, SFLAG]= GFLINEQ(A, B, P)
+//
+//   Description
+//   [X, SFLAG] = GFLINEQ(A, B) returns a particular solution (X) of AX=B in GF(2).
+//   If the equation has no solution, then X is empty and SFLAG = 0 else SFLAG = 1. 
+//
+//   [X, SFLAG]= GFLINEQ(A, B, P) returns a particular solution of the linear
+//   equation A X = B in GF(P) and SFLAG=1.
+//   If the equation has no solution, then X is empty and SFLAG = 0.
+//
+
+//   Examples
+//   A=[1 0 1; 1 1 0; 1 1 1]
+//   p=3
+//   [x,vld] = gflineq(A,[1;0;1],p)
+//   disp(A,'A=')
+//   disp(x,'x=');
+//   if(vld)
+//     disp('Linear equation has solution x')
+//   else
+//     disp('Linear equation has no solution and x is empty')
+//   end
+//   disp( pmodulo(A*x,p),'B =')
+
+//   See also
+//   gfadd, gfconv, gfdiv, gfrank, gfroots
+
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+
+//*************************************************************************************************************************************//
+ 
+// Check number of input arguments
+[out_a,inp_a]=argn(0)        
+
+if inp_a >3 | out_a> 2 | inp_a <2 then
+    error('comm:gflineq: Invalid number of arguments')
+end
+
+// Error checking .
+if inp_a < 3
+    p = 2;
+elseif ( isempty(p) | length(p)~=1 | abs(p)~=p | ceil(p)~=p | length(factor(p))~=1  )
+    error('comm:gflineq:Input argument 3 must be a positive prime integer.');
+end;
+
+[row_a, col_a] = size(a);
+[row_b, col_b] = size(b);
+
+// Error checking - A & B.
+if ( isempty(a) | ndims(a) > 2 )
+    error('comm:gflineq:Input argument 1 must be a two-dimensional matrix.');
+end
+if ( isempty(b) | ndims(b) > 2 | col_b > 1 )
+    error('comm:gflineq:Invalid dimensions of input argument 2 .');
+end
+if ( row_a ~= row_b )
+    error('comm:gflineq:Dimensions of A and B are not compatible');
+end
+if (( or( abs(a)~=a | floor(a)~=a | a>=p )) | ( or( abs(b)~=b | floor(b)~=b | b>=p )) )
+    error('comm:gflineq:Elements of input matrices should be integers between 0 and P-1.');
+end
+    
+// Solution is found by using row reduction (Reducing it to echelon form) 
+ab = [a b]; // Composite matrix
+[row_ab, col_ab] = size(ab);
+
+row_i = 1;
+col_i = 1;
+row = [];
+col = [];
+
+while (row_i <= row_ab) & (col_i < col_ab)
+
+    // Search for a non zero element in current column
+    while (ab(row_i,col_i) == 0) & (col_i < col_ab)
+
+        idx = find( ab(row_i:row_ab, col_i) ~= 0 );
+
+        if isempty(idx)
+            col_i = col_i + 1; // No non zero element
+        else
+            // Swap the current row with a row containing a non zero element 
+            // (preferably with the row with value 1).
+            idx = [ find(ab(row_i:row_ab, col_i) == 1) idx ];
+            idx = idx(1);
+            temp_row = ab(row_i,:)
+            ab(row_i,:) = ab(row_i+idx-1,:)
+            ab(row_i+idx-1,:) = temp_row
+
+        end
+    end
+    
+    if ( ( ab(row_i,col_i) ~= 0 ) & ( col_i < col_ab ) )
+        
+        // Set major element to 1.
+        if (ab(row_i,col_i) ~= 1)
+           ab(row_i,:) = pmodulo( field_inv( ab(row_i,col_i),p ) * ab(row_i,:), p );
+        end
+
+        // The current element is a major element.
+        row = [row row_i];
+        col = [col col_i];
+        
+        // Find the other elements in the column that must be cleared,
+        idx = find(ab(:,col_i)~=0);
+
+        for i = idx
+            if i ~= row_i
+                ab(i,:) = pmodulo( ab(i,:) + ab(row_i,:) * (p - ab(i,col_i)), p );
+            end
+        end
+
+        col_i = col_i + 1;
+
+    end
+
+
+    row_i = row_i + 1;
+
+
+end
+
+if ( rank(ab) > rank( ab(:,1:col_a) ) )
+    disp('comm:gflineq:Solution does not exist');
+    x = [];
+    sflag = 0;
+else
+    x = zeros(col_a, 1);
+    x(col,1) = ab(row,col_ab);
+    sflag = 1;
+end
+
+endfunction
+
+
+function [x] = field_inv(a,n)
+    t = 0;
+    newt = 1;
+    r = n;
+    newr = a;
+    
+    while newr ~= 0
+        quotient = floor(r / newr);
+        
+        temp = t;
+        t = newt;
+        newt = temp -quotient*newt;
+        
+        temp = r;
+        r = newr;  
+        newr = temp - quotient*newr;      
+    end
+    
+    if r>1
+    [x c] = find( pmodulo( (1:(p-1)).' * (1:(p-1)) , p ) == 1 );
+    end
+    
+    if t<0
+        t = t + n;
+     end
+     x = t;
+     
+endfunction
diff --git a/macros/gfrepcov.sci b/macros/gfrepcov.sci
new file mode 100644
index 0000000..6e16717
--- /dev/null
+++ b/macros/gfrepcov.sci
@@ -0,0 +1,60 @@
+function q = gfrepcov(p)
+//   This function represents a binary polynomial in standard ascending order format.
+
+//   Calling Sequence
+//   Q = GFREPCOV(P)
+
+//   Description
+//   Q = GFREPCOV(P) converts vector (P) to standard ascending
+//   order format vector (Q), which is a vector that lists the coefficients in 
+//   order of ascending exponents,  if P represents a binary polynomial 
+//   as a vector of exponents with non-zero coefficients.
+
+//   Examples
+//   The matrix below represents the binary polynomial $1 + s + s^2 + s^4$
+//   Implies output vector should be [1 1 1 0 1]
+//   A=[0 1 2 4 ]
+//   B=gfrepcov(A)
+//   disp(B)
+//   Also try A=[1 2 3 4 4] which is incorrect way of representing binary polynomial
+
+//   See also
+//   gfpretty
+
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+//*************************************************************************************************************************************//
+
+//Input arguments 
+[out_a,inp_a]=argn(0)
+
+[row_p, col_p] = size(p);
+
+// Error checking 
+if ( isempty(p) | ndims(p)>2 | row_p > 1 )
+    error('comm:gfrepcov: P should be a vector');
+end
+for j=1:col_p
+if (floor(p(1,j))~=p(1,j) | abs(p(1,j))~=p(1,j)  )
+    error('comm:gfrepcov: Elements of the vector should be non negative integers');
+end
+end
+
+// Check if the given vector is in ascending order format, if not convert //
+if max(p) > 1
+
+// Check if P has any repetative elements.
+    
+    if (length(unique(p))~=length(p))
+        error('comm:gfrepcov: Repeated elements in Vector');
+    end
+    q = zeros(1,max(p)+1);
+    q(p+1) = 1;
+
+else
+    
+    q = p;
+    
+end
+endfunction
diff --git a/macros/gftrunc.sci b/macros/gftrunc.sci
new file mode 100644
index 0000000..de0e6e6
--- /dev/null
+++ b/macros/gftrunc.sci
@@ -0,0 +1,49 @@
+function at=gftrunc(a)
+//This function is used to truncate the higher order zeroes in the given polynomial equation
+
+//Calling Sequence
+//AT=GFTRUNC(A)
+    
+//Description
+//A is considered to be matrix that gives the coefficients of polynomial GF(p) in ascending order powers
+//A = [1 2 3] denotes 1 + 2 x + 3 x^2    
+//AT=GFTRUNC(A) returns a matrix which gives the polynomial GF(p) truncating the input matrix 
+//that is if A(i)=0, where i > d + 1, where d is the degree of the polynomial, that zero is removed
+    
+//Examples
+//A= [ 0 0 1 4 0 0] returns [0 0 1 4]
+//c = gftrunc([0 0 1 2 3 0 0 0 4 5 0 1 0 0])
+
+    
+//See also
+//gfadd, gfconv, gfdeconv, gfsub, gftuple
+
+//Authors
+//Pola Lakshmi Priyanka, IIT Bombay//
+
+//*************************************************************************************************************************************//
+// Check number of input arguments
+[outa,inpa]=argn(0) 
+
+if inpa~=1 then
+    error('comm:gftrunc: Invalid number of input arguments')
+end
+
+[row_a,col_a]=size(a);
+if row_a~=1 then
+    error('comm:gftrunc: Input argument should be a row vector')
+end
+    for j=1:col_a
+        if ( abs(a(1,j))~=a(1,j) | floor(a(1,j))~=a(1,j) )
+                error('comm:gftrunc:Elements Of A should be Positive Integers');
+        end
+    end
+at=a;
+for i=col_a:-1:1
+    if a(1,i)~=0 then
+        break;
+    else
+        at=[at(1:i-1)]
+    end
+end
+endfunction
diff --git a/macros/iqcoef2imbal.sci b/macros/iqcoef2imbal.sci
new file mode 100644
index 0000000..2c08c2b
--- /dev/null
+++ b/macros/iqcoef2imbal.sci
@@ -0,0 +1,86 @@
+function [Amp_Imb_DB, Ph_Imb_Deg] = iqcoef2imbal(Comp_Coef)
+//   This function returns the amplitude imbalance and phase imbalance 
+//   that a given compensator coefficient will correct.
+
+//   Calling sequence
+//   [AMP_IMB_DB, PH_IMB_DEG] = IQCOEF2IMBAL(COMP_COEF)
+
+//   Description
+//   [AMP_IMB_DB, PH_IMB_DEG] = IQCOEF2IMBAL(COMP_COEF) returns
+//   the amplitude imbalance and phase imbalance 
+//   that a given compensator coefficient will correct.
+//   Comp_Coef is a scalar or a vector of complex numbers.
+//   AMP_IMB_DB and PH_IMB_DEG are the amplitude imbalance in dB
+//   and the phase imbalance in degrees.
+
+//   Examples
+//   [a_imb_db,ph_imb_deq] = iqcoef2imbal([4 2 complex(-0.1145,0.1297) complex(-0.0013,0.0029)])
+//   disp(a_imb_db,'amplitude imbalance in dB =')
+//   disp(ph_imb_deq,'phase imbalance in degrees=')
+
+//   Bibliography 
+//   http://in.mathworks.com/help/comm/ref/iqcoef2imbal.html
+
+//   See also
+//   iqimbal2coef
+
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+    
+//*************************************************************************************************************************************//
+
+//Input argument check
+[out_a,inp_a]=argn(0);
+
+if (inp_a > 1) | (out_a >2) then
+    error('comm:iqcoef2imbal: Invalid number of arguments')
+end
+
+
+if (or(Comp_Coef==%nan) | or(Comp_Coef==%inf))
+      error('comm:iqcoef2imbal: Input arguments should be finte')
+end
+
+Amp_Imb_DB = zeros(size(Comp_Coef));
+Ph_Imb_Deg = zeros(size(Comp_Coef));
+
+for i = 1:length(Comp_Coef)
+    if imag(Comp_Coef(i)) == 0 // To avoid division by zero
+        c = real(Comp_Coef(i));
+        if abs(c) <= 1
+            Amp_Imb_DB(i) = 20*log10((1-c)/(1+c));
+            Ph_Imb_Deg(i) = 0;
+        else
+            Amp_Imb_DB(i) = 20*log10((c+1)/(c-1));
+            Ph_Imb_Deg(i) = 180;
+        end
+    else
+        R11 = 1 + real(Comp_Coef(i));
+        R22 = 1 - real(Comp_Coef(i));
+        R21 = imag(Comp_Coef(i));
+        R12 = imag(Comp_Coef(i));
+
+        K0 = [R22 -R21; -R12 R11]; 
+        
+        if R11 == 1 
+            a = 0; 
+        else
+            C1 = -K0(1,1)*K0(1,2) + K0(2,2)*K0(2,1);
+            C2 = K0(1,2)^2 + K0(2,1)^2 - K0(1,1)^2 - K0(2,2)^2;
+
+            
+            if abs(Comp_Coef(i)) <= 1
+                a = (-C2 - sqrt(C2^2 + 4*C1^2))/(2*C1);
+            else
+                a = (-C2 + sqrt(C2^2 + 4*C1^2))/(2*C1);
+            end
+        end
+        
+        K = K0 * [1 -a; a 1];
+        
+        Amp_Imb_DB(i) = 20*log10(K(1,1)/K(2,2));
+        Ph_Imb_Deg(i) = -2*atan(K(2,1)/K(1,1))/%pi*180;
+    end
+end
+
+endfunction
diff --git a/macros/iqimbal2coef.sci b/macros/iqimbal2coef.sci
new file mode 100644
index 0000000..2116ca8
--- /dev/null
+++ b/macros/iqimbal2coef.sci
@@ -0,0 +1,68 @@
+function Comp_Coef = iqimbal2coef(Amp_Imb_dB, Ph_Imb_Deg)
+//   This function returns the I/Q imbalance compensator coefficient for given amplitude and phase imbalance.
+
+//   Calling Sequence
+//   COMP_COEF = IQIMBAL2COEF(AMP_IMB_DB, PH_IMB_DEG)
+
+//   Description
+//   COMP_COEF = IQIMBAL2COEF(AMP_IMB_DB, PH_IMB_DEG) returns the I/Q imbalance 
+//   compensator coefficient for given amplitude and phase imbalance.
+//   Comp_Coef is a scalar or a vector of complex numbers.
+//   AMP_IMB_DB and PH_IMB_DEG are the amplitude imbalance in dB
+//   and the phase imbalance in degrees and should be of same size.
+
+//   Examples
+//   [a_imb_db,ph_imb_deg] = iqcoef2imbal([4 2 complex(-0.1145,0.1297) complex(-0.0013,0.0029)])
+//   disp(a_imb_db,'amplitude imbalance in dB =')
+//   disp(ph_imb_deg,'phase imbalance in degrees=')
+//   Comp_Coef = iqimbal2coef(a_imb_db, ph_imb_deg)
+//   disp(Comp_Coef,'Compensator Coefficients=')
+
+//   Bibliography
+//   http://in.mathworks.com/help/comm/ref/iqimbal2coef.html
+
+//   See also
+//   iqcoef2imbal
+
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+//*************************************************************************************************************************************//
+
+//Input argument check
+[out_a,inp_a]=argn(0);
+
+if (inp_a > 2) | (out_a > 1) then
+    error('comm:iqimbal2coef: Invalid number of arguments')
+end
+
+if ( or(Comp_Coef==%nan) | or(Comp_Coef==%inf))
+      error('comm:iqimbal2coef: Input arguments should be finte')
+end
+
+if ( size(Amp_Imb_dB) ~= size(Ph_Imb_Deg) ) then\
+      error('comm:iqimbal2coef: Input arguments should be of same size')
+end
+
+Comp_Coef = complex(zeros(size(Amp_Imb_dB)));
+
+for i = 1:length(Amp_Imb_dB)
+
+    Igain = 10^(0.5*Amp_Imb_dB(i)/20);
+    Qgain = 10^(-0.5*Amp_Imb_dB(i)/20);
+    angle_i = -0.5*Ph_Imb_Deg(i)*%pi/180;
+    angle_q = %pi/2 + 0.5*Ph_Imb_Deg(i)*%pi/180;
+    K = [Igain*cos(angle_i) Qgain*cos(angle_q); ...
+         Igain*sin(angle_i) Qgain*sin(angle_q)];
+
+    R = inv(K);
+
+    w1r = (R(1,1)+R(2,2))/2;
+    w1i = (R(2,1)-R(1,2))/2;
+    w2r = (R(1,1)-R(2,2))/2;
+    w2i = (R(2,1)+R(1,2))/2;
+    w1 = w1r + complex(0,1) * w1i;
+    w2 = w2r + complex(0,1) * w2i;
+
+    Comp_Coef(i) = w2/w1;
+end
diff --git a/macros/iscatastrophic.sci b/macros/iscatastrophic.sci
new file mode 100644
index 0000000..9c54f19
--- /dev/null
+++ b/macros/iscatastrophic.sci
@@ -0,0 +1,91 @@
+function result = iscatastrophic(trellis)
+
+//   This function  determines if a convolutional code is catastrophic or not
+
+//   Calling Sequence
+//   RESULT = ISCATASTROPHIC(TRELLIS)
+
+//   Description
+//   RESULT = ISCATASTROPHIC(TRELLIS) returns 1 if the specified 
+//   trellis corresponds to a catastrophic convolutional code, else 0.
+
+//   Examples 
+//   eg_1.numInputSymbols = 4;
+//   eg_1.numOutputSymbols = 4;
+//   eg_1.numStates = 3;
+//   eg_1.nextStates = [0 1 2 1;0 1 2 1; 0 1 2 1];
+//   eg_1.outputs = [0 0 1 1;1 1 2 1; 1 0 1 1];
+//   res_t_eg_1=istrellis(eg_1)
+//   res_c_eg_1=iscatastrophic(eg_1)
+//   if (res_c_eg_1) then
+//       disp('Example 1 is catastrophic')
+//   else
+//       disp('Example 1 is not catastrophic')
+//   end
+
+//   eg_2.numInputSymbols = 2;
+//   eg_2.numOutputSymbols = 4;
+//   eg_2.numStates = 2;
+//   eg_2.nextStates = [0 0; 1 1 ]
+//   eg_2.outputs = [0 0; 1 1];
+//   res_t_eg_2=istrellis(eg_2)
+//   res_c_eg_2=iscatastrophic(eg_2)
+//   if (res_c_eg_2) then
+//       disp('Example 2 is catastrophic')
+//   else
+//       disp('Example 2 is not catastrophic')
+//   end 
+
+
+//   See also
+//   istrellis
+
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+//*************************************************************************************************************************************//
+
+// Check number of input arguments
+[out_a,inp_a]=argn(0)
+
+if inp_a~=1 then
+    error('comm:iscatastrophic: Invalid number of input arguments')
+end
+
+
+if out_a>1 then
+    error('comm:iscatastrophic: Invalid number of output arguments')
+end
+
+// Check if the input is a valid trellis
+if ~istrellis(trellis),
+    error('comm:iscatastrophic:Input should be a valid trellis structure.');  
+end
+
+result = 0;
+
+// Find indices of zeros in trellis structure
+[r_idx,c_idx] = find(trellis.outputs==0);
+
+//Find Connectivity matrix and check if it is catastrophic
+A = zeros(trellis.numStates,trellis.numStates);
+for k = 2:length(r_idx)
+    A(r_idx(k),trellis.nextStates(r_idx(k),c_idx(k))+1)=1;
+end
+
+
+test = A;
+for i = 1:trellis.numStates
+    for j = 1:trellis.numStates
+        if test(j,j)==1
+            result = 1
+        end
+    end
+    if result==1
+        break
+    else
+        test = test*A;
+    end
+end
+
+endfunction
diff --git a/macros/istrellis.sci b/macros/istrellis.sci
new file mode 100644
index 0000000..da8be0c
--- /dev/null
+++ b/macros/istrellis.sci
@@ -0,0 +1,163 @@
+function [isOk, status]  = istrellis(S)
+
+//    This function checks if the given input is of trellis structure
+
+//    Calling Sequence
+//    [ISOK, STATUS]  = ISTRELLIS(S)
+// 
+//    Description
+//    [ISOK, STATUS]  = ISTRELLIS(S) returns [T,''] if the given input is valid trellis structure. Otherwise ISOK is F and STATUS 
+//    indicates the reason for invalidity
+
+//    Fields in trellis structure are 
+//    numInputSymbols,  (number of input symbols)
+//    numOutputSymbols, (number of output symbols)
+//    numStates,        (number of states)
+//    nextStates,       (next state matrix)
+//    outputs,          (output matrix)
+
+//    Properties of the fields are as follows
+//    numInputSymbols and numOutputSymbols should be a power of 2 (as data is represented in bits).
+//    The 'nextStates' and 'outputs' fields are matrices of size 'numStates' x 'numInputSymbols' .
+//    Each element in the 'nextStates' matrix and 'output' matrix is an integer value between zero and (numStates-1).
+//    The (r,c) element of the 'nextStates' matrix and 'output' matrix,denotes the next state and output respectively when
+//    the starting state is (r-1) and the input bits have decimal representation (c-1).
+
+//    To convert to decimal value, use the first input bit as the most significant bit (MSB).
+
+//    Examples
+//    Valid trellis structure
+//    trellis.numInputSymbols = 4;
+//    trellis.numOutputSymbols = 4;
+//    trellis.numStates = 3;
+//    trellis.nextStates = [0 1 2 1;0 1 2 1; 0 1 2 1];
+//    trellis.outputs = [0 0 1 1;1 1 2 1; 1 0 1 1];
+//    [isok,status] = istrellis(trellis)
+
+//    Inavlid trellis structure
+//    trellis.numInputSymbols = 3;
+//    trellis.numOutputSymbols = 3;
+//    trellis.numStates = 3;
+//    trellis.nextStates = [0 1 2 ;0 1 2 ; 0 1 2 ];
+//    trellis.outputs = [0 0 1 ;1 1 2 ; 1 0 1 ];
+//    [isok,status] = istrellis(trellis)
+
+//    See also
+//    iscatastrophic
+
+//    Authors
+//    Pola Lakshmi Priyanka, IIT Bombay//
+
+//*************************************************************************************************************************************//
+
+// Check arguments
+[out_a,inp_a]=argn(0)
+
+if (inp_a~=1 | out_a>2) then
+    error('comm:istrellis:Invalid number of arguments')
+end
+
+status = '';
+
+// Check if input is a structure
+isOk = isstruct(S);
+if ~isOk then
+  status = string(('comm:trellis:Input is not a structure'));
+  return;
+end
+
+// Check for valid field names
+names  = fieldnames(S);
+numf = size(names, 1);
+
+actual = {'numInputSymbols';'numOutputSymbols';'numStates';'nextStates';'outputs'}
+
+isOk = (numf == 5) & isequal(gsort(names),gsort(actual));
+if ~isOk then
+  status = string(('comm:trellis:Structure is not of trellis type'));
+  return;
+end
+
+// Check for number of Input Symbols.
+numInputSymbols = S.numInputSymbols;
+isOk = isequal(length(numInputSymbols), 1);
+if ~isOk then
+  status = string(('comm:trellis:Number of input symbols is not scalar'));
+  return;
+end
+
+power  = log2(numInputSymbols);
+isOk = isequal(power, double(int32(power)));
+if ~isOk then
+  status = string(('comm:trellis:Number of input symbols is not power of 2'));
+  return;
+end
+
+// Check for number of Output Symbols.
+numOutputSymbols = S.numOutputSymbols;
+isOk = isequal(length(numOutputSymbols), 1);
+if ~isOk then
+  status = string(('comm:trellis:Number of output symbols is not scalar'));
+  return;
+end
+
+power  = log2(numOutputSymbols);
+isOk = isequal(power, double(int32(power)));
+
+if ~isOk then
+  status = string(('comm:trellis:Number of input symbols is not power of 2'));
+  return;
+end
+
+//Check Number of states
+isOk = isequal(length(S.numStates), 1) & ...
+       isequal(S.numStates, double(int32(S.numStates)));
+if isOk then
+  isOk =  S.numStates > 0;
+end
+
+if ~isOk then
+  status = string(('comm:trellis:Number of states is not scalar positive integer'));
+  return;
+end
+
+// Check nextStates
+nextStates = S.nextStates;
+numStates  = S.numStates;
+
+// Check size of nextStates
+s = size(nextStates);
+isOk =  ((length(s) == 2) & (s(1) == numStates) & (s(2) == numInputSymbols) & isequal(double(uint32(nextStates)), nextStates));
+if ~isOk then
+  status = string(('comm:trellis:Next states field matrix size is incorrect'));
+  return;
+end
+
+//Check values of nextState
+isOk = isempty(find(nextStates >= numStates));
+if ~isOk then
+  status = string(('comm:trellis: Elements in next state must be in between 0 to numStates - 1'));
+  return;
+end
+
+//  Check outputs
+outputs  = S.outputs;
+
+// check size of outputs
+s  = size(outputs);
+isOk =  ((length(s) == 2) & (s(1) == numStates) & (s(2)== numInputSymbols) );
+if ~isOk then
+  status = string(('comm:trellis:Outputs field matrix size is incorrect'));
+  return;
+end
+
+
+//Check values of output
+decOutputs = oct2dec(string(outputs));
+isOk = isempty(find(decOutputs >= numOutputSymbols));
+if ~isOk then
+  status = string(('comm:trellis: Elements in output must be in between 0 to numStates - 1 '));
+  return;
+end
+
+endfunction
diff --git a/macros/lteZadoffChuSeq.sci b/macros/lteZadoffChuSeq.sci
new file mode 100644
index 0000000..9f4a6c1
--- /dev/null
+++ b/macros/lteZadoffChuSeq.sci
@@ -0,0 +1,58 @@
+function seq = lteZadoffChuSeq(R, N)
+//   This function generates root Zadoff-Chu sequence of complex symbols as per LTE specifications.
+//
+//   Calling sequence
+//   SEQ = LTEZADOFFCHUSEQ(R, N)
+//
+//   Description
+//   SEQ = LTEZADOFFCHUSEQ(R, N) generates the Rth root Zadoff-Chu sequence (SEQ)
+//   of length N.
+
+//   Examples
+//   seq1 = lteZadoffChuSeq(2, 3)
+//   disp(seq1,'seq1')
+//   Error should occur because inputs are not co primes
+//   seq2 = lteZadoffChuSeq(25, 5)
+//   disp(seq2,'seq2')
+
+
+//   Bibliography 
+//   3rd Generation Partnership Project, Technical Specification Group Radio
+//   Access Network, Evolved Universal Terrestrial Radio Access (E-UTRA),
+//   Physical channels and modulation, Release 10, 3GPP TS 36.211, v10.0.0,
+//   2010-12.
+
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+//*************************************************************************************************************************************//
+
+//Input argument check
+[out_a,inp_a]=argn(0)
+
+if inp_a~=2 then
+    error('comm:lteZadoffChuSeq:Invalid number of input arguments')
+end
+
+if out_a>1 then
+    error('comm:lteZadoffChuSeq:Invalid number of output arguments')
+end
+
+// Error Check for input arguments
+if(~isscalar(R) | ~or(type(R)==[1 5 8]) | ~isreal(R) | R==%inf | R~=floor(R) | R<=0) then
+    error('comm:lteZadoffChuSeq: Input argument 1 should be a finite positive integer')
+end
+
+if(~isscalar(N) | ~or(type(N)==[1 5 8]) | ~isreal(N) | N==%inf | N~=floor(N) | N<=0 | modulo(N,2)==0) then
+    error('comm:lteZadoffChuSeq:Input argument 2 is invalid')
+end
+
+if(gcd(uint8([R N])) ~= 1) then
+    error('comm:lteZadoffChuSeq:Both the input arguments should be co primes')
+end
+
+
+m = (0:N-1).';
+seq = exp( -complex(0,1) * %pi * R * m.*(m+1) / N );
+
+endfunction
diff --git a/macros/ssbdemod.sci b/macros/ssbdemod.sci
new file mode 100644
index 0000000..7271cf9
--- /dev/null
+++ b/macros/ssbdemod.sci
@@ -0,0 +1,160 @@
+function z = ssbdemod(y, Fc, Fs, varargin)
+    
+//   This function performs Single Side Band Amplitude Demodulation
+
+//   Calling Sequence
+//   Z = SSBDEMOD(Y,Fc,Fs) 
+//   Z = SSBDEMOD(Y,Fc,Fs,INI_PHASE) 
+//   Z = SSBDEMOD(Y,Fc,Fs,INI_PHASE,NUM,DEN) 
+
+//   Description
+//   Z = SSBDEMOD(Y,Fc,Fs) 
+//   demodulates the single sideband amplitude modulated signal Y 
+//   with the carrier frequency Fc (Hz).
+//   Sample frequency Fs (Hz). zero initial phase (ini_phase).
+//   The modulated signal can be an upper or lower sideband signal. 
+//   A lowpass butterworth filter is used in the demodulation.  
+// 
+//   Z = SSBDEMOD(Y,Fc,Fs,INI_PHASE) 
+//   adds an extra argument the initial phase (rad) of the modulated signal.
+// 
+//   Z = SSBDEMOD(Y,Fc,Fs,INI_PHASE,NUM,DEN) 
+//   adds extra arguments about the filter specifications 
+//   i.e., the numerator and denominator of the lowpass filter.
+//
+//   Fs must satisfy Fs >2*(Fc + BW), where BW is the bandwidth of the
+//   modulating signal.
+ 
+
+//   Examples
+
+//   Fs =200;
+//   t = [0:2*Fs+1]'/Fs;
+//   ini_phase = 5;
+//   Fc = 20;
+//   fm1= 2;
+//   fm2= 3
+//   x =sin(2*fm1*%pi*t)+sin(2*fm2*%pi*t);
+//   y = ssbmod(x,Fc,Fs,ini_phase);
+//   o = ssbdemod(y,Fc,Fs,ini_phase);
+//   z = fft(y);
+//   zz =abs(z(1:length(z)/2+1 ));
+//   axis = (0:Fs/length(zz):Fs -(Fs/length(zz)))/2;
+//
+//   figure
+//   subplot(3,1,1); plot(x);
+//   title(' Message signal');
+//   subplot(3,1,2); plot(y);
+//   title('Amplitude modulated signal');
+//   subplot(3,1,3); plot(axis,zz);
+//   title('Spectrum of amplitude modulated signal');
+//   z1 =fft(o);
+//   zz1 =abs(z1(1:length(z1)/2+1 ));
+//   axis = (0:Fs/length(zz1):Fs -(Fs/length(zz1)))/2;
+//   figure
+//   subplot(3,1,1); plot(y);
+//   title(' Modulated signal');
+//   subplot(3,1,2); plot(o);
+//   title('Demodulated signal');
+//   subplot(3,1,3); plot(axis,zz1);
+//   title('Spectrum of Demodulated signal');
+
+//   See also 
+//   ssbmod
+
+//   Authors
+//   Pola Lakshmi Priyanka, IIT Bombay//
+
+//*************************************************************************************************************************************//
+
+//   Check number of input arguments
+[outa,inpa]=argn(0) 
+if(inpa > 6)
+    error("comm:ssbdemod:Too Many Input Arguments");
+end
+//funcprot(0) //to protect the function 
+
+//Check y,Fc, Fs.
+if(~isreal(y)| ~or(type(y)==[1 5 8]) )
+    error("comm:ssbdemod: Y must be real");
+end
+
+if(~isreal(Fc) | ~isscalar(Fc) | Fc<=0 | ~or(type(Fc)==[1 5 8]) )
+    error("comm:ssbdemod:Fc must be Real, scalar, positive");
+end
+
+if(~isreal(Fs) | ~isscalar(Fs) | Fs<=0 | ~or(type(Fs)==[1 5 8]) )
+    error("comm:ssbdemod:Fs must be Real, scalar, positive");
+end
+
+// Check if Fs is greater than 2*Fc
+if(Fs<=2*Fc)
+    error("comm:ssbdemod:Fs<2Fc:Nyquist criteria");
+end
+
+// Check initial phase
+
+if(inpa<4 )
+    ini_phase = 0;
+else 
+    ini_phase = varargin(1);
+end
+if(~isreal(ini_phase) | ~isscalar(ini_phase)| ~or(type(ini_phase)==[1 5 8]) )
+    error("comm:ssbdemod:Initial phase shoould be Real");
+end
+
+// Filter specifications
+if(inpa<5)  
+    H = iir(5,'lp','butt',[Fc/Fs,0],[0,0]); 
+    
+    num = coeff(H(2));
+    den = coeff(H(3));
+    num = num(length(num):-1:1);
+    den = den(length(den):-1:1);
+    
+// Check that the numerator and denominator are valid, and come in a pair
+elseif( (inpa == 5) )
+    error("comm:ssbdemod:NumDenPair: Filter error :Two arguments required");
+
+// Check to make sure that both num and den have values
+elseif( bitxor( isempty(varargin(1)), isempty(varargin(2))))
+    error(message('comm:ssbdemod:Filter specifications'));
+elseif(  isempty(varargin(1)) & isempty(varargin(2)) ) 
+    H = iir(7,'lp','butt',[Fc/Fs*2*%pi,0],[0,0]); 
+    
+    num = coeff(H(2));
+    den = coeff(H(3));
+    num = num(length(num):-1:1);
+    den = den(length(den):-1:1);
+else 
+    num = varargin(1);
+    den = varargin(2);
+end
+
+
+
+// check if Y is one dimensional 
+wid = size(y,1);
+if(wid ==1)
+    y = y(:);
+end
+
+// Demodulation
+t = (0 : 1/Fs :(size(y,1)-1)/Fs)';
+t = t(:, ones(1, size(y, 2)));
+z = y .* cos(2*%pi * Fc * t + ini_phase);
+for i = 1 : size(z, 2)
+    z(:, i) = filter(num, den, z(:, i)) ;
+    z=z(length(z):-1:1)
+    z(:, i) = filter(num, den, z(:, i)) ;
+    z=z*-2;
+end;
+
+// restore the output signal to the original orientation 
+if(wid == 1)
+    z = z';
+end
+
+endfunction
+
+// End of function