diff options
Diffstat (limited to 'macros')
| -rw-r--r-- | macros/alignsignals.sci | 166 | ||||
| -rw-r--r-- | macros/arithdeco.sci | 170 | ||||
| -rw-r--r-- | macros/arithenco.sci | 166 | ||||
| -rw-r--r-- | macros/build_help.sce | 5 | ||||
| -rw-r--r-- | macros/finddelay.sci | 187 | ||||
| -rw-r--r-- | macros/gfcosets.sci | 93 | ||||
| -rw-r--r-- | macros/gflineq.sci | 170 | ||||
| -rw-r--r-- | macros/gfrepcov.sci | 60 | ||||
| -rw-r--r-- | macros/gftrunc.sci | 49 | ||||
| -rw-r--r-- | macros/iqcoef2imbal.sci | 86 | ||||
| -rw-r--r-- | macros/iqimbal2coef.sci | 68 | ||||
| -rw-r--r-- | macros/iscatastrophic.sci | 91 | ||||
| -rw-r--r-- | macros/istrellis.sci | 163 | ||||
| -rw-r--r-- | macros/lteZadoffChuSeq.sci | 58 | ||||
| -rw-r--r-- | macros/ssbdemod.sci | 160 |
15 files changed, 1692 insertions, 0 deletions
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 |