code by jitendra

13 files changed, 323 insertions, 147 deletions

diff --git a/macros/diric.sci b/macros/diric.sci
index 077bb13..f60fee8 100644
--- a/macros/diric.sci
+++ b/macros/diric.sci
@@ -7,7 +7,7 @@ function [y]= diric(x,n)
 //	x: Real valued vector or matrix
 //	n: Real positive integer or complex integer
 // Description
-//	This is an Octave function
+
 //	y=diric(x,n) returns the dirichlet function values of parameter x.
 // Examples
 // 1. 	diric([1 2 3],3)
@@ -17,8 +17,14 @@ function [y]= diric(x,n)
 
 funcprot(0);
 rhs=argn(2);
-if (rhs~=2) then
+if (argn(2)~=2) then
 	error ("Wrong number of input arguments.")
-else y= callOctave("diric",x,n)
+elseif (n <= 0 | ceil(n) ~= n) then
+          error("n must be an positive integer.");
+          
+   else       
+    y = sin(n.*x./2)./(n.*sin(x./2));
+  y(pmodulo(x,2*%pi)==0) = (-1).^((n-1).*x(pmodulo(x,2*%pi)==0)./(2.*%pi));
 end
+
 endfunction
diff --git a/macros/filtic.sci b/macros/filtic.sci
index efd48ba..56367b7 100644
--- a/macros/filtic.sci
+++ b/macros/filtic.sci
@@ -1,5 +1,4 @@
-function zf = filtic (b, a, y, x)
-
+function zf = filtic(b,a,y,x)
 //This function finds the initial conditions for the delays in the transposed direct-form II filter implementation
 //Calling Sequence
 //zf = filtic (b, a, y)
@@ -18,19 +17,36 @@ function zf = filtic (b, a, y, x)
 //    0.00000 - 22.60000i
 //    2.40000 +  0.00000i
 //    0.00000 +  0.00000i
-//This function is being called from Octave
 
-funcprot(0);
-rhs = argn(2)
+  if (argn(2)>4 | argn(2)<3) | (argn(1)>1)
+    error("Wrong number of input agruments.")
+  end
+  if argn(2) < 4, x = []; end
+
+  nz = max(length(a)-1,length(b)-1);
+  zf=zeros(nz,1);
 
-if(rhs>4 | rhs<3)
-	error("Wrong number of input agruments.")
-end
+  
+  if length(a)<(nz+1)
+    a(length(a)+1:nz+1)=0;
+  end
+  if length(b)<(nz+1)
+    b(length(b)+1:nz+1)=0;
+  end
+ 
+  if length(x) < nz
+    x(length(x)+1:nz)=0;
+  end
+  if length(y) < nz
+    y(length(y)+1:nz)=0;
+  end
+
+  for i=nz:-1:1
+    for j=i:nz-1
+      zf(j) = b(j+1)*x(i) - a(j+1)*y(i)+zf(j+1);
+    end
+    zf(nz)=b(nz+1)*x(i)-a(nz+1)*y(i);
+  end
 
-select(rhs)
-case 3 then
-zf = callOctave("filtic",b,a,y)
-case 4 then
-zf = callOctave("filtic",b,a,y,x)
-end
 endfunction
+
diff --git a/macros/fractdiff.sci b/macros/fractdiff.sci
index 979a079..fda26be 100644
--- a/macros/fractdiff.sci
+++ b/macros/fractdiff.sci
@@ -3,15 +3,44 @@ function y= fractdiff(x,d)
 //Calling Sequence
 // fractdiff (X, D)
 //Description
-//This is an Octave function.
 //Compute the fractional differences (1-L)^d x where L denotes the lag-operator and d is greater than -1.
-	funcprot(0);
-	rhs= argn(2);
-	if(rhs < 2 | rhs >2)
-		error("Wrong number of input arguments");
-	end
-	select(rhs)
-	case 2 then
-		y= callOctave("fractdiff",x,d);
-	end
-endfunction
\ No newline at end of file
+
+if(argn(2)~=2)
+	error("Wrong number of input arguments");
+end
+
+
+N = 100;
+
+  if (~ isvector (x) | ~isscalar(d))
+    error ("X must be a vector and D must be a scalar");
+  end
+  
+  if (d<= -1) then
+      error ("D must be > -1");
+  end
+  
+  if (d >= 1)
+    for k = 1 : d
+      x = x(2 : length (x)) - x(1 : length (x) - 1);
+    end
+  end
+//disp(x)
+  if (d >-1)
+       d=d-fix(d./1).*1;
+//disp(d)
+    if (d ~= 0)
+      n = (0 : N)';
+      w = real (gamma (-d+n) ./ gamma (-d) ./ gamma (n+1));
+    //  disp(w); disp(x)
+      y = fftfilt (w, x);
+      y = y(1 : length (x));
+    else
+             // disp(x)
+      y = x;
+    end
+
+  end
+  
+  
+endfunction
diff --git a/macros/gauspuls.sci b/macros/gauspuls.sci
index ad853ee..621670b 100644
--- a/macros/gauspuls.sci
+++ b/macros/gauspuls.sci
@@ -1,5 +1,4 @@
 function [y]=gauspuls(t,fc,bw)
-
 // Generates Gaussian-modulated sinusoidal pulses
 // Calling Sequence
 //	[y]=gauspuls(t,fc,bw)
@@ -10,7 +9,6 @@ function [y]=gauspuls(t,fc,bw)
 //	fc: Real non negative number or complex number
 //	bw: Real positive number or complex number
 // Description
-//	This is an Octave function
 //	This function returns a Gaussian RF pulse of unity amplitude at the times indicated in array t.
 // Examples
 // 1.	gauspuls(1,2,3)
@@ -18,14 +16,21 @@ function [y]=gauspuls(t,fc,bw)
 // 2.	gauspuls([1 2 3],1,1)
 //	ans= 0.0281016    0.0000006    1.093D-14
 
-funcprot(0);
-rhs=argn(2);
-if ( rhs<1 ) then
- 	error ("Wrong number of input arguments.")
-elseif (rhs==1)
-	y= callOctave("gauspuls",t)
-elseif (rhs==2)
-	y= callOctave("gauspuls",t,fc)
-else y= callOctave("gauspuls",t,fc,bw)
-end
+if ( argn(2)<1 | argn(2)>3 ) then
+ error ("Wrong number of input arguments.")
+elseif (argn(2)==1)
+    fc = 1e3; bw = 0.5;
+elseif (argn(2)==2)
+     bw = 0.5;
+ end
+ 
+ 
+ if (~isscalar(fc) | ~isreal(fc) | fc<0) then
+     error('fc must be non-negative real scalar')
+ end
+  if (~isscalar(bw) | ~isreal(bw) | bw<=0) then
+     error('bw must be positive real scalar')
+ end
+  y = exp (-t .* t / (2*(1 / (4*%pi^2 * (-(bw^2 * fc^2) / (8 * log (10 ^ (-6/20)))))))) .* cos (2*%pi*fc * t);
+
 endfunction
diff --git a/macros/kaiserord.sci b/macros/kaiserord.sci
index 47f4105..5f46193 100644
--- a/macros/kaiserord.sci
+++ b/macros/kaiserord.sci
@@ -10,7 +10,6 @@ function [n, Wn, beta, ftype] = kaiserord (f, m, dev, fs)
 //dev: Deviation of the filter.
 //fs: Sampling rate.
 //Description
-//This is an Octave function.
 //The vector f contains pairs of frequency band edges in the range [0,1]. The vector m specifies the magnitude response for each band. The values of m must be zero for all stop bands and must have the
 //same magnitude for all pass bands. The deviation of the filter dev can be specified as a scalar or a vector of the same length as m. The optional sampling rate fs can be used to indicate that f is in
 //Hz in the range [0,fs/2].
@@ -28,43 +27,70 @@ function [n, Wn, beta, ftype] = kaiserord (f, m, dev, fs)
 //    n = (A-8)/(2.285 dw)
 //Examples
 //[n, w, beta, ftype] = kaiserord ([1000, 1200], [1, 0], [0.05, 0.05], 11025)
-//n =  1
-//w =  1100
+//n =  70
+//w =  0.199
 //beta =  1.5099
 //ftype = low
 
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 3 | rhs > 4)
+if (argn(2) < 3 | argn(2) > 4)
 error("Wrong number of input arguments.")
 end
 
-select(rhs)
-	
-	case 3 then
-		if(lhs==1)
-		n = callOctave("kaiserord",f, m, dev)
-		elseif(lhs==2)
-		[n, Wn] = callOctave("kaiserord",f, m, dev)
-		elseif(lhs==3)
-		[n, Wn, beta] = callOctave("kaiserord",f, m, dev)
-		elseif(lhs==4)
-		[n, Wn, beta, ftype] = callOctave("kaiserord",f, m, dev)
-		else
-		error("Wrong number of output argments.")
-	       	end
-	case 4 then
-		if(lhs==1)
-		n = callOctave("kaiserord",f, m, dev, fs)
-		elseif(lhs==2)
-		[n, Wn] = callOctave("kaiserord",f, m, dev, fs)
-		elseif(lhs==3)
-		[n, Wn, beta] = callOctave("kaiserord",f, m, dev, fs)
-		elseif(lhs==4)
-		[n, Wn, beta, ftype] = callOctave("kaiserord",f, m, dev, fs)
-		else
-		error("Wrong number of output argments.")
-	       	end
-	end
+
+  if argn(2)<4 then 
+      fs=2; 
+      end
+
+  if length(f)~=(2*length(m)-2) then
+    error("kaiserord must have one magnitude for each frequency band");
+  end
+  
+  if or(m(1:length(m)-2)~=m(3:length(m))) then
+    error("kaiserord pass and stop bands must be strictly alternating");
+  end
+  if length(dev)~=length(m) & length(dev)~=1 then
+    error("kaiserord must have one deviation for each frequency band");
+  end
+  dev = min(dev);
+  if dev <= 0 then
+      error("kaiserord must have dev>0"); 
+      end
+  Wn = (f(1:2:length(f))+f(2:2:length(f)))/fs;
+
+  if length(Wn) == 1 then
+    if m(1)>m(2) then
+        ftype='low'; 
+        else 
+        ftype='high';
+         end
+  elseif length(w) == 2 then
+    if m(1)>m(2) then 
+        ftype='stop'; 
+        else 
+        ftype='pass'; 
+        end
+  else
+    if m(1)>m(2),
+         ftype='DC-1'; 
+        else 
+        ftype='DC-0'; 
+        end
+  end
+
+  A = -20*log10(dev);
+  if (A > 50)
+    beta = 0.1102*(A-8.7);
+  elseif (A >= 21)
+    beta = 0.5842*(A-21)^0.4 + 0.07886*(A-21);
+  else
+    beta = 0.0;
+  end
+
+  dw = 2*%pi*min(f(2:2:length(f))-f(1:2:length(f)))/fs;
+  n = max(1,ceil((A-8)/(2.285*dw)));
+  
+  if ((m(1)>m(2)) == ((length(Wn)-fix(length(Wn)./2).*2)==0)) & (n-fix(n./2).*2)==1 then
+      n = n+1; 
+      end
+
 endfunction
diff --git a/macros/meyeraux.sci b/macros/meyeraux.sci
index 657e829..7cde1b9 100644
--- a/macros/meyeraux.sci
+++ b/macros/meyeraux.sci
@@ -1,12 +1,10 @@
 function [y]=meyeraux(x)
-
 // Returns value of Meyer Wavelet Auxiliary function
 // Calling Sequence
 //	[y]=meyeraux(x)
 // Parameters
 //	x: Real or complex valued vector or matrix
 // Description
-//	This is an Octave function.
 //	This function returns values of the auxiliary function used for Meyer wavelet generation.
 // Examples
 // 1.	meyeraux([1 2 3])
@@ -14,10 +12,11 @@ function [y]=meyeraux(x)
 // 2.	meyeraux([1 2 3;4 5 6])
 //	ans=  [1      -208    -10287  ;	 -118016   -709375  -2940624 ]
 
-funcprot(0);
-rhs=argn(2);
-if (rhs~=1) then
+if (argn(2)~=1) then
 	error ("Wrong number of input arguments.")
-else y=callOctave("meyeraux",x)
+else 
+  y = 35.*x.^4-84.*x.^5+70.*x.^6-20.*x.^7;
 end
 endfunction
+
+
diff --git a/macros/rceps.sci b/macros/rceps.sci
index 8b4d89a..f0433a3 100644
--- a/macros/rceps.sci
+++ b/macros/rceps.sci
@@ -12,18 +12,37 @@ function [y, xm]= rceps(x)
 // s(1:Fs/f0:length(s)) = 1;              # Impulse glottal wave
 // x = filter (1, a, s);                  # Speech signal in x
 // [y, xm] = rceps (x .* hanning (1024)); # cepstrum and min phase reconstruction
-funcprot(0)
-lhs= argn(1)
-rhs= argn(2)
 
-if(rhs <1 | rhs> 1 )
+if(argn(2)~= 1 )
 error("Wrong number of Input Arguments");
 end
 
-if(lhs<2 | lhs>2)
+if(argn(1)>2)
 error("Wrong number of Output Arguments")
 end
-
-	[y,xm]= callOctave("rceps",x);
-
+  f = abs(fft1(x));
+  if (or(f == 0))
+    error ("The spectrum of x contains zeros, unable to compute real cepstrum");
+  end
+   
+  y = real(ifft1(log(f)));
+  
+  if argn(1) == 2 then
+    n=length(x);
+    if size(x,1)==1 then
+      if (n-fix(n./2).*2) ==1 then
+        xm = [y(1), 2*y(2:n/2+1), zeros(1,n/2)];
+      else
+        xm = [y(1), 2*y(2:n/2), y(n/2+1), zeros(1,n/2-1)];
+      end
+    else
+      if (n-fix(n./2).*2)==1
+        xm = [y(1,:); 2*y(2:n/2+1,:); zeros(n/2,size(y,2))];
+      else
+        xm = [y(1,:); 2*y(2:n/2,:); y(n/2+1,:); zeros(n/2-1,size(y,2))];
+      end
+    end
+    xm = real(ifft1(exp(fft1(xm)')));
+  end
+  
 endfunction
diff --git a/macros/rectpuls.sci b/macros/rectpuls.sci
index 2796ef0..27138d0 100644
--- a/macros/rectpuls.sci
+++ b/macros/rectpuls.sci
@@ -1,5 +1,4 @@
 function [y]=rectpuls(t,w)
-
 // Generates a Rectangular pulse based on the width and sampling times.
 // Calling Sequence
 //	[y]=rectpuls(t)
@@ -8,7 +7,6 @@ function [y]=rectpuls(t,w)
 //	t: Real or complex valued vector or matrix
 //	w: Real or complex valued vector or matrix
 // Description
-//	This is an Octave function
 //	y = rectpuls(t) returns a continuous, aperiodic, unity-height rectangular pulse depending upon input t, centered about t=0 and having default width of 1.
 //	y = rectpuls(t,w) generates a rectangle of width w.
 // Examples
@@ -17,12 +15,17 @@ function [y]=rectpuls(t,w)
 // 2.	rectpuls([1000 1000 100 100])
 //	ans =   0   0   0   0
 
-funcprot(0);
-rhs=argn(2);
-if (rhs<1) then
+if (argn(2)<1 | argn(2)>2) then
 	error ("Wrong number of input arguments.")
-elseif (rhs==1)
-	y=callOctave("rectpuls",t)
-else y=callOctave("rectpuls",t,w)
+elseif (argn(2)==1) then
+	w=1;
 end
+
+  if (~ isreal (w) | ~ isscalar (w))
+    error ("W must be a real scalar");
+  end
+  
+y = zeros (size(t,1), size(t,2));
+idx = find ((t >= -w/2) & (t < w/2));
+y(idx) = 1;
 endfunction
diff --git a/macros/sawtooth.sci b/macros/sawtooth.sci
index 50d093d..8c3d689 100644
--- a/macros/sawtooth.sci
+++ b/macros/sawtooth.sci
@@ -8,7 +8,6 @@ function [y]=sawtooth (t,width)
 //	t: Real valued vector or matrix
 //	width: Real number between 0 and 1
 // Description
-//	This is an Octave function
 //	This function returns a sawtooth wave with period 2*pi with +1/-1 as the maximum and minimum values for elements of t. If width is specified, it determines where the maximum is in the interval [0,2*pi].
 // Examples
 // 1.	sawtooth([1 2 3 4 5],0.5)
@@ -16,12 +15,21 @@ function [y]=sawtooth (t,width)
 // 2.	sawtooth([1 2; 4 5])
 //	ans =  [-0.68169  -0.36338;   0.27324   0.59155]
 
-funcprot(0);
-rhs=argn(2);
-if (rhs<1) then
+if (argn(2)<1 | argn(2)>2) then
 	error ("Wrong number of input arguments.")
-elseif (rhs==1)
-	y=callOctave("sawtooth",t)
-else	y=callOctave("sawtooth",t,width)
+elseif (argn(2)==1)
+	width=1
 end
+  if (width < 0 | width > 1 | ~ isreal (width))
+      error ("width must be a real value between 0 and 1");
+    end  
+ t = pmodulo (t / (2 * %pi), 1);
+  y = zeros (size (t,1), size(t,2));
+
+  if (width ~= 0)
+    y (t < width) = 2 * t (t < width) / width - 1;
+  end
+  if (width ~= 1)
+    y( t >= width) = -2 * (t (t >= width) - width) / (1 - width) + 1;
+  end 
 endfunction
diff --git a/macros/sgolay.sci b/macros/sgolay.sci
index 1156d40..09062a7 100644
--- a/macros/sgolay.sci
+++ b/macros/sgolay.sci
@@ -10,7 +10,6 @@ function F = sgolay (p, n, m, ts)
 //m: positive integer less than 2^31 or logical
 //ts: real or complex value
 //Description
-//This is an Octave function.
 //This function computes the filter coefficients for all Savitzsky-Golay smoothing filters of order p for length n (odd). 
 //m can be used in order to get directly the mth derivative; ts is a scaling factor.
 //Examples
@@ -20,21 +19,32 @@ function F = sgolay (p, n, m, ts)
 //   0.33333   0.33333   0.33333
 //  -0.16667   0.33333   0.83333
 
-funcprot(0);
-rhs = argn(2)
-
-if(rhs<2 | rhs>4)
+if(argn(2)<2 | argn(2)>4) then
 error("Wrong number of input arguments.")
+elseif ((n-fix(n./2).*2)~=1)  then
+error ("sgolay needs an odd filter length n");
+elseif (p>=n) then
+     error ("sgolay needs filter length n larger than polynomial order p");
+ end
+
+if (argn(2)==2) then
+    m=0; ts=1;
+end
+if (argn(2)==3) then
+    ts=1;
 end
 
-	select(rhs)
-	case 2 then
-	F = callOctave("sgolay",p,n)
-	case 3 then
-	F = callOctave("sgolay",p,n,m)
-	case 4 then
-	F = callOctave("sgolay",p,n,m,ts)
-	end
-endfunction
+if length(m) > 1, error("weight vector unimplemented"); end
+
+ F = zeros (n, n);
+ k = floor (n/2);
+    for row = 1:k+1
+      C = ( [(1:n)-row]'*ones(1,p+1) ) .^ ( ones(n,1)*[0:p] );   
+      A = pinv(C);   
+      F(row,:) = A(1+m,:);
+    end
 
+  F(k+2:n,:) = (-1)^m*F(k:-1:1,n:-1:1);
+  F =  F * ( prod(1:m) / (ts^m) );
 
+endfunction
diff --git a/macros/sos2tf.sci b/macros/sos2tf.sci
index d38a4f7..f674c85 100644
--- a/macros/sos2tf.sci
+++ b/macros/sos2tf.sci
@@ -1,4 +1,4 @@
-function [B,A] = sos2tf(sos, g)
+function [A,B] = sos2tf(sos, g)
 //This function converts series second-order sections to direct H(z) = B(z)/A(z) form.
 //Calling Sequence
 //[B] = sos2tf(sos)
@@ -8,7 +8,6 @@ function [B,A] = sos2tf(sos, g)
 //sos: matrix of real or complex numbers
 //g: real or complex value, default value is 1
 //Description
-//This is an Octave function.
 //This function converts series second-order sections to direct H(z) = B(z)/A(z) form.
 //The input is the sos matrix and the second parameter is the overall gain, default value of which is 1. 
 //The output is a vector.
@@ -19,18 +18,40 @@ function [B,A] = sos2tf(sos, g)
 //  -2   1  2  4  1
 //b =
 //   1 10  0 -10 -1
-funcprot(0);
-rhs = argn(2)
-if(rhs<1 | rhs>2)
+if(argn(2)<1 | argn(2)>2)
 error("Wrong number of input arguments.")
 end
 
+if argn(2)==1 then
+    g=1;
+end
+
+  [N,M] = size(sos);
+
+  if M~=6
+    error('sos matrix must be N by 6');
+  end
+
+  A = 1;
+  B = 1;
+
+  for i=1:N
+    A = conv(A, sos(i,1:3));
+    B = conv(B, sos(i,4:6));
+  end
+
+  nB = length(A);
+  
+  while nB & A(nB)==0
+    A=A(1:nB-1);
+    nB=length(A);
+  end
 
+  nA = length(B);
+  while nA & B(nA)==0
+    B=B(1:nA-1);
+    nA=length(B);
+  end
+  A = A * g;
 
-	select(rhs)
-	case 1 then
-	[B,A] = callOctave("sos2tf",sos)
-	case 2 then
-	[B,A] = callOctave("sos2tf",sos,g)
-	end
 endfunction
diff --git a/macros/sos2zp.sci b/macros/sos2zp.sci
index 98abc36..1f0f335 100644
--- a/macros/sos2zp.sci
+++ b/macros/sos2zp.sci
@@ -11,7 +11,6 @@ function [z,p,k] = sos2zp (sos, g)
 //z: column vector
 //p: column vector
 //Description
-//This is an Octave function. 
 //This function converts series second-order sections to zeros, poles, and gains (pole residues).
 //The input is the sos matrix and the second parameter is the overall gain, default value of which is 1.
 //The outputs are z, p, k. z and p are column vectors containing zeros and poles respectively, and k is the overall gain. 
@@ -24,18 +23,32 @@ function [z,p,k] = sos2zp (sos, g)
 //  -0.6250 + 1.0533i
 //  -0.6250 - 1.0533i
 //c =  1
-
-funcprot(0);
-rhs = argn(2)
-if(rhs<1 | rhs>2)
+if(argn(2)<1 | argn(2)>2)
 error("Wrong number of input arguments.")
 end
+if argn(2)==1 then
+    g=1;
+end
+  gns = sos(:,1); 
+  k = prod(gns)*g; 
+  if k==0 then
+     error('one or more section gains is zero'); 
+     end
+  sos(:,1:3) = sos(:,1:3)./ [gns gns gns];
 
-	select(rhs)
-	case 1 then
-	[z,p,k] = callOctave("sos2zp",sos)
-	case 2 then
-	[z,p,k] = callOctave("sos2zp",sos,g)
-	end
-endfunction
+  [N,m] = size(sos);
+  if m~=6 then
+      error('sos matrix should be N by 6'); 
+      end
 
+  z = zeros(2*N,1);
+  p = zeros(2*N,1);
+  for i=1:N
+    ndx = [2*i-1:2*i];
+    zi = roots(sos(i,1:3));
+    z(ndx) = zi;
+    pi = roots(sos(i,4:6));
+    p(ndx) = pi;
+  end
+
+endfunction
diff --git a/macros/welchwin.sci b/macros/welchwin.sci
index a871a59..9b5fac1 100644
--- a/macros/welchwin.sci
+++ b/macros/welchwin.sci
@@ -23,15 +23,36 @@ function w = welchwin (m, opt)
 
 
 rhs = argn(2)
-if(rhs<1 | rhs>2)
+if(argn(2)<1 | argn(2)>2)
 error("Wrong number of input arguments.")
 end
 
-	select(rhs)
-	case 1 then
-	w = callOctave("welchwin",m)
-	case 2 then
-	w = callOctave("welchwin",m,opt)
-	end
+ if ((~isscalar (m) & (m == fix (m)) & (m > 0))) then
+    error (" M must be a positive integer");
+  end
+
+
+  N = (m - 1) / 2;
+  mmin = 3;
+if argn(2)==2 then
+       
+    select (opt)
+      case "periodic"
+        N = m / 2;
+        mmin = 2;
+      case "symmetric"
+        N = (m - 1) / 2;
+      else
+        error ("window type must be either periodic or symmetric");
+    end
+  end
+  
+
+  if (m < mmin)
+    error ('M must be an integer greater than mmin');
+  end
+
+  n = [0:m-1]';
+  w = 1 - ((n-N)./N).^2;
 
 endfunction