31 files changed, 2721 insertions, 1081 deletions
diff --git a/macros/arch_fit.sci b/macros/arch_fit.sci
index 3e431a6..2fa462c 100644
--- a/macros/arch_fit.sci
+++ b/macros/arch_fit.sci
@@ -1,35 +1,85 @@
-function [A, B] = arch_fit(Y, varargin)
-//This functions fits an ARCH regression model to the time series Y using the scoring algorithm in Engle's original ARCH paper.
-//Calling Sequence
-//[A, B] = arch_fit (Y, X, P, ITER, GAMMA, A0, B0)
-//Parameters
-//Description
-//Fit an ARCH regression model to the time series Y using the scoring algorithm in Engle's original ARCH paper.
-//
-//The model is
-//
-//          y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t),
-//          h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(p+1) * e(t-p)^2
-//
-//in which e(t) is N(0, h(t)), given a time-series vector Y up to time t-1 and a matrix of (ordinary) regressors X up to t. The order of the regression of the residual variance is specified by P.
-//
-//If invoked as 'arch_fit (Y, K, P)' with a positive integer K, fit an ARCH(K, P) process, i.e., do the above with the t-th row of X given by
-//
-//          [1, y(t-1), ..., y(t-k)]
-//
-//Optionally, one can specify the number of iterations ITER, the updating factor GAMMA, and initial values a0 and b0 for the scoring algorithm.
-funcprot(0);
-rhs = argn(2);
-lhs=argn(1);
-if(rhs<7 | rhs>7)
-error("Wrong number of input arguments.");
-end
-if (lhs<2 | lhs>2)
-    error("Wrong number of output arguments.");
-end
-
-	select(rhs)
-	case 7 then
-	[A, B] = callOctave("arch_fit",Y, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6));
-	end
-endfunction
+/*
+Dependencies : ols, autoreg_matrix
+Calling Sequence
+      [a, b] = arch_fit (y, x, p)
+      [a, b] = arch_fit (y, x, p, iter, gamma, a0, b0)
+Parameters
+      y(vector) : A time-series data vector up to time t-1 .
+      x (Matrix): A matrix of (ordinary) regressors x up to t.
+      p (scalar): The order of the regression of the residual variance.
+      iter (scaler) : Number of iterations
+      gamma (real number) : updating factor
+      a0 b0 (real numbers) : Initial values for the scoring algorithm
+Description:        
+        Fit an ARCH regression model to the time series y using the scoring algorithm in Engle’s original ARCH paper.
+        The model is
+        y(t) = b(1) * x(t,1) + … + b(k) * x(t,k) + e(t),
+        h(t) = a(1) + a(2) * e(t-1)^2 + … + a(p+1) * e(t-p)^2
+        in which e(t) is N(0, h(t)), given a time-series vector y up to time t-1 and a matrix of (ordinary) regressors x upto t. The order of the regression of the residual variance is specified by p.
+        If invoked as arch_fit (y, k, p) with a positive integer k, fit an ARCH(k, p) process, i.e., do the above with the t-th row of x given by
+        [1, y(t-1), …, y(t-k)]
+        Optionally, one can specify the number of iterations iter, the updating factor gamma, and initial values a0 and b0 for the scoring algorithm.
+*/
+function [a, b] = arch_fit (y, x, p, iter, gamma, a0, b0)
+    nargin = argn(2)
+    if (nargin < 3 || nargin == 6)
+      error("invalid inputs");
+    end
+    if (~ (isvector (y)))
+      error ("arch_fit: Y must be a vector");
+    end
+    T = max(size(y));
+    y = matrix (y, T, 1);
+    [rx, cx] = size (x);
+    if ((rx == 1) && (cx == 1))
+      x = autoreg_matrix (y, x);
+    elseif (~ (rx == T))
+      error ("arch_fit: either rows (X) == length (Y), or X is a scalar");
+    end
+    [T, k] = size (x);
+    if (nargin == 7)
+      a = a0;
+      b = b0;
+      e = y - x * b;
+    else
+      [b, v_b, e] = ols (y, x);
+      zer = zeros(1,p);
+      a = [v_b zer]';
+      if (nargin < 5)
+        gamma = 0.1;
+        if (nargin < 4)
+          iter = 50;
+        end
+      end
+    end
+    esq = e.^2;
+    Z = autoreg_matrix (esq, p);
+    for i = 1 : iter
+      h   = Z * a;
+      tmp = esq ./ h.^2 - 1 ./ h;
+      s   = 1 ./ h(1:T-p);
+      for j = 1 : p
+        s = s - a(j+1) * tmp(j+1:T-p+j);
+      end
+      r = 1 ./ h(1:T-p);
+      for j = 1:p
+        r = r + 2 * h(j+1:T-p+j).^2 .* esq(1:T-p);
+      end
+      r = sqrt (r);
+      X_tilde = x(1:T-p, :) .* (r * ones (1,k));
+      e_tilde = e(1:T-p) .*s ./ r;
+      delta_b = inv (X_tilde' * X_tilde) * X_tilde' * e_tilde;
+      b  = b + gamma * delta_b;  
+      e   = y - x * b;
+      esq = e .^ 2;
+      if isempty(esq) then
+          esq = zeros(size(y))
+      end
+      Z   = autoreg_matrix (esq, p);
+      h   = Z * a;
+      f   = esq ./ h - ones (T,1);
+      Z_tilde = Z ./ (h * ones (1, p+1));
+      delta_a = inv (Z_tilde' * Z_tilde) * Z_tilde' * f;
+      a = a + gamma * delta_a;
+    end
+  endfunction
diff --git a/macros/arch_test.sci b/macros/arch_test.sci
index 3c53fc5..5cb8ed8 100644
--- a/macros/arch_test.sci
+++ b/macros/arch_test.sci
@@ -1,46 +1,58 @@
-function [PVAL, LM]= arch_test(Y,X,P)
-// perform a Lagrange Multiplier (LM) test of thenull hypothesis of no conditional heteroscedascity against the alternative of CH(P)
-//Calling Sequence
-//arch_test(Y,X,P)
-//PVAL = arch_test(Y,X,P)
-//[PVAL, LM]= arch_test(Y,X,P)
-//Parameters 
-//P: Degrees of freedom
-//PVAL:PVAL is the p-value (1 minus the CDF of this distribution at LM) of the test
-//Description
-//perform a Lagrange Multiplier (LM) test of thenull hypothesis of no conditional heteroscedascity against the alternative of CH(P).
-//
-//I.e., the model is
-//
-//          y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t),
-//
-//given Y up to t-1 and X up to t, e(t) is N(0, h(t)) with
-//
-//          h(t) = v + a(1) * e(t-1)^2 + ... + a(p) *e(t-p)^2, and the null is a(1) == ... == a(p) == 0.
-//
-//If the second argument is a scalar integer, k,perform the sametest in a linear autoregression model of orderk, i.e., with
-//
-//          [1, y(t-1), ..., y(t-K)] as the t-th row of X.
-//
-// Under the null, LM approximatel has a chisquare distribution with P degrees of freedom and PVAL is the p-value (1 minus the CDF of this distribution at LM) of the test.
-//
-// If no output argument is given, the p-value is displayed.
-	funcprot(0)
-	rhs= argn(2);
-	lhs= argn(1);
-	if(rhs<3 | rhs>3)
-		error("Wrong number of input arguments");
-	end
-	if(lhs<1 | lhs>2)
-		error("Wrong number of output arguments");
-	end
-	select(rhs)
-	case 3 then
-		select(lhs)
-		case 1 then
-		PVAL= callOctave("arch_test", Y, X, P);
-		case 2 then
-		[PVAL,LM]= callOctave("arch_test", Y, X, P);
-		end
-	end
-endfunction
-\ No newline at end of file
+/*
+Description:
+            Perform a Lagrange Multiplier (LM) test for conditional heteroscedasticity.
+            For a linear regression model
+            y = x * b + e
+            perform a Lagrange Multiplier (LM) test of the null hypothesis of no conditional heteroscedascity against the alternative of CH(p).
+            I.e., the model is
+            y(t) = b(1) * x(t,1) + … + b(k) * x(t,k) + e(t),
+            given y up to t-1 and x up to t, e(t) is N(0, h(t)) with
+            h(t) = v + a(1) * e(t-1)^2 + … + a(p) * e(t-p)^2,
+            and the null is a(1) == … == a(p) == 0.
+            If the second argument is a scalar integer, k, perform the same test in a linear autoregression model of order k, i.e., with
+            [1, y(t-1), …, y(t-k)]
+            as the t-th row of x.
+            Under the null, LM approximately has a chisquare distribution with p degrees of freedom and pval is the p-value (1 minus the CDF of this distribution at LM) of the test.
+            If no output argument is given, the p-value is displayed.
+        Calling Sequence
+            [pval, lm] = arch_test (y, x, p)
+        Parameters
+            y: Array-like. Dependent variable of the regression model.
+            x: Array-like. Independent variables of the regression model. If x is a scalar integer k, it represents the order of autoregression.
+            p : Integer. Number of lagged squared residuals to include in the heteroscedasticity model.
+        Returns:
+            pval: Float. p-value of the LM test.
+            lm: Float. Lagrange Multiplier test statistic.*/
+        Dependencies : ols, autoreg_matrix
+//helper function
+function cdf = chi2cdf ( X, n)
+    df = resize_matrix ( n , size (X) , "" , n);
+    [cdf,Q] = cdfchi ( "PQ" , X ,df);
+endfunction
+//main function
+function [pval, lm] = arch_test (y, x, p)
+  nargin = argn(2)
+  if (nargin ~= 3)
+    error ("arch_test: 3 input arguments required");
+  end
+  if (~ (isvector (y)))
+    error ("arch_test: Y must be a vector");
+  end
+  T = max(size(y));
+  y = matrix (y, T, 1);
+  [rx, cx] = size (x);
+  if ((rx == 1) && (cx == 1))
+    x = autoreg_matrix (y, x);
+  elseif (~ (rx == T))
+    error ("arch_test: either rows (X) == length (Y), or X is a scalar");
+  end
+  if (~ (isscalar (p) && (modulo (p, 1) == 0) && (p > 0)))
+    error ("arch_test: P must be a positive integer");
+  end
+  [b, v_b, e] = ols (y, x);
+  Z    = autoreg_matrix (e.^2, p);
+  f    = e.^2 / v_b - ones (T, 1);
+  f    = Z' * f;
+  lm   = f' * inv (Z'*Z) * f / 2;
+  pval = 1 - chi2cdf (lm, p);
+endfunction
diff --git a/macros/cohere.sci b/macros/cohere.sci
index ba013ab..12a35bc 100644
--- a/macros/cohere.sci
+++ b/macros/cohere.sci
@@ -1,28 +1,38 @@
-function [Pxx,freqs] = cohere(x,y,Nfft,Fs,win,overlap,ran,plot_type,detrends)
-//Estimate (mean square) coherence of signals "x" and "y"
-//Calling Sequence
-// [Pxx,freqs] = cohere(x,y,Nfft,Fs,win,overlap,ran,plot_type,detrends)
-//Parameters 
-//x: [non-empty vector] system-input time-series data
-//y: [non-empty vector] system-output time-series data
-//win:[real vector] of window-function values between 0 and 1; the data segment has the same length as the window. Default window shape is Hamming. [integer scalar] length of each data segment.  The default value is window=sqrt(length(x)) rounded up to the nearest integer power of 2; see 'sloppy' argument.
-//overlap:[real scalar] segment overlap expressed as a multiple of window or segment length.   0 <= overlap < 1, The default is overlap=0.5 .
-//Nfft:[integer scalar] Length of FFT.  The default is the length of the "window" vector or has the same value as the scalar "window" argument.  If Nfft is larger than the segment length, "seg_len", the data segment is padded with "Nfft-seg_len" zeros.  The default is no padding. Nfft values smaller than the length of the data segment (or window) are ignored silently.
-//Fs:[real scalar] sampling frequency (Hertz); default=1.0 
-//range:'half',  'onesided' : frequency range of the spectrum is zero up to but not including Fs/2.  Power from negative frequencies is added to the positive side of the spectrum, but not at zero or Nyquist (Fs/2) frequencies.  This keeps power equal in time and spectral domains.  See reference [2]. 'whole', 'twosided' : frequency range of the spectrum is-Fs/2 to Fs/2, with negative frequenciesstored in "wrap around" order after the positivefrequencies; e.g. frequencies for a 10-point 'twosided'spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1 'shift', 'centerdc' : same as 'whole' but with the first half of the spectrum swapped with second half to put the zero-frequency value in the middle. (See "help fftshift". If data (x and y) are real, the default range is 'half', otherwise default range is 'whole'.
-//plot_type: 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db': specifies the type of plot.  The default is 'plot', which means linear-linear axes. 'squared' is the same as 'plot'. 'dB' plots "10*log10(psd)".  This argument is ignored and a spectrum is not plotted if the caller requires a returned value.
-//detrends:'no-strip', 'none' -- do NOT remove mean value from the data'short', 'mean' -- remove the mean value of each segment from each segment of the data. 'linear',-- remove linear trend from each segment of the data.'long-mean'-- remove the mean value from the data before splitting it into segments. This is the default.
-//Description
-//Estimate (mean square) coherence of signals "x" and "y".
-//
-//Use the Welch (1967) periodogram/FFT method.
-	rhs= argn(2);
-	lhs= argn(1);
-	if(rhs < 10 | rhs > 10)
-		error("Wrong number of input arguments");
-	end
-	select(rhs)
-	case 10 then
-		[Pxx,freqs] = callOctave("cohere",x,y,Nfft,Fs,win,overlap,ran,plot_type,detrends);
-	end
-endfunction
-\ No newline at end of file
+/*
+Calling Sequence
+      [Pxx, freq] = cohere(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
+Estimate (mean square) coherence of signals "x" and "y".
+Use the Welch (1967) periodogram/FFT method.
+Compatible with Matlab R11 cohere and earlier.
+See "help pwelch" for description of arguments, hints and references — especially hint (7) for Matlab R11 defaults. */
+function varargout = cohere(varargin)
+      if ( nargin<2 )
+        error( 'cohere: Need at least 2 args. Use help cohere.' );
+      end
+      nvarargin = length(varargin);
+      // remove any pwelch RESULT args and add 'trans'
+      for iarg=1:nvarargin
+        arg = varargin(iarg);
+        if ( ~isempty(arg) && (type(arg)== 10) && ( ~strcmp(arg,'power') || ...
+               ~strcmp(arg,'cross') || ~strcmp(arg,'trans') || ...
+               ~strcmp(arg,'coher') || ~strcmp(arg,'ypower') ))
+          varargin(iarg) = [];
+        end
+      end
+      varargin(nvarargin+1) = 'coher';
+      saved_compatib = pwelch('R11-');
+      if ( nargout==0 )
+        pwelch(varargin(:));
+      elseif ( nargout==1 )
+        Pxx = pwelch(varargin(:));
+        varargout(1) = Pxx;
+      elseif ( nargout>=2 )
+        [Pxx,f] = pwelch(varargin(:));
+        varargout(1) = Pxx;
+        varargout(2) = f;
+      end
+      pwelch(saved_compatib);
+      saved_compatib = 0;
+    endfunction
+
+   
diff --git a/macros/cplxpair.sci b/macros/cplxpair.sci
new file mode 100644
index 0000000..16fceac
--- /dev/null
+++ b/macros/cplxpair.sci
@@ -0,0 +1,107 @@
+/*Description
+        Sort the numbers z into complex conjugate pairs ordered by increasing real part.
+        The negative imaginary complex numbers are placed first within each pair. All real numbers (those with abs (imag (z) / z) < tol) are placed after the complex pairs.
+        and the resulting tolerance for a given complex pair is 100 * eps (abs (z(i))).
+        Signal an error if some complex numbers could not be paired. Signal an error if all complex numbers are not exact conjugates (to within tol).
+        Note that there is no defined order for pairs with identical real parts but differing imaginary parts
+Calling Sequence
+         cplxpair (z)
+         cplxpair (z, tol)
+         cplxpair (z, tol, dim)
+Parameters:
+         z is a matrix or vector.
+         tol is a weighting factor which determines the tolerance of matching. The default value is 100.
+         By default the complex pairs are sorted along the first non-singleton dimension of z. If dim is specified, then the complex pairs are sorted along this dimension.
+Dependencies: ipermute
+Example:
+        The following code
+        [ cplxpair(exp(2*%i*%pi*[0:4]'/5)), exp(2*%i*%pi*[3; 2; 4; 1; 0]/5) ]
+        Produces the following output
+                ans =
+             -0.80902 - 0.58779i  -0.80902 - 0.58779i
+             -0.80902 + 0.58779i  -0.80902 + 0.58779i
+              0.30902 - 0.95106i   0.30902 - 0.95106i
+              0.30902 + 0.95106i   0.30902 + 0.95106i
+              1.00000 + 0.00000i   1.00000 + 0.00000i
+*/
+function zsort = cplxpair (z, tol, dim)
+  if (nargin < 1)
+    error("Invalid inputs");
+  end
+  // default  double
+  realmin = 2.2251e-308
+  if (isempty (z))
+    zsort = zeros (size (z,1) , size (z,2));
+    return;
+  end
+  if (nargin < 2 || isempty (tol))
+    tol = 100* %eps;
+  elseif (~ isscalar (tol) || tol < 0 || tol >= 1)
+    error ("cplxpair: TOL must be a scalar number in the range 0 <= TOL < 1");
+  end
+  nd = ndims (z);
+  if (nargin < 3)
+    // Find the first singleton dimension.
+    sz = size (z);
+    dim = find (sz > 1, 1);
+    if isempty(dim) 
+        dim  = 1;
+    end
+  else
+    dim = floor (dim);
+    if (dim < 1 || dim > nd)
+      error ("cplxpair: invalid dimension DIM");
+    end
+  end
+  // Move dimension to analyze to first position, and convert to a 2-D matrix.
+  perm = [dim:nd, 1:dim-1];
+  z = permute (z, perm);
+  sz = size (z);
+  n = sz(1);
+  m = prod (sz) / n;
+  z = matrix (z, n, m);
+  // Sort the sequence in terms of increasing real values.
+  [temp, idx] = gsort (real (z), 1 , "i");
+  z = z(idx + n * ones (n, 1) * [0:m-1]);
+  // Put the purely real values at the end of the returned list.
+  [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin) <= tol);
+  // Force values detected to be real within tolerance to actually be real.
+  z(idxi + n*(idxj-1)) = real (z(idxi + n*(idxj-1)));
+  //if idxi and idxj are not scalers
+  if ~isscalar(idxi) then
+      v = ones(size(idxi,1),size(idxi,2));
+  else
+      v = 1 ;
+  end
+  q = sparse ([idxi' idxj'], v, [n m]);
+  nr = sum (q, 1);
+  [temp, idx] = gsort (q, 'r','i');
+  midx = idx + size (idx,1) * ones (size (idx,1), 1) * [0:size(idx,2)-1];
+  z = z(midx);
+  zsort = z;
+  // For each remaining z, place the value and its conjugate at the start of
+  // the returned list, and remove them from further consideration.
+  for j = 1:m
+    p = n - nr(j);
+    for i = 1:2:p
+      if (i+1 > p)
+        error ("cplxpair: could not pair all complex numbers");
+      end
+      [v, idx] = min (abs (z(i+1:p,j) - conj (z(i,j))));
+      if (v >= tol * abs (z(i,j)))
+        error ("cplxpair: could not pair all complex numbers");
+      end
+      // For pairs, select the one with positive imaginary part and use it and
+      // it's conjugate, but list the negative imaginary pair first.
+      if (imag (z(i,j)) > 0)
+        zsort([i, i+1],j) = [conj(z(i,j)), z(i,j)];
+      else
+        zsort([i, i+1],j) = [conj(z(idx+i,j)), z(idx+i,j)];
+      end
+      z(idx+i,j) = z(i+1,j);
+    end
+  end
+  // Reshape the output matrix.
+  zsort = ipermute (matrix (zsort, sz), perm);
+endfunction
+
diff --git a/macros/cplxreal.sci b/macros/cplxreal.sci
index a552998..b06fac3 100644
--- a/macros/cplxreal.sci
+++ b/macros/cplxreal.sci
@@ -1,48 +1,51 @@
-function [zc, zr] = cplxreal (z, thresh)
-//Function to divide vector z into complex and real elements, removing the one of each complex conjugate pair.
-//Calling Sequence
-//[zc, zr] = cplxreal (z, thresh)
-//[zc, zr] = cplxreal (z)
-//zc = cplxreal (z, thresh)
-//zc = cplxreal (z)
-//Parameters 
-//z: vector of complex numbers.
-//thresh: tolerance for comparisons.
-//zc: vector containing the elements of z that have positive imaginary parts.
-//zr: vector containing the elements of z that are real.
-//Description
-//This is an Octave function.
-//Every complex element of z is expected to have a complex-conjugate elsewhere in z. From the pair of complex-conjugates, the one with the negative imaginary part is removed.
-//If the magnitude of the imaginary part of an element is less than the thresh, it is declared as real.  
-//Examples
-//[zc, zr] = cplxreal([1 2 3+i 4 3-i 5])
-//zc =  3 + 1i
-//zr =
-//   1   2   4   5
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 1 | rhs > 2)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-	case 1 then
-		if(lhs==1)
-		zc = callOctave("cplxreal",z)
-		elseif (lhs==2)
-		[zc, zr] = callOctave("cplxreal",z)
-		else
-		error("Wrong number of output argments.")
-		end
-
-	case 2 then
-		if(lhs==1)
-		zc = callOctave("cplxreal",z, thresh)
-		elseif (lhs==2)
-		[zc, zr] = callOctave("cplxreal",z, thresh)
-		else
-		error("Wrong number of output argments.")
-		end
-	end
+/*Description
+        Sort the numbers z into complex-conjugate-valued and real-valued elements.
+        The positive imaginary complex numbers of each complex conjugate pair are returned in zc and the real numbers are returned in zr.
+        Signal an error if some complex numbers could not be paired.
+        Signal an error if all complex numbers are not exact conjugates (to within tol).
+        Note that there is no defined order for pairs with identical real parts but differing imaginary parts.
+Calling Sequence:
+        [zc, zr] = cplxreal (z)
+        [zc, zr] = cplxreal (z, tol)
+        [zc, zr] = cplxreal (z, tol, dim)
+Parameters
+        Inputs
+             z  -  A vector of numbers or Matrix
+            tol - tol is a weighting factor in the range [0, 1) which determines the tolerance of the matching.
+            The default value is 100 * eps and the resulting tolerance for a given complex pair is tol * abs (z(i))).
+            dim - By default the complex pairs are sorted along the first non-singleton dimension of z. If dim is specified, then the complex pairs are sorted along this dimension.
+        Outputs
+            zc - complex conjugate pair
+            zr -  real numbers
+Example:
+        with 2 real zeros, one of them equal to 0
+        [zc, zr] = cplxreal (roots ([1, 0, 0, 1, 0])) */
+function [zc, zr] = cplxreal (z, tol, dim)
+  if (nargin < 1 || nargin > 3)
+    error("invalid inputs");
+  end
+  if (isempty (z))
+    zc = zeros (size (z,1),size(z,2));
+    zr = zeros (size (z,1),size(z,2));
+    return;
+  end
+  if (nargin < 2 || isempty (tol))
+    tol = 100 * %eps ;
+  end
+  if (nargin >= 3)
+    zcp = cplxpair(z,tol,dim);
+  else
+    zcp = cplxpair (z , tol);
+  end
+  nz = max(size (z) );
+  idx = nz;
+  while ((idx > 0) && (zcp(idx) == 0 || (abs (imag (zcp(idx))) ./ abs (zcp(idx))) <= tol))
+    zcp(idx) = real (zcp(idx));
+    idx = idx - 1;
+  end
+  if (pmodulo (idx, 2) ~= 0)
+    error ("cplxreal: odd number of complex values was returned from cplxpair");
+  end
+  zc = zcp(2:2:idx);
+  zr = zcp(idx+1:nz);
 endfunction
diff --git a/macros/cpsd.sci b/macros/cpsd.sci
index 5e56d7a..0d8fcce 100644
--- a/macros/cpsd.sci
+++ b/macros/cpsd.sci
@@ -1,50 +1,39 @@
-function [PXX, FREQ] = cpsd(X, Y, varargin) 
-//This function estimates cross power spectrum of data x and y by the Welch (1967) periodogram/FFT method.
-//Calling Sequence
-//[PXX, FREQ] = cpsd(X, Y)
-//[...] = cpsd(X, Y, WINDOW)
-//[...] = cpsd(X, Y, WINDOW, OVERLAP)
-//[...] = cpsd(X, Y, WINDOW, OVERLAP, NFFT)
-//[...] = cpsd(X, Y, WINDOW, OVERLAP, NFFT, FS)
-//[...] = cpsd(X, Y, WINDOW, OVERLAP, NFFT, FS, RANGE)
-//cpsd(...)
-//Parameters
-//X, Y: Matrix or integer
-//Description
-//Estimate cross power spectrum of data X and Y by the Welch (1967) periodogram/FFT method.
-//Examples
-// [a, b] = cpsd([1,2,3],[4,5,6])
-//ans = 
-//     b  = 
-//          0.    
-//          0.25  
-//          0.5   
-//     a  = 
-//          2.7804939               
-//          4.4785583 + 1.0743784i  
-//          0.7729851               
-funcprot(0);
-rhs=argn(2);
-lhs=argn(1);
-if(rhs<2 | rhs>7) then
-    error("Wrong number of input arguments.");
-end
-if (lhs<2 | lhs>2)
-    error("Wrong number of output arguments.");
-end
-select(rhs)
-case 2 then
-    [PXX, FREQ] = callOctave("cpsd",X, Y);
-case 3 then
-    [PXX, FREQ] = callOctave("cpsd",X, Y, varargin(1));
-case 4 then 
-    [PXX, FREQ] = callOctave("cpsd",X, Y, varargin(1), varargin(2));
-case 5 then
-    [PXX, FREQ] = callOctave("cpsd",X, Y, varargin(1), varargin(2), varargin(3));
-case 6 then
-    [PXX, FREQ] = callOctave("cpsd",X, Y, varargin(1), varargin(2), varargin(3), varargin(4));
-case 7 then 
-    [PXX, FREQ] = callOctave("cpsd",X, Y, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5)); 
-end
+/*Calling Sequence
+   [Pxx, freq] = cpsd (x, y)
+   […] = cpsd (x, y, window)
+   […] = cpsd (x, y, window, overlap)
+   […] = cpsd (x, y, window, overlap, Nfft)
+   […] = cpsd (x, y, window, overlap, Nfft, Fs)
+   […] = cpsd (x, y, window, overlap, Nfft, Fs, range)
+   cpsd (…)
+Estimate cross power spectrum of data x and y by the Welch (1967) periodogram/FFT method.
+See "help pwelch" for description of arguments, hints and references
+*/
+function varargout = cpsd(varargin)
+    // Check fixed argument
+    if (nargin < 2 || nargin > 7)
+      error( "Invalid number of inputs" );
+    end
+    nvarargin = length(varargin);
+    // remove any pwelch RESULT args and add 'cross'
+    for iarg=1:nvarargin
+      arg = varargin(iarg);
+      if ( ~isempty(arg) && (type(arg) == 10 ) && ( ~strcmp(arg,'power') || ...
+             ~strcmp(arg,'cross') || ~strcmp(arg,'trans') || ...
+             ~strcmp(arg,'coher') || ~strcmp(arg,'ypower') ))
+        varargin(iarg) = [];
+      end
+    end
+    varargin(nvarargin+1) = 'cross';
+    if ( nargout==0 )
+      pwelch(varargin(:));
+    elseif ( nargout==1 )
+      Pxx = pwelch(varargin(:));
+      varargout(1) = Pxx;
+    elseif ( nargout>=2 )
+      [Pxx,f] = pwelch(varargin(:));
+      varargout(1) = Pxx;
+      varargout(2) = f;
+    end
+  endfunction
 
-endfunction
diff --git a/macros/czt.sci b/macros/czt.sci
index 84a0253..f7c4138 100644
--- a/macros/czt.sci
+++ b/macros/czt.sci
@@ -1,41 +1,65 @@
-function y = czt(x, varargin)
-//Chirp Z Transform
-//Calling Sequence
-//czt (x)
-//czt (x, m)
-//czt (x, m, w)
-//czt (x, m, w, a)
-//Parameters 
-//x: Input scalar or vector
-//m: Total Number of steps
-//w: ratio between points in each step
-//a: point in the complex plane 
-//Description
-//This is an Octave function.
-//Chirp z-transform. Compute the frequency response starting at a and stepping by w for m steps. a is a point in the complex plane, and w is the ratio between points in each step (i.e., radius increases exponentially, and angle increases linearly).
-//Examples
-// m = 32;                               ## number of points desired
-// w = exp(-j*2*pi*(f2-f1)/((m-1)*Fs));  ## freq. step of f2-f1/m
-// a = exp(j*2*pi*f1/Fs);                ## starting at frequency f1
-// y = czt(x, m, w, a);
-
-funcprot(0);
-lhs= argn(1);
-rhs= argn(2);
-
-if(rhs<1 | rhs > 4)
-error("Wrong number of input arguments")
-end
-
-select (rhs)
-	case 1 then
-		y= callOctave("czt", x);
-	case 2 then
-		y= callOctave("czt", x, varargin(1));
-	case 3 then
-		y= callOctave("czt", x, varargin(1), varargin(2));
-	case 4 then
-		y= callOctave("czt", x, varargin(1), varargin(2), varargin(3));
-end
+/*Description
+        Chirp z-transform. Compute the frequency response starting at a and stepping by w for m steps. a is a point in the complex plane,
+        and w is the ratio between points in each step (i.e., radius increases exponentially, and angle increases linearly).
+  Calling Sequence
+        czt (x)
+        czt (x, m)
+        czt (x, m, w)
+        czt (x, m, w, a)
+  Parameters
+        x: Input scalar or vector
+        m: Total Number of steps
+        w: ratio between points in each step
+        a: point in the complex plane
+  Examples: This example uses the czt function to determine the frequency components of a signal, as shown in the following
+    t=linspace(0,50,1000); 
+    f=linspace(0,3,1000);    
+    x_t=sin(t) + cos(t*2*%pi);  
+    x_f=czt(x_t);   
+    plot(f,abs(x_f)); */
+function y = czt(x, m, w, a)
+    funcprot(0);
+    nargin=argn(2);
+    if nargin < 1 || nargin > 4 then
+        error("Please input valid number of arguments");
+    end
+    [row, col] = size(x);
+    if row == 1 then
+        x = x(:); col = 1;
+    end
+    if nargin < 2 || isempty(m) then
+        m = max(size(x(:,1)));
+    end
+    if max(size(m) ) > 1 then
+        error("czt: m must be a single element\n");
+    end
+    if nargin < 3 || isempty(w) then
+        w = exp(-2*%i*%pi/m);
+    end
+    if nargin < 4 || isempty(a) then
+        a = 1;
+    end
+    if max(size(w)) > 1 then
+        error("czt: w must be a single element\n");
+    end
+    if max(size(a)) > 1 then
+        error("czt: a must be a single element\n");
+    end
+    // indexing to make the statements a little more compact
+    n = max(size(x(:,1)));
+    N = [0:n-1]'+n;
+    NM = [-(n-1):(m-1)]'+n;
+    M = [0:m-1]'+n;
+    nfft = 2^nextpow2(n+m-1); // fft pad
+    W2 = w.^(([-(n-1):max(m-1,n-1)]'.^2)/2); // chirp
+    for idx = 1:col
+        fg = fft1(x(:,idx).*(a.^-(N-n)).*W2(N), nfft);
+        fw = fft1(1./W2(NM), nfft);
+        gg = ifft1(fg.*fw, nfft);
+        y(:,idx) = gg(M).*W2(M);
+    end
+    if row == 1, y = y.';
+    end
+    y = clean ( y ) ;
 endfunction
 
diff --git a/macros/fft1.sci b/macros/fft1.sci
index a0a6438..f8a27fc 100644
--- a/macros/fft1.sci
+++ b/macros/fft1.sci
@@ -1,49 +1,63 @@
-function res =  fft1 (x, n, dim)
-//Calculates the discrete Fourier transform of a matrix using Fast Fourier Transform algorithm.
-//Calling Sequence
-//fft (x, n, dim)
-//fft (x, n)
-//fft (x)
-//Parameters 
-//x: input matrix
-//n: Specifies the number of elements of x to be used
-//dim: Specifies the dimention of the matrix along which the FFT is performed
-//Description
-//This is an Octave function.
-//The FFT is calculated along the first non-singleton dimension of the array. Thus, FFT is computed for each column of x.
-//
-//n is an integer specifying the number of elements of x to use. If n is larger than dimention along. which the FFT is calculated, then x is resized and padded with zeros.
-//Similarly, if n is smaller, then x is truncated.
-//
-//dim is an integer specifying the dimension of the matrix along which the FFT is performed.
-//Examples
-//x = [1 2 3; 4 5 6; 7 8 9]
-//n = 3
-//dim = 2
-//fft1 (x, n, dim)
-//ans =
-//
-//    6.0000 +  0.0000i   -1.5000 +  0.8660i   -1.5000 -  0.8660i
-//   15.0000 +  0.0000i   -1.5000 +  0.8660i   -1.5000 -  0.8660i
-//   24.0000 +  0.0000i   -1.5000 +  0.8660i   -1.5000 -  0.8660i
-
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 1 | rhs > 3)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-	
-	case 1 then
-	res = callOctave("fft", x)
-
-	case 2 then
-	res = callOctave("fft", x, n)
-
-	case 3 then
-	res = callOctave("fft", x, n, dim)
-
-	end
-endfunction
+/*Description
+        Calculates the discrete Fourier transform of a matrix using Fast Fourier Transform algorithm.
+        The FFT is calculated along the first non-singleton dimension of the array. Thus, FFT is computed for each column of D.
+        The variable 'N' is an integer that determines the number of elements of 'D' to use.
+        If 'N' is larger than the dimension along which the FFT is calculated,
+        then 'D' is resized and padded with zeros to match the required size.On the other hand,
+        if 'N' is smaller than the size of 'D', then 'D' is truncated to match the required size.
+        DIM is an integer specifying the dimension of the matrix along which the FFT is performed.
+Calling Sequence
+       fft1 (D)
+       fft1 (D, N)
+       fft1 (D, N, DIM)
+Parameters
+        D: input matrix
+        N: Specifies the number of elements of x to be used
+        DIM: Specifies the dimention of the matrix along which the FFT is performed
+Examples
+        D = [1 2 3; 4 5 6; 7 8 9]
+        N = 3
+        DIM = 2
+        fft1 (D,N,DIM)
+        ans =
+            6.0000 +  0.0000i   -1.5000 +  0.8660i   -1.5000 -  0.8660i
+            15.0000 +  0.0000i   -1.5000 +  0.8660i   -1.5000 -  0.8660i
+            24.0000 +  0.0000i   -1.5000 +  0.8660i   -1.5000 -  0.8660i */
+function res =  fft1 (D, N, DIM)
+    funcprot(0);
+    lhs = argn(1)
+    rhs = argn(2)
+    if (rhs < 1 | rhs > 3)
+        error("Wrong number of input arguments.")
+    end
+    dimension = size(D);
+    // first non-singleton dimension
+    nsdim = find(dimension >1,1)
+    if isempty(nsdim) then
+        nsdim = 1 // default to 1 to avoid error while calling fft
+    end    
+    select(rhs)
+    case 1 then
+        res = fft(D, -1, nsdim);
+    case 2 then
+        if isempty(N) then
+            n = size(D, nsdim);
+        else
+            n = N;
+        end
+        new_size = size(D);
+        new_size(nsdim) = n;
+        D = resize_matrix(D, new_size);
+        res = fft(D, -1, nsdim);
+    case 3 then
+        if isempty(N) then
+            n = size(D, DIM);
+        else
+            n = N;
+        end
+        new_size = size(D);
+        new_size(DIM) = n;
+        D = resize_matrix(D, new_size);
+        res= fft(D, -1, DIM);
+        end
+endfunction
diff --git a/macros/fft21.sci b/macros/fft21.sci
index 2f04200..ee73a6b 100644
--- a/macros/fft21.sci
+++ b/macros/fft21.sci
@@ -1,45 +1,57 @@
-function res = fft21 (A, m, n)
-//Calculates the two-dimensional discrete Fourier transform of A using a Fast Fourier Transform algorithm.
-//Calling Sequence
-//fft2 (A, m, n)
-//fft2 (A)
-//Parameters 
-//A: input matrix
-//m: number of rows of A to be used
-//n: number of columns of A to be used
-//Description
-//This is an Octave function.
-//It performs two-dimentional FFT on the matrix A. m and n may be used specify the number of rows and columns of A to use. If either of these is larger than the size of A, A is resized and padded with zeros.
-//If A is a multi-dimensional matrix, each two-dimensional sub-matrix of A is treated separately.
-//Examples
-//x = [1 2 3; 4 5 6; 7 8 9]
-//m = 4
-//n = 4
-//fft21 (A, m, n)
-//ans =
-//
-//   45 +  0i   -6 - 15i   15 +  0i   -6 + 15i
-//  -18 - 15i   -5 +  8i   -6 -  5i    5 -  4i
-//   15 +  0i   -2 -  5i    5 +  0i   -2 +  5i
-//  -18 + 15i    5 +  4i   -6 +  5i   -5 -  8i
-
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 1 | rhs > 3)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-	
-	case 1 then
-	res = callOctave("fft2", A)
-
-	case 2 then
-	error("Wrong number of input arguments.")
-
-	case 3 then
-	res = callOctave("fft2", A, m, n)
-
-	end
-endfunction
+/* Description
+        Calculates the two-dimensional discrete Fourier transform of A using a Fast Fourier Transform algorithm.
+        It performs two-dimentional FFT on the matrix A. You can use the variables m and n to specify the number of rows and columns
+        of A that you want to use. If either of these variables is larger than the size of A,
+        then A will be resized, and zeros will be added as padding.
+        If A is a multi-dimensional matrix, the function will treat each two-dimensional sub-matrix of A separately.
+    Calling Sequence
+        fft21 (A)
+        fft21 (A, m, n)
+    Parameters 
+        A: input matrix
+        m: number of rows of A to be used
+        n: number of columns of A to be used
+    Examples
+        A = [1 2 3; 4 5 6; 7 8 9]
+        m = 4
+        n = 4
+        fft21 (A, m, n)
+        ans =
+            45 +  0i   -6 - 15i   15 +  0i   -6 + 15i
+            -18 - 15i   -5 +  8i   -6 -  5i    5 -  4i
+            15 +  0i   -2 -  5i    5 +  0i   -2 +  5i
+            -18 + 15i    5 +  4i   -6 +  5i   -5 -  8i */
+function res = fft21 (A, m, n)
+    funcprot(0);
+    lhs = argn(1)
+    rhs = argn(2)
+    if (rhs < 1 | rhs > 3)
+        error("Wrong number of input arguments.")
+    end
+    siz=size(A)
+    len=length(siz)
+    select(rhs)
+    case 1 then
+        if len>2 then
+            last_dim=siz(len)
+        else
+            last_dim=1
+        end
+        res=[]
+        for i=1:last_dim
+            res(:,:,i)=fft(A(:,:,i),-1)
+        end
+    case 2 then
+        error("Wrong number of input arguments.")
+    case 3 then
+        if len>2 then
+            last_dim=siz(len)
+        else
+            last_dim=1
+        end
+        res=[]
+        for i=1:last_dim
+            res(:,:,i)=fft(resize_matrix(A(:,:,i),m,n),-1)
+        end;
+    end
+endfunction
diff --git a/macros/fftconv.sci b/macros/fftconv.sci
index d39441d..5028359 100644
--- a/macros/fftconv.sci
+++ b/macros/fftconv.sci
@@ -1,26 +1,48 @@
-function y = fftconv(X, Y, varargin)
-//Convolve two vectors using the FFT for computation.
-//Calling Sequence
-//Y = fftconv(X, Y)
-//Y = fftconv(X, Y, N)
-//Parameters
-//X, Y: Vectors 
-//Description
-//Convolve two vectors using the FFT for computation. 'c' = fftconv (X, Y)' returns a vector of length equal to 'length(X) + length (Y) - 1'.  If X and Y are the coefficient vectors of two polynomials, the returned value is the coefficient vector of the product polynomial.
-//Examples
-//fftconv([1,2,3], [3,4,5])
-//ans = 
-//    3.    10.    22.    22.    15. 
-funcprot(0);
-rhs = argn(2)
-if(rhs<2 | rhs>3)
-error("Wrong number of input arguments.");
-end
-
-	select(rhs)
-	case 2 then
-	y = callOctave("fftconv", X, Y);
-	case 3 then
-	y = callOctave("fftconv",X, Y, varargin(1));
-	end
-endfunction
+/*Description
+        Convolve two vectors using the FFT for computation. c = fftconv (X, Y) returns a vector of length equal to 'length(X) + length (Y) - 1'.  If X and Y are         the coefficient vectors of two polynomials, the returned value is the coefficient vector of the product polynomial.
+        If the optional argument n is specified, an N-point FFT is used.
+    Calling Sequence
+        Y = fftconv(X, Y)
+        Y = fftconv(X, Y, n)
+    Parameters
+        X, Y: Vectors
+    Examples
+        fftconv([1,2,3], [3,4,5])
+        ans = 
+            3.    10.    22.    22.    15. 
+ */
+ function y = fftconv(X, Y, n)
+     funcprot(0);
+     rhs = argn(2);
+     if(rhs<2 | rhs>3)
+         error("Wrong number of input arguments.");
+     end
+     shape_x = size(X);
+     shape_y = size(Y);
+     if (shape_x(1) ~= 1 || length(shape_x) ~= 2 || shape_y(1) ~= 1 || length(shape_y) ~= 2)
+         error('The inputs must be a vector');
+     end
+     nx=length(X);
+     ny=length(Y);
+     select(rhs)
+     case 2 then
+         n=nx + ny;
+         X=resize_matrix(X,1,n);
+         Y=resize_matrix(Y,1,n);
+         fftX=fft(X);
+         fftY=fft(Y);
+         y=fft(fftX.*fftY,1);
+         y=y(1:nx+ny-1);
+     case 3 then
+         n = 2^(fix(log2(nx+ny))+1);
+         X=resize_matrix(X,1,n);
+         Y=resize_matrix(Y,1,n);
+         fftX=fft(X);
+         fftY=fft(Y);
+         y=fft(fftX.*fftY,1);
+         y=y(1:nx+ny-1);
+     end
+     y=clean(y);
+ endfunction
+  
+
diff --git a/macros/fftn.sci b/macros/fftn.sci
index b8b4170..c90b096 100644
--- a/macros/fftn.sci
+++ b/macros/fftn.sci
@@ -1,26 +1,57 @@
+/*Description:
+        This function computes the N-dimensional discrete Fourier transform of A using a Fast Fourier Transform (FFT) algorithm.
+        The optional vector argument SIZE may be used to specify the dimensions of the array to be used.
+        If an element of SIZE is smaller than the corresponding dimension of A, then the dimension of A is truncated prior to performing the FFT.
+        Otherwise, if an element of SIZE is larger than the corresponding dimension, then A is resized and padded with zeros.
+  Calling sequence:
+         Y = fftn(A)
+         Y = fftn(A, SIZE)
+  Parameters:
+        A: Matrix, the input data for which the FFT is computed.
+        SIZE: Optional vector specifying the dimensions of the output array. If provided, the dimensions of A are adjusted accordingly.
+  Examples:
+        fftn([6 9 7 ;2 9 9 ;0 3 1],[2 2])
+        ans  =
+        26.  -10.
+        4.    4. 
+ */
 function y = fftn(A, SIZE)
-//This function computes the N-dimensional discrete Fourier transform of A using a Fast Fourier Transform (FFT) algorithm.
-//Calling Sequence
-//Y = fftn(A)
-//Y = fftn(A, size)
-//Parameters
-//A: Matrix 
-//Description
-//This function computes the N-dimensional discrete Fourier transform of A using a Fast Fourier Transform (FFT) algorithm. The optional vector argument SIZE may be used specify the dimensions of the array to be used.  If an element of SIZE is smaller than the corresponding dimension of A, then the dimension of A is truncated prior to performing the FFT.  Otherwise, if an element of SIZE is larger than the corresponding dimension then A is resized and padded with zeros.
-//Examples
-//fftn([2,3,4])
-//ans = 
-//     9.  - 1.5 + 0.8660254i  - 1.5 - 0.8660254i  
-funcprot(0);
-rhs = argn(2)
-if(rhs<1 | rhs>2)
-error("Wrong number of input arguments.");
-end
-
-	select(rhs)
-	case 1 then
-	y = callOctave("fftn",A);
-	case 2 then
-	y = callOctave("fftn",A, SIZE);
-	end
+    funcprot(0);
+    // Get the number of input arguments
+    rhs = argn(2);
+    // Check if the number of input arguments is valid
+    if rhs < 1 | rhs > 2 then
+        error("Wrong number of input arguments.");
+    end
+    // Perform action based on the number of input arguments
+    switch rhs
+    case 1 // If only A is provided
+        y = fft(A);
+    case 2 // If both A and SIZE are provided
+        // Check if A needs resizing
+        if size(A) == SIZE then
+            // No resizing needed
+            break;
+        elseif length(size(A)) ~= length(SIZE) then
+            error("Output size must have at least Ndims");
+        else
+            // Resize A using the resize_matrix function
+            A = resize_matrix(A, SIZE);
+        end
+        y = fft(A);
+    end
 endfunction
+/*
+//test cases
+i=%i;tol = 0.00001; pi = %pi ;
+//test 1
+Y = [41.0000 + 0*i ,-2.5000 +  2.5981*i ,-2.5000 -  2.5981*i ;2.0000 +  8.6603*i  ,  2.0000 +  3.4641*i  , -8.5000 + 11.2583*i ;2.0000 -  8.6603*i  , -8.5000 - 11.2583*i ,   2.0000 -  3.4641*i]
+assert_checkalmostequal(fftn([3 7 5;0 1 7;9 5 4]),Y,tol)
+// test
+assert_checkalmostequal(fftn([100 278;334 985]),[1697,-829;-941,473],tol)
+assert_checkalmostequal(fftn([2:5;8 7 5 1;7 3 4 5;0 0 1 6],[2 2]),[20,0;-10,-2],tol)
+assert_checkalmostequal(fftn([[i 2+3*i 2+pi;pi i+2 pi+i;1 5+2*5*i pi+10*i]],[2,2]),[7.1416 + 5.0000*i , -0.8584 - 3.0000*i; -3.1416 + 3.0000*i , -3.1416 - 1.0000*i],tol)
+//error
+fftn([2:5;8 7 5 1;7 3 4 5;0 0 1 6],[2]);
+*/
+
diff --git a/macros/fht.sci b/macros/fht.sci
index f7e5c8e..9bd6b38 100644
--- a/macros/fht.sci
+++ b/macros/fht.sci
@@ -1,29 +1,61 @@
-function y=fht(d,n,dim)
-//The Function calculates the Fast Hartley Transform of real input.
-//Calling Sequence
-//M = fht (D)
-//M = fht (D, N)
-//M = fht (D, N, DIM)
-//Parameters 
-//Description
-//This function calculates the Fast Hartley transform of real input D. If D is a matrix, the Hartley transform is calculated along the columns by default.
-//Examples
-//fht(1:4)
-//ans =
-//   10   -4   -2   0  
-//This function is being called from Octave.
-funcprot(0);
-rhs=argn(2);
-if(rhs<1 | rhs>3)
-    error("Wrong number of input arguments.")
-end
-select(rhs)
-case 1 then
-    y=callOctave("fht",d)
-case 2 then
-    y=callOctave("fht",d,n)
-case 3 then
-    y=callOctave("fht",d,n,dim)
-end
-
-endfunction
+/* Description
+        This function calculates the Fast Hartley transform of real input D.
+        If D is a matrix, the Hartley transform is calculated along the columns by default.
+        If N is specified, the first N elements along the specified dimension are used for the transform.
+        If DIM is specified, the transform is calculated along the specified dimension.
+   Calling Sequence
+        M = fht (D)
+        M = fht (D, N)
+        M = fht (D, N, DIM)
+   Parameters 
+        D: Input data (real matrix or vector).
+        N: Number of elements of D to be used for the transform (optional).
+        DIM: Dimension along which the transform is to be computed (optional).
+   Examples
+        fht(1:4)
+        ans =
+       10   -4   -2   0  
+*/
+function M = fht(D, N, DIM)
+    funcprot(0);
+    rhs = argn(2);
+    if rhs < 1 | rhs > 3 then
+        error("Wrong number of input arguments.")
+    end
+    // The fht will be calculated along the first non-singleton dimension of the array i.e along the columns by default.
+    dimension = size(D);
+    nsdim = 1;
+    for i = 1:length(dimension)
+        if dimension(i) ~= 1 then
+            nsdim = i;
+            break;
+        end
+    end
+    // Process input arguments
+    select(rhs)
+    case 1 then
+        M = fft(D, -1, nsdim);
+    case 2 then
+        if isempty(N) then
+            n = size(D, nsdim);
+        else
+            n = N;
+        end
+        new_size = size(D);
+        new_size(nsdim) = n;
+        D = resize_matrix(D, new_size);
+        M = fft(D, -1, nsdim);
+    case 3 then
+        if isempty(N) then
+            n = size(D, DIM);
+        else
+            n = N;
+        end
+        new_size = size(D);
+        new_size(DIM) = n;
+        D = resize_matrix(D, new_size);
+        M = fft(D, -1, DIM);
+    end
+    // Return real part of the result
+    M = real(M) - imag(M);
+endfunction
diff --git a/macros/grpdelay.sci b/macros/grpdelay.sci
index a7ffdad..dec9823 100644
--- a/macros/grpdelay.sci
+++ b/macros/grpdelay.sci
@@ -1,51 +1,157 @@
-function [gd,w] = grpdelay(b,a,nfft,whole,Fs)
-//This function computes the group delay of a filter.
-//Calling Sequence
-//gd = grpdelay(b)
-//gd = grpdelay(b, a)
-//gd = grpdelay(b, a, nfft)
-//gd = grpdelay(b, a, nfft, whole)
-//gd = grpdelay(b, a, nfft, whole, Fs)
-//[gd, w] = grpdelay(...)
+/*Description:
+        Compute the group delay of a filter.
+        Theory: group delay, g(w) = -d/dw [arg{H(e^jw)}], is the rate of change of phase with respect to frequency. It can be computed as:
+                      d/dw H(e^-jw)
+              g(w) = -------------
+                        H(e^-jw)
+        where
+                H(z) = B(z)/A(z) = sum(b_k z^k)/sum(a_k z^k).
+        By the quotient rule,
+                          A(z) d/dw B(z) - B(z) d/dw A(z)
+              d/dw H(z) = -------------------------------
+                                      A(z) A(z)
+        Substituting into the expression above yields:
+                      A dB - B dA
+              g(w) =  ----------- = dB/B - dA/A
+                          A B
+        Note that,
+              d/dw B(e^-jw) = sum(k b_k e^-jwk)
+              d/dw A(e^-jw) = sum(k a_k e^-jwk)
+        which is just the FFT of the coefficients multiplied by a ramp.
+Calling Sequence
+        [g, w] = grpdelay (b) : returns the group delay g of the FIR filter with coefficients b. The response is evaluated at 512 angular frequencies between 0 and pi. w is a vector containing the 512 frequencies. The group delay is in units of samples. It can be converted to seconds by multiplying by the sampling period (or dividing by the sampling rate fs).
+        [g, w] = grpdelay (b, a) : returns the group delay of the rational IIR filter whose numerator has coefficients b and denominator coefficients a.
+        [g, w] = grpdelay (…, n) : returns the group delay evaluated at n angular frequencies. For fastest computation n should factor into a small number of small primes.
+        [g, w] = grpdelay (…, n, "whole") : evaluates the group delay at n frequencies between 0 and 2*pi.
+        [g, f] = grpdelay (…, n, Fs) : evaluates the group delay at n frequencies between 0 and Fs/2.
+        [g, f] = grpdelay (…, n, "whole", Fs) : evaluates the group delay at n frequencies between 0 and Fs.
+        [g, w] = grpdelay (…, w) : evaluates the group delay at frequencies w (radians per sample).
+        [g, f] = grpdelay (…, f, Fs) : evaluates the group delay at frequencies f (in Hz).
+        grpdelay(...) plots the group delay vs. frequency.
+         If the denominator of the computation becomes too small, the group delay is set to zero. (The group delay approaches infinity when there are poles or zeros very close to the unit circle in the z plane.)
+Parameters:
+          b : Coefficients of the numerator polynomial of the filter (FIR or IIR).
+          a : (optional): Coefficients of the denominator polynomial of the filter (only for IIR filters).
+          n : (optional): Number of angular frequencies at which to evaluate the group delay. For efficient computation, n should be a product of small primes.
+          "whole" : (optional): Specifies frequency range. If specified: Frequencies range from 0 to 2π radians per sample.
+          Fs : (numeric value): Frequencies range from 0 to Fs Hz or Fs radians per sample.
+          w : (optional): Vector of specific frequencies (in radians per sample) at which to evaluate the group delay.
+          f : (optional): Vector of specific frequencies (in Hz) at which to evaluate the group delay.
+          Fs : (optional): Sampling rate in Hz. Used to specify the frequency range when f or "whole" is specified.
+Dependencies: fft1
+*/
+function [gd,w] = grpdelay (b, a, n, whole, Fs)
+    lhs=argn(1);
+    rhs= argn(2);
+    if (rhs < 1 | rhs > 5) then
+        error("Invalid number of inputs")
+    end
+    HzFlag= %F;
+    
+    select rhs
+     case 1 then
+         Fs=1; whole = "" ; n = 512 ; a = 1 ; 
+     case 2 then
+         Fs=1; whole = "" ; n = 512 ;  
+     case 3 then
+         Fs=1; whole = "" ;
+     case 4 then
+         Fs=1;
+    end
+   if (max(size(n)) > 1)
+    if (rhs > 4)
+      error("invalid inputs");
+    elseif (rhs > 3)
+      // grpdelay (B, A, F, Fs)
+      Fs     = whole;
+      HzFlag = %T;
+    else
+      // grpdelay (B, A, W)
+      Fs = 1;
+    end
+    w     = 2*%pi*n/Fs;
+    n  = max(size (w) ) * 2;
+    whole = "";
+  else
+    if (rhs < 5)
+      Fs = 1; // return w in radians per sample
+      if (rhs < 4)
+          whole = "" ;
+      elseif (~type(whole)==10)
+        Fs      = whole;
+        HzFlag  = %T;
+        whole   = "";
+      end
+      if (rhs < 3)
+        n = 512;
+      end
+      if (rhs < 2)
+        a = 1;
+      end
+    else
+      HzFlag = %T;
+    end
 
-//Parameters 
-//b:
-//a:
-//nfft:
-
-//Description
-//This is an Octave function.
-//This function computes the group delay of a filter.
-//[g, w] = grpdelay(b) returns the group delay g of the FIR filter with coefficients b. The response is evaluated at 512 angular frequencies between 0 and pi. w is a vector containing the 512 frequencies. 
-//[g, w] = grpdelay(b, a) returns the group delay of the rational IIR filter whose numerator has coefficients b and denominator coefficients a.
-//[g, w] = grpdelay(b, a, n) returns the group delay evaluated at n angular frequencies. 
-//[g, w] = grpdelay(b, a, n, ’whole’) evaluates the group delay at n frequencies between 0 and 2*pi.
-//[g, f] = grpdelay(b, a, n, Fs) evaluates the group delay at n frequencies between 0 and Fs/2.
-//[g, f] = grpdelay(b, a, n, ’whole’, Fs) evaluates the group delay at n frequencies between 0 and Fs.
-//[g, w] = grpdelay(b, a, w) evaluates the group delay at frequencies w (radians per sample).
-//[g, f] = grpdelay(b, a, f, Fs) evaluates the group delay at frequencies f (in Hz).
-//Examples
-//grpdelay(1,2,3)
-//ans =
-//   0.
-//   0.
-//   0.
-
-funcprot(0);
-rhs = argn(2)
-if(rhs<1 | rhs>5)
-error("Wrong number of input arguments.")
-end
-	select(rhs)
-	case 1 then
-	[gd,w] = callOctave("grpdelay",b)
-	case 2 then
-	[gd,w] = callOctave("grpdelay",b,a)
-	case 3 then
-	[gd,w] = callOctave("grpdelay",b,a,nfft)
-	case 4 then
-	[gd,w] = callOctave("grpdelay",b,a,nfft,whole)
-	case 5 then
-	[gd,w] = callOctave("grpdelay",b,a,nfft,whole,Fs)
-	end
+    if (isempty (n))
+      n = 512;
+    end
+    if ( type(whole) == 10 && strcmp (whole, "whole"))
+      n = 2*n;
+    end
+    
+    w = Fs*[0:n-1]/n;
+  end
+  if (~ HzFlag)
+    w = w * 2 * %pi;
+  end
+    a = a(:).';
+    b = b(:).';
+    oa = max(size(a)) -1;     // order of a(z)
+    if (oa < 0)             // a can be []
+      a  = 1;
+      oa = 0;
+    end
+    ob = max(size(b)) -1;     // order of b(z)
+    if (ob < 0)             // b can be [] as well
+      b  = 1;
+      ob = 0;
+    end
+    oc = oa + ob;           // order of c(z)
+    c   = conv (b, flipdim(conj (a),2));  // c(z) = b(z)*conj(a)(1/z)*z^(-oa)
+    cr  = c.*(0:oc);                    // cr(z) = derivative of c wrt 1/z
+    num = fft1 (cr, n);
+    den = fft1 (c, n);
+    minmag    = 10*%eps;
+    polebins  = find (abs (den) < minmag);
+  for b = polebins
+      warning ("signal:grpdelay-singularity", "grpdelay: setting group delay to 0 at singularity");
+      num(b) = 0;
+      den(b) = 1;
+      //// try to preserve angle:
+      //// db = den(b);
+      //// den(b) = minmag*abs(num(b))*exp(j*atan2(imag(db),real(db)));
+      //// warning(sprintf('grpdelay: den(b) changed from %f to %f',db,den(b)));
+    end
+    gd = real (num ./ den) - oa;
+    if ( type(whole) == 10 && strcmp (whole, "whole"))
+      ns = n/2; // Matlab convention ... should be n/2 + 1
+      gd = gd(1:ns);
+      w  = w(1:ns);
+    else
+      ns = n; // used in plot below
+    end
+    //// compatibility
+    gd = gd(:);
+    w  = w(:);
+    if (lhs == 1)
+      if (HzFlag)
+        funits = "Hz";
+      else
+        funits = "radian/sample";
+      end
+      disp(" Plotting ")
+      plot (w(1:ns), gd(1:ns));
+      xlabel (["Frequency (" funits ")"]);
+      ylabel ("Group delay (samples)");
+  end
 endfunction
diff --git a/macros/hilbert1.sci b/macros/hilbert1.sci
index fbcb136..476c00c 100644
--- a/macros/hilbert1.sci
+++ b/macros/hilbert1.sci
@@ -1,39 +1,106 @@
-function h= hilbert1(f, varargin)
-//Analytic extension of real valued signal.
-//Calling Sequence
-// h= hilbert1(f)
-// h= hilbert1(f,N)
-// h= hilbert1(f,N,dim)
-//Parameters 
-//f: real or complex valued scalar or vector
-//N: The result will have length N
-//dim : It analyses the input in this dimension
-//Description
-//h = hilbert1 (f) computes the extension of the real valued signal f to an analytic signal. If f is a matrix, the transformation is applied to each column. For N-D arrays, the transformation is applied to the first non-singleton dimension.
-//
-//real (h) contains the original signal f. imag (h) contains the Hilbert transform of f.
-//
-//hilbert1 (f, N) does the same using a length N Hilbert transform. The result will also have length N.
-//
-//hilbert1 (f, [], dim) or hilbert1 (f, N, dim) does the same along dimension dim.
-//Examples
-//## notice that the imaginary signal is phase-shifted 90 degrees
-// t=linspace(0,10,256);
-// z = hilbert1(sin(2*pi*0.5*t));
-// grid on; plot(t,real(z),';real;',t,imag(z),';imag;');
-
-funcprot(0);
-rhs= argn(2);
-if(rhs<1 | rhs>3)
-	error("Wrong number of Input Arguments")
-end
-
-select(rhs)
-	case 1 then
-		h= callOctave("hilbert", f);
-	case 2 then
-		h= callOctave("hilbert", f, varargin(1));
-	case 3 then
-		h= callOctave("hilbert", f, varargin(1), varargin(2));
-end
+/*Calling Sequence
+          h = hilbert1 (f)  
+          h = hilbert1 (f, N)
+          h = hilbert1 (f, N, dim)
+Description
+          Analytic extension of real valued signal.
+          h = hilbert (f) computes the extension of the real valued signal f to an analytic signal.
+          If f is a matrix, the transformation is applied to each column.
+          For N-D arrays, the transformation is applied to the first non-singleton dimension.
+          real (h) contains the original signal f. imag (h) contains the Hilbert transform of f.
+          hilbert (f, N) does the same using a length N Hilbert transform. The result will also have length N.
+          hilbert (f, [], dim) or hilbert (f, N, dim) does the same along dimension dim.
+Dependencies
+          fft1, ifft1, ipermute
+Example
+            //the magnitude of the hilbert transform eliminates the carrier
+            t=linspace(0,10,1024);
+            x=5*cos(0.2*t).*sin(100*t);
+            plot(t,x,t,abs(hilbert(x)));
+          */
+function f=hilbert1(f, N ,dim )
+  // ------ PRE: initialization and dimension shifting ---------
+  nargin = argn(2); 
+  if (nargin<1 || nargin>3)
+    error("Please enter valid number of inputs")
+  end
+  if ~isreal(f)
+    warning ('HILBERT: ignoring imaginary part of signal');
+    f = real (f);
+  end
+  D=ndims(f);
+  select nargin
+  case 1 then
+    N=[];
+    dim=[];
+  case 2 then
+    dim=[]
+  end
+  // Dummy assignment.
+  order=1;
+  if isempty(dim)
+    dim=1;
+    if sum(size(f)>1)==1
+      // We have a vector, find the dimension where it lives.
+      dim=find(size(f)>1);
+    end
+  else
+    if (length(dim)~=1 || ~or(type(dim)==[1 5 8]))
+      error('HILBERT: dim must be a scalar.');
+    end
+    if modulo(dim,1)~=0
+      error('HILBERT: dim must be an integer.');
+    end
+    if (dim<1) || (dim>D)
+      error('HILBERT: dim must be in the range from 1 to %d.',D);
+    end
+  end
+  if (length(N)>1 || ~or(type(N)==[1 5 8]))
+    error('N must be a scalar.');
+  elseif (~isempty(N) && modulo(N,1)~=0)
+    error('N must be an integer.');
+  end
+  if dim>1
+    order=[dim, 1:dim-1,dim+1:D];
+    // Put the desired dimension first.
+    f=permute(f,order);
+  end
+  Ls=size(f,1);
+  // If N is empty it is set to be the length of the transform.
+  if isempty(N)
+    N=Ls;
+  end
+  // moduloember the exact size for later and modify it for the new length
+  permutedsize=size(f);
+  permutedsize(1)=N;
+  // Reshape f to a matrix.
+  f=resize_matrix(f,size(f,1),length(f)/size(f,1));
+  W=size(f,2);
+  if ~isempty(N)
+    siz=size(f);
+    siz(1)=N;
+    f=resize_matrix(f,siz);
+  end
+  // ------- actual computation -----------------
+  if N>2
+    f=fft1(f);
+    if modulo(N,2)==0
+      f=[f(1,:);
+         2*f(2:N/2,:);
+         f(N/2+1,:);
+         zeros(N/2-1,W)];
+    else
+      f=[f(1,:);
+         2*f(2:(N+1)/2,:);
+         zeros((N-1)/2,W)];
+    end
+    f=ifft1(f);
+  end
+  // ------- POST: Restoration of dimensions ------------
+  // Restore the original, permuted shape.
+  f=matrix(f,permutedsize);
+  if dim>1
+    // Undo the permutation.
+    f=ipermute(f,order);
+  end
 endfunction
diff --git a/macros/idct1.sci b/macros/idct1.sci
index 5015187..3f007fa 100644
--- a/macros/idct1.sci
+++ b/macros/idct1.sci
@@ -1,27 +1,39 @@
+/*Description
+    This function computes the inverse discrete cosine transform of input X.
+    If N is given, then X is padded or trimmed to length N before computing the transform.
+    If X is a matrix, compute the transform along the columns of the the matrix.
+    The transform is faster if X is real-valued and even length.
+Calling Sequence
+    Y = idct1(X)
+    Y = idct1(X, N)
+Parameters
+    X: Matrix or integer
+    N: If N is given, then X is padded or trimmed to length N before computing the transform.
+Examples
+    idct1([1,3,6])
+    ans =
+        5.1481604  - 4.3216292    0.9055197
+*/
 function y = idct1(x,n)
-//Compute the inverse discrete cosine transform of input.
-//Calling Sequence
-//Y = idct1(X)
-//Y = idct1(X, N)
-//Parameters
-//X: Matrix or integer
-//N: If N is given, then X is padded or trimmed to length N before computing the transform.
-//Description
-// This function computes the inverse discrete cosine transform of input X.  If N is given, then X is padded or trimmed to length N before computing the transform.  If X is a matrix, compute the transform along the columns of the the matrix.  The transform is faster if X is real-valued and even length.
-//Examples
-//idct1([1,3,6])
-//ans = 
-//     5.1481604  - 4.3216292    0.9055197 
-funcprot(0);
-rhs=argn(2);
-if (rhs<1 | rhs>2) then
-    error("Wrong number of input arguments.");
-end
-select(rhs)
-case 1 then
-    y=callOctave("idct",x);
-case 2 then 
-    y=callOctave("idct",x,n);
-end
-
+    funcprot(0);
+    rhs=argn(2);
+    if (rhs<1 | rhs>2) then
+        error("Wrong number of input arguments.");
+    end
+    nsdim=1;
+    siz=size(x);
+    len=length(siz);
+    for i=1:len
+        if siz(i) ~= 1 then 
+            nsdim=i//calculating along non-singlton dimension
+            break;
+        end;
+    end;
+    select(rhs)
+    case 1 then
+        y=idct(x,nsdim);
+    case 2 then
+        siz(nsdim)=n;
+        y=idct(resize_matrix(x,siz),nsdim)
+    end;
 endfunction
diff --git a/macros/idct2.sci b/macros/idct2.sci
index c48980a..890927f 100644
--- a/macros/idct2.sci
+++ b/macros/idct2.sci
@@ -1,29 +1,31 @@
-function y = idct2(x,varargin)
-//This function computes the inverse 2-D discrete cosine transform of input matrix.
-//Calling Sequence
-//Y = idct2(X)
-//Y = idct2(X, M, N)
-//Y = idct2(X, [M, N])
-//Parameters
-//X: Matrix or integer
-//M, N: If specified Matrix X is padded with M rows and N columns.
-//Description
-// This function computes the inverse 2-D discrete cosine transform of matrix X. If M and N are specified, the input is either padded or truncated to have M rows and N columns.
-//Examples
-//idct2(3, 4, 6)
-//ans = 
-//     2.811261   0.612372  -0.525856   0.250601   0.612372  -0.086516
-funcprot(0);
-rhs=argn(2);
-if (rhs<1 | rhs>3) then
-    error("Wrong number of input arguments.");
-end
-select(rhs)
-case 1 then
-    y=callOctave("idct2",x)
-case 2 then 
-    y=callOctave("idct2",x,varargin(1))
-case 3 then
-    y=callOctave("idct2",x,varargin(1),varargin(2))
-end
+/*Description
+    This function computes the inverse 2-D discrete cosine transform of matrix X. If M and N are specified, the input is either padded or truncated to have M rows and N columns.
+Calling Sequence
+        Y = idct2(X)
+        Y = idct2(X, M, N)
+        Y = idct2(X, [M, N])
+Parameters
+        X: Matrix or integer
+        M, N: If specified Matrix X is padded with M rows and N columns.
+Examples
+     idct2(3, 4, 6)
+     ans =
+         2.811261   0.612372  -0.525856   0.250601   0.612372  -0.086516 */
+function y = idct2 (x, m, n)
+   funcprot(0);
+   rhs=argn(2);
+   select (rhs)
+   case 1 then
+       [m,n]=size(x);
+   case 2 then    
+        n=m(2);
+        m=m(1);
+   end
+   if m==1 then
+       y=idct1(x.',n).';
+   elseif n==1 then
+       y=idct1(x,m);
+   else
+       y=idct1(idct1(x,m).',n).';        
+   end
 endfunction
diff --git a/macros/idst1.sci b/macros/idst1.sci
index 98fc874..f7fdc6d 100644
--- a/macros/idst1.sci
+++ b/macros/idst1.sci
@@ -1,27 +1,26 @@
-function y = idst1(x,varargin)
-//This function computes the inverse type I discrete sine transform.
-//Calling Sequence
-//Y = idst(X)
-//Y = idst(X, N)
-//Parameters
-//X: Matrix or integer
-//N: If N is given, then X is padded or trimmed to length N before computing the transform.
-//Description
-//This function computes the inverse type I discrete sine transform of Y. If N is given, then Y is padded or trimmed to length N before computing the transform. If Y is a matrix, compute the transform along the columns of the the matrix.
-//Examples
-//idst([1,3,6])
-//ans = 
-//     3.97487  -2.50000   0.97487 
-funcprot(0);
-rhs=argn(2);
-if(rhs<1 | rhs>2) then
-    error("Wrong number of input arguments.");
-end
-select(rhs)
-case 1 then
-    y=callOctave("idst",x);
-case 2 then
-    y=callOctave("idst",x,varargin(1));
-end
-
+/*Description
+        This function computes the inverse type I discrete sine transform of X If N is given,
+        then X is padded or trimmed to length N before computing the transform.
+        If X is a matrix, compute the transform along the columns of the the matrix.
+Calling Sequence
+        Y = idst1(X)
+        Y = idst1(X, N)
+Parameters
+        X: Matrix or integer
+        N: If N is given, then X is padded or trimmed to length N before computing the transform.
+Examples
+        idst1([1,3,6])
+    ans = 
+         3.97487  -2.50000   0.97487 */
+function x = idst1 (y, n)
+  nargin=argn(2)
+  if (nargin < 1 || nargin > 2)
+    error("invalid input arguments")
+  end
+  if nargin == 1,
+    n = size(y,1);
+    if n==1, n = size(y,2); end
+  end
+  x = dst1(y, n) * 2/(n+1);
 endfunction
+
diff --git a/macros/ifft1.sci b/macros/ifft1.sci
index 3e94934..7c66b19 100644
--- a/macros/ifft1.sci
+++ b/macros/ifft1.sci
@@ -1,51 +1,63 @@
-function res =  ifft1 (x, n, dim)
-//Calculates the inverse discrete Fourier transform of a matrix using Fast Fourier Transform algorithm.
-//Calling Sequence
-//ifft (x, n, dim)
-//ifft (x, n)
-//ifft (x)
-//Parameters 
-//x: input matrix
-//n: Specifies the number of elements of x to be used
-//dim: Specifies the dimention of the matrix along which the inverse FFT is performed
-//Description
-//This is an Octave function.
-//Description
-//This is an Octave function.
-//The inverse FFT is calculated along the first non-singleton dimension of the array. Thus, inverse FFT is computed for each column of x.
-//
-//n is an integer specifying the number of elements of x to use. If n is larger than dimention along. which the inverse FFT is calculated, then x is resized and padded with zeros.
-//Similarly, if n is smaller, then x is truncated.
-//
-//dim is an integer specifying the dimension of the matrix along which the inverse FFT is performed.
-//Examples
-//x = [1 2 3; 4 5 6; 7 8 9]
-//n = 3
-//dim = 2
-//ifft1 (x, n, dim)
-//ans =
-//
-//   2.00000 + 0.00000i  -0.50000 - 0.28868i  -0.50000 + 0.28868i
-//   5.00000 + 0.00000i  -0.50000 - 0.28868i  -0.50000 + 0.28868i
-//   8.00000 + 0.00000i  -0.50000 - 0.28868i  -0.50000 + 0.28868i
-
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 1 | rhs > 3)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-	
-	case 1 then
-	res = callOctave("ifft", x)
-
-	case 2 then
-	res = callOctave("ifft", x, n)
-
-	case 3 then
-	res = callOctave("ifft", x, n, dim)
-
-	end
-endfunction
+/* Description
+    Calculates the inverse discrete Fourier transform of a matrix using Fast Fourier Transform algorithm.
+    The inverse FFT is calculated along the first non-singleton dimension of the array. Thus, inverse FFT is computed for each column of x.
+    n is an integer specifying the number of elements of x to use. If n is larger than dimention along. which the inverse FFT is calculated, then x is resized and padded with zeros.
+    Similarly, if n is smaller, then x is truncated.
+    dim is an integer specifying the dimension of the matrix along which the inverse FFT is performed.
+Calling Sequence
+        ifft1 (x)
+        ifft1 (x, n)
+        ifft1 (x, n, dim)
+Parameters
+        x: input matrix
+        n: Specifies the number of elements of x to be used
+        dim: Specifies the dimention of the matrix along which the inverse FFT is performed
+Examples
+        x = [1 2 3; 4 5 6; 7 8 9]
+        n = 3
+        dim = 2
+        ifft1 (x, n, dim)
+        ans =
+        2.00000 + 0.00000i  -0.50000 - 0.28868i  -0.50000 + 0.28868i
+        5.00000 + 0.00000i  -0.50000 - 0.28868i  -0.50000 + 0.28868i
+        8.00000 + 0.00000i  -0.50000 - 0.28868i  -0.50000 + 0.28868i */
+function res =  ifft1 (x, n, dim)
+    funcprot(0);
+    lhs = argn(1)
+    rhs = argn(2)
+    if (rhs < 1 | rhs > 3)
+        error("Wrong number of input arguments.")
+    end
+    dimension = size(x);
+    nsdim = 1;
+    for i = 1:length(dimension)
+        if dimension(i) ~= 1 then
+            nsdim = i;
+            break;
+        end
+    end
+    select(rhs)
+    case 1 then
+        res=fft(x,1,nsdim)
+    case 2 then
+        if isempty(n) then
+            res=fft(x,1,nsdim)
+        else
+            dimension(nsdim)=n;
+            res=fft(resize_matrix(x,dimension),1,nsdim)
+        end
+    case 3 then
+        if isempty(n) then
+            res=fft(x,1,dim)
+        else
+            if (length(dimension) <dim )then
+                error("ifft1: DIM must be a valid dimension along which to perform FFT")
+            end
+            dimension(dim)=n;
+           res=resize_matrix(x,dimension);
+            res=fft(res,1,dim);
+        end
+    end
+endfunction
+
+    
diff --git a/macros/ifft2.sci b/macros/ifft2.sci
index 8ecb1a5..762631d 100644
--- a/macros/ifft2.sci
+++ b/macros/ifft2.sci
@@ -1,45 +1,55 @@
-function res = ifft2 (A, m, n)
-//Calculates the inverse two-dimensional discrete Fourier transform of A using a Fast Fourier Transform algorithm.
-//Calling Sequence
-//ifft2 (A, m, n)
-//ifft2 (A)
-//Parameters 
-//A: input matrix
-//m: number of rows of A to be used
-//n: number of columns of A to be used
-//Description
-//This is an Octave function.
-//It performs inverse two-dimensional FFT on the matrix A. m and n may be used specify the number of rows and columns of A to use. If either of these is larger than the size of A, A is resized and padded with zeros.
-//If A is a multi-dimensional matrix, each two-dimensional sub-matrix of A is treated separately.
-//Examples
-//x = [1 2 3; 4 5 6; 7 8 9]
-//m = 4
-//n = 4
-//ifft2 (A, m, n)
-//ans =
-//
-//   2.81250 + 0.00000i  -0.37500 + 0.93750i   0.93750 + 0.00000i  -0.37500 - 0.93750i
-//  -1.12500 + 0.93750i  -0.31250 - 0.50000i  -0.37500 + 0.31250i   0.31250 + 0.25000i
-//   0.93750 + 0.00000i  -0.12500 + 0.31250i   0.31250 + 0.00000i  -0.12500 - 0.31250i
-//  -1.12500 - 0.93750i   0.31250 - 0.25000i  -0.37500 - 0.31250i  -0.31250 + 0.50000i
-
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 1 | rhs > 3)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-	
-	case 1 then
-	res = callOctave("ifft2", A)
-
-	case 2 then
-	error("Wrong number of input arguments.")
-
-	case 3 then
-	res = callOctave("ifft2", A, m, n)
-
-	end
-endfunction
+Author: Abinash Singh <abinashsinghlalotra@gmail.com>
+*/
+/*Description
+    Calculates the inverse two-dimensional discrete Fourier transform of A using a Fast Fourier Transform algorithm.
+    It performs inverse two-dimensional FFT on the matrix A. m and n may be used specify the number of rows and columns of A to use. If either of these is larger
+    than the size of A, A is resized and padded with zeros.
+    If A is a multi-dimensional matrix, each two-dimensional sub-matrix of A is treated separately.
+Calling Sequence
+     ifft2 (A)
+     ifft2 (A, m, n)
+Parameters
+     A: input matrix
+     m: number of rows of A to be used
+     n: number of columns of A to be used
+Examples
+     A= [1 2 3; 4 5 6; 7 8 9]
+     m = 4
+     n = 4
+     ifft2 (A, m, n) --functionCall
+     ans =
+        2.81250 + 0.00000i  -0.37500 + 0.93750i   0.93750 + 0.00000i  -0.37500 - 0.93750i
+        -1.12500 + 0.93750i  -0.31250 - 0.50000i  -0.37500 + 0.31250i   0.31250 + 0.25000i
+        0.93750 + 0.00000i  -0.12500 + 0.31250i   0.31250 + 0.00000i  -0.12500 - 0.31250i
+        -1.12500 - 0.93750i   0.31250 - 0.25000i  -0.37500 - 0.31250i  -0.31250 + 0.50000i
+*/
+function res = ifft2(A, m, n)
+    funcprot(0);
+    lhs = argn(1)
+    rhs = argn(2)
+    if (rhs < 1 | rhs > 3)
+        error("Wrong number of input arguments.")
+    end
+    size_x=size(A)
+    len=length(size_x)
+    //to figure out size of multidimensional array
+    if len>2 then
+        last_dim=size_x(len)
+    else
+        last_dim=1
+    end    
+    select(rhs)
+        case 1 then
+        res=[]
+        for i=1:last_dim //treating each submatrix seperately
+            res(:,:,i)=fft(A(:,:,i),1)
+        end
+    case 2 then
+        error("Wrong number of input arguments.")
+    case 3 then
+        res=[]
+        for i=1:last_dim 
+            res(:,:,i)=fft(resize_matrix(A(:,:,i),m,n),1)
+        end
+    end
+endfunction
diff --git a/macros/ifftn.sci b/macros/ifftn.sci
index 3d26c04..941dc67 100644
--- a/macros/ifftn.sci
+++ b/macros/ifftn.sci
@@ -1,27 +1,42 @@
-function y = ifftn(A, varargin)
-//Compute the inverse N-dimensional discrete Fourier transform of A using a Fast Fourier Transform (FFT) algorithm.
-//Calling Sequence
-//Y = ifftn(A)
-//Y = ifftn(A, size)
-//Parameters
-//A: Matrix 
-//Description
-//Compute the inverse N-dimensional discrete Fourier transform of A using a Fast Fourier Transform (FFT) algorithm. The optional vector argument SIZE may be used specify the dimensions of the array to be used.  If an element of SIZE is smaller than the corresponding dimension of A, then the dimension of A is truncated prior to performing the inverse FFT. Otherwise, if an element of SIZE is larger than the corresponding dimension then A is resized and padded with zeros.
-//Examples
-//ifftn([2,3,4])
-//ans = 
-//    3.  - 0.5 - 0.2886751i  - 0.5 + 0.2886751i 
-funcprot(0);
-funcprot(0);
-rhs = argn(2)
-if(rhs<1 | rhs>2)
-error("Wrong number of input arguments.");
-end
-
-	select(rhs)
-	case 1 then
-	y = callOctave("ifftn",A);
-	case 2 then
-	y = callOctave("ifftn",A, varargin(1));
-	end
-endfunction
+*/
+/*Description
+    Compute the inverse N-dimensional discrete Fourier transform of A using a Fast Fourier Transform (FFT) algorithm.
+    The optional vector argument SIZE may be used specify the dimensions of the matrix to be used.
+    If an element of SIZE is smaller than the corresponding dimension of A, then the dimension of A is truncated prior to performing the inverse FFT.
+    Otherwise, if an element of SIZE is larger than the corresponding dimension then A is resized and padded with zeros.
+Calling Sequence
+     Y = ifftn(A)
+     Y = ifftn(A, size)
+Parameters
+     A: Matrix
+     SIZE : (optional) dimension of matrix to be used
+Examples
+     ifftn([2,3,4])
+     ans =
+            3.  - 0.5 - 0.2886751i  - 0.5 + 0.2886751i */
+function y = ifftn(A, SIZE)
+    funcprot(0);
+    funcprot(0);
+    rhs = argn(2)
+    if(rhs<1 | rhs>2)
+        error("Wrong number of input arguments.");
+    end
+    select(rhs)
+    case 1 then
+        y=fft(A,1);
+    case 2 then
+        // Check if A needs resizing
+        if size(A) == SIZE then
+            // No resizing needed
+            break;
+        elseif length(size(A)) ~= length(SIZE) then
+            error("Output size must have at least Ndims");
+        else
+            // Resize A using the resize_matrix function
+            A = resize_matrix(A, SIZE);
+        end
+        y = fft(A,1);
+    end
+    y = clean( y ) ;
+endfunction
+
diff --git a/macros/ipermute.sci b/macros/ipermute.sci
new file mode 100644
index 0000000..ed78383
--- /dev/null
+++ b/macros/ipermute.sci
@@ -0,0 +1,25 @@
+/*
+Description
+    The inverse of the permute function.
+    The expression
+        ipermute (permute (A, perm), perm)
+    returns the original array A.
+Calling Sequence
+    ipermute (A, iperm)
+*/
+function B = ipermute(A, perm)
+    // ipermute : Inverse permute the dimensions of a matrix A.
+    // B = ipermute(A, perm) returns the array A with dimensions inverted
+    // according to the permutation vector `perm`.
+    // Validate the permutation vector
+    if max(size(perm)) ~= ndims(A) || or(gsort(perm, "g", "i") ~= 1:ndims(A))
+        error('Permutation vector must contain unique integers from 1 to ndims(A).');
+    end
+    // Compute the inverse permutation vector
+    invPerm = zeros(size(perm,1),size(perm , 2));
+    for i = 1:max(size(perm))
+        invPerm(perm(i)) = i;
+    end
+    // Use the permute function with the inverse permutation
+    B = permute(A, invPerm);
+endfunction
diff --git a/macros/mscohere.sci b/macros/mscohere.sci
index 357da58..b9a45f8 100644
--- a/macros/mscohere.sci
+++ b/macros/mscohere.sci
@@ -1,82 +1,42 @@
-function [PXX, FREQ] = mscohere (X, Y, WINDOW, OVERLAP, NFFT, FS, RANGE)
-//It estimate (mean square) coherence of signals x and y. 
-//Calling Sequence
-//[Pxx, freq] = mscohere (x, y)
-//[Pxx, freq] = mscohere (x, y, window)
-//[Pxx, freq] = mscohere (x, y, window, overlap)
-//[Pxx, freq] = mscohere (x, y, window, overlap, Nfft)
-//[Pxx, freq] = mscohere (x, y, window, overlap, Nfft, Fs)
-//[Pxx, freq] = mscohere (x, y, window, overlap, Nfft, Fs, range)
-//mscohere (...)
-//Description
-//This function estimate (mean square) coherence of signals x and y. 
-//Examples
-//[Pxx, freq] = mscohere(4,5)
-//ans =
-//PXX = 
-//    Nan
-//     1
-//FREQ =
-//     0
-//     0.5
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 2 | rhs > 7)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-	
-	case 2 then
-		if(lhs==0)
-		callOctave("mscohere",X,Y)
-		elseif(lhs==2)
-		[PXX, FREQ] = callOctave("mscohere",X,Y)
-		else
-		error("Wrong number of output arguments.")
-		end
-
-	case 3 then
-		if(lhs==0)
-		callOctave("mscohere",X,Y,WINDOW)
-		elseif(lhs==2)
-		[PXX, FREQ] = callOctave("mscohere",X,Y,WINDOW)
-		else
-		error("Wrong number of output arguments.")
-	       	end
-	case 4 then
-		if(lhs==0)
-		callOctave("mscohere",X,Y,WINDOW,OVERLAP)
-		elseif(lhs==2)
-		[PXX, FREQ] = callOctave("mscohere",X,Y,WINDOW,OVERLAP)
-		else
-		error("Wrong number of output arguments.")
-	       	end
-	case 5 then
-		if(lhs==0)
-		callOctave("mscohere",X,Y,WINDOW,OVERLAP,NFFT)
-		elseif(lhs==2)
-		[PXX, FREQ] = callOctave("mscohere",X,Y,WINDOW,OVERLAP,NFFT)
-		else
-		error("Wrong number of output arguments.")
-	       	end
-	case 6 then
-		if(lhs==0)
-		callOctave("mscohere",X,Y,WINDOW,OVERLAP,NFFT,FS)
-		elseif(lhs==2)
-		[PXX, FREQ] = callOctave("mscohere",X,Y,WINDOW,OVERLAP,NFFT,FS)
-		else
-		error("Wrong number of output arguments.")
-	       	end
-	case 7 then
-		if(lhs==0)
-		callOctave("mscohere",X,Y,WINDOW,OVERLAP,NFFT,FS,RANGE,)
-		elseif(lhs==2)
-		[PXX, FREQ] = callOctave("mscohere",X,Y,WINDOW,OVERLAP,NFFT,FS,RANGE)
-		else
-		error("Wrong number of output arguments.")
-	       	end
-	end
-endfunction
+/*Description:
+Estimate (mean square) coherence of signals x and y. Use the Welch (1967) periodogram/FFT method.
+Calling Sequence:
+        [Pxx, freq] = mscohere (x, y) 
+        […] = mscohere (x, y, window) 
+        […] = mscohere (x, y, window, overlap) 
+        […] = mscohere (x, y, window, overlap, Nfft) 
+        […] = mscohere (x, y, window, overlap, Nfft, Fs) 
+        […] = mscohere (x, y, window, overlap, Nfft, Fs, range) 
+        mscohere (…)
+See "help pwelch" for description of arguments, hints and references
+Dependencies : pwelch
+*/
+function varargout = mscohere(varargin)
+    // Check fixed argument
+    if (nargin < 2 || nargin > 7)
+      error("Invalid number of arguments");
+    end
+    nvarargin = length(varargin);
+    // remove any pwelch RESULT args and add 'cross'
+    for iarg=1:nvarargin
+      arg = varargin(iarg);
+      if ( ~isempty(arg) && ( type(arg) == 10 ) && ( ~strcmp(arg,'power') || ...
+             ~strcmp(arg,'cross') || ~strcmp(arg,'trans') || ...
+             ~strcmp(arg,'coher') || ~strcmp(arg,'ypower') ))
+        varargin(iarg) = [];
+      end
+    end
+    varargin(nvarargin+1) = 'coher';
+    disp(varargin)
+    if ( nargout==0 )
+      pwelch(varargin(:));
+    elseif ( nargout==1 )
+      Pxx = pwelch(varargin(:));
+      varargout(1) = Pxx;
+    elseif ( nargout>=2 )
+      [Pxx,f] = pwelch(varargin(:));
+      varargout(1) = Pxx;
+      varargout(2) = f;
+    end
+  endfunction
 
diff --git a/macros/ols.sci b/macros/ols.sci
new file mode 100644
index 0000000..f2b0b35
--- /dev/null
+++ b/macros/ols.sci
@@ -0,0 +1,70 @@
+/*
+Description:
+            Ordinary least squares estimation.
+            OLS applies to the multivariate model y = x*b + e with mean (e) = 0 and
+            cov (vec (e)) = kron (s, I). where y is a t by p matrix, x is a t by k matrix, b is a k by p matrix, and e is a t by p matrix.
+            Each row of y and x is an observation and each column a variable.
+            The return values beta, sigma, and r are defined as follows.
+Calling Sequence:
+            [beta, sigma, r] = ols (y, x)
+Arguments:
+            beta:
+                The OLS estimator for b. beta is calculated directly via inv (x'*x) * x' * y if the matrix x'*x is of full rank. Otherwise, beta = pinv (x) * y where pinv (x) denotes the pseudoinverse of x.
+            sigma
+                The OLS estimator for the matrix s,
+                sigma = (y-x*beta)'* (y-x*beta) / (t-rank(x))
+            r
+                The matrix of OLS residuals, r = y - x*beta.
+*/
+function [bt, sigma, r] = ols (y, x)
+    function [u , p] = formatted_chol(z)
+      //p flags whether the matrix A was positive definite and chol does not fail. A zero value of p indicates that matrix A is positive definite and R gives the factorization. Otherwise, p will have a positive value.
+      ierr  = execstr (" u = chol (z);",'errcatch' )
+      if ( ierr == 29 )
+        p = %T ;
+        warning("chol: Matrix is not positive definite.")
+        u = %nan
+      else
+        p = %F ;
+      end
+    endfunction
+    nargin = argn ( 2 )
+    if (nargin ~= 2)
+      error ("null");
+    end
+    if (~ (or (type(x) == [ 1 5 8] ) && or (type(y) == [1 5 8])))
+      error ("ols: X and Y must be numeric matrices or vectors");
+    end
+    if (ndims (x) ~= 2 || ndims (y) ~= 2)
+      error ("ols: X and Y must be 2-D matrices or vectors");
+    end
+    [nr, nc] = size (x);
+    [ry, cy] = size (y);
+    if (nr ~= ry)
+      error ("ols: number of rows of X and Y must be equal");
+    end
+    if (type(x) == 8)
+      x = double (x);
+    end
+    if ( type(y) == 8 )
+      y = double (y);
+    end
+    // Start of algorithm
+    z = x' * x;
+    [u, p] = formatted_chol (z);
+    if (p)
+      bt = pinv (x) * y;
+    else
+      bt = u \ (u' \ (x' * y));
+    end
+    if (nargout > 1)
+      r = y - x * bt;
+      // z is of full rank, avoid the SVD in rnk
+      if (p == 0)
+        rnk = size (z,2);
+      else
+        rnk = rank (z);
+      end
+      sigma = r' * r / (nr - rnk);
+    end
+endfunction
diff --git a/macros/pwelch.sci b/macros/pwelch.sci
index 3e4aa8d..e1c0e17 100644
--- a/macros/pwelch.sci
+++ b/macros/pwelch.sci
@@ -1,67 +1,819 @@
-function [spectra, Pxx_ci, freqn] = pwelch(x, varargin)
-//Estimate power spectral density of data "x" by the Welch (1967) periodogram/FFT method.
-//Calling Sequence
-//[spectra,freq] = pwelch(x, window, overlap, Nfft, Fs, range, plot_type, detrend, sloppy)
-//[spectra,freq] = pwelch(x, y, window, overlap, Nfft, Fs, range, plot_type, detrend, sloppy, results)
-//[spectra,Pxx_ci,freq] = pwelch(x, window, overlap, Nfft, Fs, conf, range, plot_type, detrend, sloppy)
-//[spectra,Pxx_ci,freq] = pwelch(x, y, window, overlap, Nfft, Fs, conf, range, plot_type, detrend, sloppy, results)
-//Parameters
-//x: [non-empty vector] system-input time-series data
-//
-//y: [non-empty vector] system-output time-series data
-//
-//window: [real vector] of window-function values between 0 and 1; the data segment has the same length as the window. Default window shape is Hamming. [integer scalar] length of each data segment.  The default value is window=sqrt(length(x)) rounded up to the nearest integer power of 2; see 'sloppy' argument.
-//
-//overlap: [real scalar] segment overlap expressed as a multiple of window or segment length. 0 <= overlap < 1, The default is overlap=0.5 .
-//
-//Nfft: [integer scalar] Length of FFT.  The default is the length of the "window" vector or has the same value as the scalar "window" argument.  If Nfft is larger than the segment length, "seg_len", the data segment is padded with "Nfft-seg_len" zeros.  The default is no padding. Nfft values smaller than the length of the data segment (or window) are ignored silently.
-//
-//Fs:[real scalar] sampling frequency (Hertz); default=1.0
-//
-//conf: [real scalar] confidence level between 0 and 1.  Confidence intervals of the spectral density are estimated from scatter in the periodograms and are returned as Pxx_ci. Pxx_ci(:,1) is the lower bound of the confidence interval and Pxx_ci(:,2) is the upper bound.  If there are three return values, or conf is an empty matrix, confidence intervals are calculated for conf=0.95 . If conf is zero or is not given, confidence intervals are not calculated. Confidence intervals can be obtained only for the power spectral density of x; nothing else.
-//
-//range: 
-//'half',  'onesided' : frequency range of the spectrum is zero up to but not including Fs/2.  Power from negative frequencies is added to the positive side of the spectrum, but not at zero or Nyquist (Fs/2) frequencies.  This keeps power equal in time and spectral domains.  See reference [2]. 
-//'whole', 'twosided' : frequency range of the spectrum is -Fs/2 to Fs/2, with negative frequencies stored in "wrap around" order after the positive frequencies; e.g. frequencies for a 10-point 'twosided' spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
-//'shift', 'centerdc' : same as 'whole' but with the first half of the spectrum swapped with second half to put the zero-frequency value in the middle. (See "help fftshift". If data (x and y) are real, the default range is 'half', otherwise default range is 'whole'.
-//
-//plot_type:  
-//'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db': specifies the type of plot.  The default is 'plot', which means linear-linear axes. 'squared' is the same as 'plot'. 'dB' plots "10*log10(psd)".  This argument is ignored and a spectrum is not plotted if the caller requires a returned value.
-//
-//detrend: 'no-strip', 'none' -- do NOT remove mean value from the data 'short', 'mean' -- remove the mean value of each segment from each segment of the data.
-//'linear', -- remove linear trend from each segment of the data.
-//'long-mean' -- remove the mean value from the data before splitting it into segments.  This is the default.
-//
-//sloppy: FFT length is rounded up to the nearest integer power of 2 by zero padding.  FFT length is adjusted after addition of padding by explicit Nfft argument. The default is to use exactly the FFT and window.
-//Description
-//Estimate power spectral density of data "x" by the Welch (1967) periodogram/FFT method. The data is divided into segments.  If "window" is a vector, each segment has the same length as "window" and is multiplied by "window" before (optional) zero-padding and calculation of its periodogram. If "window" is a scalar, each segment has a length of "window" and a Hamming window is used. The spectral density is the mean of the periodograms, scaled so that area under the spectrum is the same as the mean square of the data.  This equivalence is supposed to be exact, but in practice there is a mismatch of up to 0.5% when comparing area under a periodogram with the mean square of the data.
-funcprot(0);
-rhs=argn(2);
-lhs=argn(1)
-if(rhs<9 | rhs>9) then
-    error("Wrong number of input arguments.");
-end
-if(rhs<11 | rhs>11) then
-    error("Wrong number of input arguments.");
-end
-if(lhs<2 | lhs>3) then
-    error("Wrong number of output arguments.");
-end
-select(rhs)
-case 9 then
-    select(lhs)
-        case 2 then
-            [spectra, freqn] = callOctave("pwelch", x, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6), varargin(7), varargin(8));
-        case 3 then
-            [spectra, Pxx_ci, freqn] = callOctave("pwelch", x, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6), varargin(7), varargin(8));
+/*
+Dependencies : fft1
+Description:
+            Estimate power spectral density of data "x" by the Welch (1967) periodogram/FFT method.
+            All arguments except "x" are optional.
+Calling Sequence:
+            The data is divided into segments. If "window" is a vector, each segment has the same length as "window" and is multiplied by "window" before (optional) zero-padding and calculation of its periodogram. If "window" is a scalar, each segment has a length of "window" and a Hamming window is used.
+            The spectral density is the mean of the periodograms, scaled so that area under the spectrum is the same as the mean square of the data. This equivalence is supposed to be exact, but in practice there is a mismatch of up to 0.5% when comparing area under a periodogram with the mean square of the data.
+            [spectra,freq] = pwelch(x,y,window,overlap,Nfft,Fs, range,plot_type,detrend,sloppy,results)
+                    Two-channel spectrum analyser. Estimate power spectral density, cross- spectral density, transfer function and/or coherence functions of time- series input data "x" and output data "y" by the Welch (1967) periodogram/FFT method.
+                    pwelch treats the second argument as "y" if there is a control-string argument "cross", "trans", "coher" or "ypower"; "power" does not force the 2nd argument to be treated as "y". All other arguments are optional. All spectra are returned in matrix "spectra".
+            [spectra,Pxx_ci,freq] = pwelch(x,window,overlap,Nfft,Fs,conf, range,plot_type,detrend,sloppy)
+            [spectra,Pxx_ci,freq] = pwelch(x,y,window,overlap,Nfft,Fs,conf, range,plot_type,detrend,sloppy,results)
+                    Estimates confidence intervals for the spectral density.
+                    See Hint (7) below for compatibility options.
+                    Confidence level "conf" is the 6th or 7th numeric argument. If "results" control-string arguments are used, one of them must be "power" when the "conf" argument is present; pwelch can estimate confidence intervals only for the power spectrum of the "x" data. It does not know how to estimate confidence intervals of the cross-power spectrum, transfer function or coherence; if you can suggest a good method, please send a bug report.
+            ARGUMENTS
+                    All but the first argument are optional and may be empty, except that the "results" argument may require the second argument to be "y".
+                    x : [non-empty vector] system-input time-series data
+                    y : [non-empty vector] system-output time-series data
+                    window : [real vector] of window-function values between 0 and 1; the data segment has the same length as the window. Default window shape is Hamming.
+                            [integer scalar] length of each data segment. The default value is window=sqrt(length(x)) rounded up to the nearest integer power of 2; see ’sloppy’ argument.
+                    overlap: [real scalar] segment overlap expressed as a multiple of window or segment length. 0 <= overlap < 1, The default is overlap=0.5 .
+                    Nfft : [integer scalar] Length of FFT. The default is the length of the "window" vector or has the same value as the scalar "window" argument. If Nfft is larger than the segment length, "seg_len", the data segment is padded with "Nfft-seg_len" zeros. The default is no padding. Nfft values smaller than the length of the data segment (or window) are ignored silently.
+                    Fs : [real scalar] sampling frequency (Hertz); default=1.0
+                    conf : [real scalar] confidence level between 0 and 1. Confidence intervals of the spectral density are estimated from scatter in the periodograms and are returned as Pxx_ci. Pxx_ci(:,1) is the lower bound of the confidence interval and Pxx_ci(:,2) is the upper bound. If there are three return values, or conf is an empty matrix, confidence intervals are calculated for conf=0.95 . If conf is zero or is not given, confidence intervals are not calculated. Confidence intervals can be obtained only for the power spectral density of x; nothing else.
+                    CONTROL-STRING ARGUMENTS – each of these arguments is a character string. Control-string arguments must be after the other arguments but can be in any order.
+                    range :
+                          ’half’, ’onesided’ : frequency range of the spectrum is zero up to but not including Fs/2. Power from negative frequencies is added to the positive side of the spectrum, but not at zero or Nyquist (Fs/2) frequencies. This keeps power equal in time and spectral domains. See reference [2].
+                          ’whole’, ’twosided’ : frequency range of the spectrum is -Fs/2 to Fs/2, with negative frequencies stored in "wrap around" order after the positive frequencies; e.g. frequencies for a 10-point ’twosided’ spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
+                          ’shift’, ’centerdc’ : same as ’whole’ but with the first half of the spectrum swapped with second half to put the zero-frequency value in the middle. (See "help fftshift".
+                          If data (x and y) are real, the default range is ’half’, otherwise default range is ’whole’.
+                    plot_type
+                            ’plot’, ’semilogx’, ’semilogy’, ’loglog’, ’squared’ or ’db’: specifies the type of plot. The default is ’plot’, which means linear-linear axes. ’squared’ is the same as ’plot’. ’dB’ plots "10*log10(psd)". This argument is ignored and a spectrum is not plotted if the caller requires a returned value.
+                            detrend
+                            ’no-strip’, ’none’ – do NOT remove mean value from the data
+                            ’short’, ’mean’ – remove the mean value of each segment from each segment of the data.
+                            ’linear’, – remove linear trend from each segment of the data.
+                            ’long-mean’ – remove the mean value from the data before splitting it into segments. This is the default.
+                            sloppy
+                            ’sloppy’: FFT length is rounded up to the nearest integer power of 2 by zero padding. FFT length is adjusted after addition of padding by explicit Nfft argument. The default is to use exactly the FFT and window/ segment lengths specified in argument list.
+                    results : specifies what results to return (in the order specified and as many as desired).
+                            ’power’ calculate power spectral density of "x"
+                            ’cross’ calculate cross spectral density of "x" and "y"
+                            ’trans’ calculate transfer function of a system with input "x" and output "y"
+                            ’coher’ calculate coherence function of "x" and "y"
+                            ’ypower’ calculate power spectral density of "y"
+                            The default is ’power’, with argument "y" omitted.
+            RETURNED VALUES:
+                    If return values are not required by the caller, the results are plotted and nothing is returned.
+                    spectra : [real-or-complex matrix] columns of the matrix contain results in the same order as specified by "results" arguments. 
+                              Each column contains one of the result vectors.
+                    Pxx_ci : [real matrix] estimate of confidence interval for power spectral density of x. First column is the lower bound. 
+                              Second column is the upper bound.
+                    freq   : [real column vector] frequency values
+            HINTS
+                    EMPTY ARGS: if you don’t want to use an optional argument you can leave it empty by writing its value as [].
+                    FOR BEGINNERS: 
+                            The profusion of arguments may make pwelch difficult to use, and an unskilled user can easily produce a meaningless result or can easily mis-interpret the result.
+                            With real data "x" and sampling frequency "Fs", the easiest and best way for a beginner to use pwelch is probably "pwelch(x,[],[],[],Fs)". 
+                            Use the "window" argument to control the length of the spectrum vector. For real data and integer scalar M, "pwelch(x,2*M,[],[],Fs)" gives an M+1 point spectrum. 
+                            Run "help pwelch".
+                    WINDOWING FUNCTIONS: 
+                            Without a window function, sharp spectral peaks can have strong sidelobes because the FFT of a data in a segment is in effect convolved with a rectangular window. 
+                            A window function which tapers off (gradually) at the ends produces much weaker sidelobes in the FFT. Hann (hanning), hamming, bartlett, blackman, flattopwin etc are available as separate Matlab/sigproc or Octave functions. 
+                            The sidelobes of the Hann window have a roll-off rate of 60dB/decade of frequency. The first sidelobe of the Hamming window is suppressed and is about 12dB lower than the first Hann sidelobe, but the roll-off rate is only 20dB/decade. You can inspect the FFT of a Hann window by plotting "abs(fft(postpad(hanning(256),4096,0)))". The default window is Hamming.
+                    ZERO PADDING: 
+                            Zero-padding reduces the frequency step in the spectrum, and produces an artificially smoothed spectrum. 
+                            For example, "Nfft=2*length(window)" gives twice as many frequency values, but adjacent PSD (power spectral density) values are not independent; adjacent PSD values are independent if "Nfft=length(window)", which is the default value of Nfft.
+                    REMOVING MEAN FROM SIGNAL: 
+                            If the mean is not removed from the signal there is a large spectral peak at zero frequency and the sidelobes of this peak are likely to swamp the rest of the spectrum. 
+                            For this reason, the default behavior is to remove the mean. However, the matlab pwelch does not do this.
+                    WARNING ON CONFIDENCE INTERVALS :
+                            Confidence intervals are obtained by measuring the sample variance of the periodograms and assuming that the periodograms have a Gaussian probability distribution. 
+                            This assumption is not accurate. If, for example, the data (x) is Gaussian, the periodogram has a Rayleigh distribution. 
+                            However, the confidence intervals may still be useful.
+                    COMPATIBILITY WITH Matlab R11, R12, etc
+                            When used without the second data (y) argument, 
+                            arguments are compatible with the pwelch of Matlab R12, R13, R14, 2006a and 2006b except that
+                                    1) overlap is expressed as a multiple of window length — effect of overlap scales with window length
+                                    2) default values of length(window), Nfft and Fs are more sensible, and
+                                    3) Goertzel algorithm is not available so Nfft cannot be an array of frequencies as in Matlab 2006b.
+                            Pwelch has four persistent Matlab-compatibility levels. Calling pwelch with an empty first argument sets the order of arguments and defaults specified above in the USAGE and ARGUMENTS section of this documentation.
+                                    prev_compat=pwelch([]);
+                                    [Pxx,f]=pwelch(x,window,overlap,Nfft,Fs,conf,...);
+                            Calling pwelch with a single string argument (as described below) gives compatibility with Matlab R11 or R12, or the R14 spectrum.welch defaults. The returned value is the PREVIOUS compatibility string.
+                            Matlab R11: For compatibility with the Matlab R11 pwelch:
+                                    prev_compat=pwelch('R11-');
+                                    [Pxx,f]=pwelch(x,Nfft,Fs,window,overlap,conf,range,units);
+                                    // units of overlap are "number of samples"
+                                    // defaults: Nfft=min(length(x),256), Fs=2*pi, length(window)=Nfft,
+                                    //           window=Hanning, do not detrend,
+                                    // N.B.  "Sloppy" is not available.
+                            Matlab R12: For compatibility with Matlab R12 to 2006a pwelch:
+                                    prev_compat=pwelch('R12+');
+                                    [Pxx,f]=pwelch(x,window,overlap,nfft,Fs,...);
+                                    // units of overlap are "number of samples"
+                                    // defaults: length(window)==length(x)/8, window=Hamming,
+                                    //           Nfft=max(256,NextPow2), Fs=2*pi, do not detrend
+                                    // NextPow2 is the next power of 2 greater than or equal to the
+                                    // window length. "Sloppy", "conf" are not available.  Default
+                                    // window length gives very poor amplitude resolution.
+                            To adopt defaults of the Matlab R14 "spectrum.welch" spectrum object associated "psd" method.
+                                    prev_compat=pwelch('psd');
+                                    [Pxx,f] = pwelch(x,window,overlap,Nfft,Fs,conf,...);
+                                    // overlap is expressed as a percentage of window length,
+                                    // defaults: length(window)==64, Nfft=max(256,NextPow2), Fs=2*pi
+                                    //           do not detrend
+                                    // NextPow2 is the next power of 2 greater than or equal to the
+                                    // window length. "Sloppy" is not available.
+                                    // Default window length gives coarse frequency resolution.
+        */
+function varargout = pwelch(x,varargin)
+    //
+    // COMPATIBILITY LEVEL
+    // Argument positions and defaults depend on compatibility level selected
+    // by calling pwelch without arguments or with a single string argument.
+    //   native:      compatib=1; prev_compat=pwelch(); prev_compat=pwelch([]);
+    //   matlab R11:  compatib=2; prev_compat=pwelch('R11-');
+    //   matlab R12:  compatib=3; prev_compat=pwelch('R12+');
+    //   spectrum.welch defaults:  compatib=4; prev_compat=pwelch('psd');
+    // In each case, the returned value is the PREVIOUS compatibility string.
+    //
+    nargin = argn(2)
+    compat_str = {[]; 'R11-'; 'R12+'; 'psd'};
+    global compatib
+    if ( isempty(compatib) || compatib<=0 || compatib>4 )
+      // legal values are 1, 2, 3, 4
+      compatib = 1;
     end
-case 11 then
-    select(lhs)
-        case 2 then
-            [spectra, freqn] = callOctave("pwelch", x, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6), varargin(7), varargin(8), varargin(9), varargin(10));
-        case 3 then
-            [spectra, Pxx_ci, freqn] = callOctave("pwelch", x, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6), varargin(7), varargin(8), varargin(9), varargin(10));
+    if ( nargin <= 0 )
+      error( 'pwelch: Need at least 1 arg. Use help pwelch.' );
+    elseif ( nargin==1 && (type(x) == 10 ) || isempty(x)) 
+      varargout(1) = compat_str{compatib};
+      if ( isempty(x) ) // native
+        compatib = 1;
+      elseif ( ~strcmp(x,'R11-') )
+        compatib = 2;
+      elseif ( ~strcmp(x,'R12+') )
+        compatib = 3;
+      elseif ( ~strcmp(x,'psd') )
+        compatib = 4;
+      else
+        error( 'pwelch: compatibility arg must be empty, R11-, R12+ or psd' );
+      end
+      // return
+    //
+    // Check fixed argument
+    elseif ( isempty(x) || ~isvector(x) )
+      error( 'pwelch: arg 1 (x) must be vector.' );
+    else
+      //  force x to be COLUMN vector
+      if ( size(x,1)==1 )
+        x=x(:);
+      end
+      //
+      // Look through all args to check if  cross PSD, transfer function or
+      // coherence is required.  If yes, the second arg is data vector "y".
+      arg2_is_y = 0;
+      x_len = max(size(x));
+      nvarargin = max(size(varargin));
+      for iarg=1:nvarargin
+        arg = varargin(iarg);
+        if ( ~isempty(arg) && (type(arg) == 10 ) && ...
+             ( ~strcmp(arg,'cross') || ~strcmp(arg,'trans') || ...
+               ~strcmp(arg,'coher') || ~strcmp(arg,'ypower') ))
+          // OK. Need "y". Grab it from 2nd arg.
+          arg = varargin(1);
+          if ( nargin<2 || isempty(arg) || ~isvector(arg) || max(size(arg))~=x_len )
+            error( 'pwelch: arg 2 (y) must be vector, same length as x.' );
+          end
+          // force  COLUMN vector
+          y = varargin(1)(:);
+          arg2_is_y = 1;
+          break;
+        end
+      end
+      //
+      // COMPATIBILITY
+      // To select default argument values, "compatib" is used as an array index.
+      // Index values are   1=native,  2=R11,  3=R12,  4=spectrum.welch
+      //
+      //  argument positions:
+      //  arg_posn = varargin index of window, overlap, Nfft, Fs and conf
+      //             args respectively, a value of zero ==>> arg does not exist
+      arg_posn = [1 2 3 4 5;  // native
+                  3 4 1 2 5;  // Matlab R11- pwelch
+                  1 2 3 4 0;  // Matlab R12+ pwelch
+                  1 2 3 4 5]; // spectrum.welch defaults
+      arg_posn  = arg_posn(compatib,:) + arg2_is_y;
+      //
+      //  SPECIFY SOME DEFAULT VALUES for (not all) optional arguments
+      //  Use compatib as array index.
+      //  Fs = sampling frequency
+      Fs        = [ 1.0 2*%pi 2*%pi 2*%pi ];
+      Fs        = Fs(compatib);
+      //  plot_type: 1='plot'|'squared'; 5='db'|'dB'
+      plot_type = [ 1 5 5 5 ];
+      plot_type = plot_type(compatib);
+      //  rm_mean: 3='long-mean'; 0='no-strip'|'none'
+      rm_mean   = [ 3 0 0 0 ];
+      rm_mean   = rm_mean(compatib);
+      // use max_overlap=x_len-1 because seg_len is not available yet
+      // units of overlap are different for each version:
+      //    fraction, samples, or percent
+      max_overlap = [ 0.95 x_len-1 x_len-1 95];
+      max_overlap = max_overlap(compatib);
+      // default confidence interval
+      //  if there are more than 2 return values and if there is a "conf" arg
+      conf      = 0.95 * (nargout>2) * (arg_posn(5)>0);
+      //
+      is_win    = 0;    // =0 means valid window arg is not provided yet
+      Nfft      = [];   // default depends on segment length
+      overlap   = [];   // WARNING: units can be //samples, fraction or percentage
+      range     = ~isreal(x) || ( arg2_is_y && ~isreal(y) );
+      is_sloppy = 0;
+      n_results = 0;
+      do_power  = 0;
+      do_cross  = 0;
+      do_trans  = 0;
+      do_coher  = 0;
+      do_ypower = 0;
+    //
+    //  DECODE AND CHECK OPTIONAL ARGUMENTS
+      end_numeric_args = 0;
+      for iarg = 1+arg2_is_y:nvarargin
+        arg = varargin(iarg);
+        if ( ( type (arg) == 10 ) )
+          // first string arg ==> no more numeric args
+          // non-string args cannot follow a string arg
+          end_numeric_args = 1;
+          //
+          // decode control-string arguments
+          if ( ~strcmp(arg,'sloppy') )
+            is_sloppy = ~is_win || is_win==1;
+          elseif ( ~strcmp(arg,'plot') || ~strcmp(arg,'squared') )
+            plot_type = 1;
+          elseif ( ~strcmp(arg,'semilogx') )
+            plot_type = 2;
+          elseif ( ~strcmp(arg,'semilogy') )
+            plot_type = 3;
+          elseif ( ~strcmp(arg,'loglog') )
+            plot_type = 4;
+          elseif ( ~strcmp(arg,'db') || ~strcmp(arg,'dB') )
+            plot_type = 5;
+          elseif ( ~strcmp(arg,'half') || ~strcmp(arg,'onesided') )
+            range = 0;
+          elseif ( ~strcmp(arg,'whole') || ~strcmp(arg,'twosided') )
+            range = 1;
+          elseif ( ~strcmp(arg,'shift') || ~strcmp(arg,'centerdc') )
+            range = 2;
+          elseif ( ~strcmp(arg,'long-mean') )
+            rm_mean = 3;
+          elseif ( ~strcmp(arg,'linear') )
+            rm_mean = 2;
+          elseif ( ~strcmp(arg,'short') || ~strcmp(arg,'mean') )
+            rm_mean = 1;
+          elseif ( ~strcmp(arg,'no-strip') || ~strcmp(arg,'none') )
+            rm_mean = 0;
+          elseif ( ~strcmp(arg, 'power' ) )
+            if ( ~do_power )
+              n_results = n_results+1;
+              do_power = n_results;
+            end
+          elseif ( ~strcmp(arg, 'cross' ) )
+            if ( ~do_cross )
+              n_results = n_results+1;
+              do_cross = n_results;
+            end
+          elseif ( ~strcmp(arg, 'trans' ) )
+            if ( ~do_trans )
+              n_results = n_results+1;
+              do_trans = n_results;
+            end
+          elseif ( ~strcmp(arg, 'coher' ) )
+            if ( ~do_coher )
+              n_results = n_results+1;
+              do_coher = n_results;
+            end
+          elseif ( ~strcmp(arg, 'ypower' ) )
+            if ( ~do_ypower )
+              n_results = n_results+1;
+              do_ypower = n_results;
+            end
+          else
+            error( 'pwelch: string arg %d illegal value: %s', iarg+1, arg );
+          end
+          // end of processing string args
+          //
+        elseif ( end_numeric_args )
+          if ( ~isempty(arg) )
+            // found non-string arg after a string arg ... oops
+            error( 'pwelch: control arg must be string' );
+          end
+        //
+        // first 4 optional arguments are numeric -- in fixed order
+        //
+        // deal with "Fs" and "conf" first because empty arg is a special default
+        // -- "Fs" arg -- sampling frequency
+        elseif ( iarg == arg_posn(4) )
+          if ( isempty(arg) )
+            Fs = 1;
+          elseif ( ~isscalar(arg) || ~isreal(arg) || arg<0 )
+            error( 'pwelch: arg %d (Fs) must be real scalar >0', iarg+1 );
+          else
+            Fs = arg;
+          end
+        //
+        //  -- "conf" arg -- confidence level
+        //    guard against the "it cannot happen" iarg==0
+        elseif ( arg_posn(5) && iarg == arg_posn(5) )
+          if ( isempty(arg) )
+            conf = 0.95;
+          elseif ( ~isscalar(arg) || ~isreal(arg) || arg < 0.0 || arg >= 1.0 )
+            error( 'pwelch: arg %d (conf) must be real scalar, >=0, <1',iarg+1 );
+          else
+            conf = arg;
+          end
+        //
+        // skip all empty args from this point onward
+        elseif ( isempty(arg) )
+          1;
+        //
+        //  -- "window" arg -- window function
+        elseif ( iarg == arg_posn(1) )
+          if ( isscalar(arg) )
+            is_win = 1;
+          elseif ( isvector(arg) )
+            is_win = max(size(arg));
+            if ( size(arg,2)>1 )  // vector must be COLUMN vector
+              arg = arg(:);
+            end
+          else
+            is_win = 0;
+          end
+          if ( ~is_win )
+            error( 'pwelch: arg %d (window) must be scalar or vector', iarg+1 );
+          elseif ( is_win==1 && ( ~isreal(arg) || fix(arg)~=arg || arg<=3 ) )
+            error( 'pwelch: arg %d (window) must be integer >3', iarg+1 );
+          elseif ( is_win>1 && ( ~isreal(arg) ) )
+            error( 'pwelch: arg %d (window) vector must be real and >=0',iarg+1);
+          end
+          window = arg;
+          is_sloppy = 0;
+        //
+        // -- "overlap" arg -- segment overlap
+        elseif ( iarg == arg_posn(2) )
+          if (~isscalar(arg) || ~isreal(arg) || arg<0 || arg>max_overlap )
+            error( 'pwelch: arg %d (overlap) must be real from 0 to %f', ...
+                   iarg+1, max_overlap );
+          end
+          overlap = arg;
+        //
+        // -- "Nfft" arg -- FFT length
+        elseif ( iarg == arg_posn(3) )
+          if ( ~isscalar(arg) || ~isreal(arg) || fix(arg)~=arg || arg<0 )
+            error( 'pwelch: arg %d (Nfft) must be integer >=0', iarg+1 );
+          end
+          Nfft = arg;
+        //
+        else
+          error( 'pwelch: arg %d  must be string', iarg+1 );
+        end
+      end
+      if ( conf>0 && (n_results && ~do_power ) )
+        error('pwelch: can give confidence interval for x power spectrum only' );
+      end
+      //
+      // end DECODE AND CHECK OPTIONAL ARGUMENTS.
+      //
+      // SETUP REMAINING PARAMETERS
+      // default action is to calculate power spectrum only
+      if ( ~n_results )
+        n_results = 1;
+        do_power = 1;
+      end
+      need_Pxx = do_power || do_trans || do_coher;
+      need_Pxy = do_cross || do_trans || do_coher;
+      need_Pyy = do_coher || do_ypower;
+      log_two = log(2);
+      nearly_one = 0.99999999999;
+      //
+      // compatibility-options
+      // provides exact compatibility with Matlab R11 or R12
+      //
+      // Matlab R11 compatibility
+      if ( compatib==2 )
+        if ( isempty(Nfft) )
+          Nfft = min( 256, x_len );
+        end
+        if ( is_win > 1 )
+          seg_len = min( max(size(window)), Nfft );
+          window = window(1:seg_len);
+        else
+          if ( is_win )
+            // window arg is scalar
+            seg_len = window;
+          else
+            seg_len = Nfft;
+          end
+          // make Hann window (don't depend on sigproc)
+          xx = seg_len - 1;
+          window = 0.5 - 0.5 * cos( (2*%pi/xx)*[0:xx].' );
+        end
+      //
+      // Matlab R12 compatibility
+      elseif ( compatib==3 )
+        if ( is_win > 1 )
+          // window arg provides window function
+          seg_len = max(size(window));
+        else
+          // window arg does not provide window function; use Hamming
+          if ( is_win )
+            // window arg is scalar
+            seg_len = window;
+          else
+            // window arg not available; use R12 default, 8 windows
+            // ignore overlap arg; use overlap=50% -- only choice that makes sense
+            // this is the magic formula for 8 segments with 50% overlap
+            seg_len = fix( (x_len-3)*2/9 );
+          end
+          // make Hamming window (don't depend on sigproc)
+          xx = seg_len - 1;
+          window = 0.54 - 0.46 * cos( (2*%pi/xx)*[0:xx].' );
+        end
+        if ( isempty(Nfft) )
+          Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) );
+        end
+      // Matlab R14 psd(spectrum.welch) defaults
+      elseif ( compatib==4 )
+        if ( is_win > 1 )
+          // window arg provides window function
+          seg_len = max(size(window));
+        else
+          // window arg does not provide window function; use Hamming
+          if ( is_win )
+            // window arg is scalar
+            seg_len = window;
+          else
+            // window arg not available; use default seg_len = 64
+            seg_len = 64;
+          end
+          // make Hamming window (don't depend on sigproc)
+          xx = seg_len - 1;
+          window = 0.54 - 0.46 * cos( (2*%pi/xx)*[0:xx].' );
+        end
+        // Now we know segment length,
+        // so we can set default overlap as number of samples
+        if ( ~isempty(overlap) )
+          overlap = fix(seg_len * overlap / 100 );
+        end
+        if ( isempty(Nfft) )
+          Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) );
+        end
+      //
+      // default compatibility level
+      else // if ( compatib==1 )
+        // calculate/adjust segment lenght, window function
+        if ( is_win > 1 )
+          // window arg provides window function
+          seg_len = max(size(window));
+        else
+          // window arg does not provide window function; use Hamming
+          if ( is_win )       // window arg is scalar
+            seg_len = window;
+          else
+            // window arg not available; use default length:
+            // = sqrt(max(size(x))) rounded up to nearest integer power of 2
+            if ( isempty(overlap) )
+              overlap=0.5;
+            end
+            seg_len = 2 ^ ceil( log(sqrt(x_len/(1-overlap)))*nearly_one/log_two );
+          end
+          // make Hamming window (don't depend on sigproc)
+          xx = seg_len - 1;
+          window = 0.54 - 0.46 * cos( (2*%pi/xx)*[0:xx].' );
+        end
+        // Now we know segment length,
+        // so we can set default overlap as number of samples
+        if ( ~isempty(overlap) )
+          overlap = fix(seg_len * overlap);
+        end
+        //
+        // calculate FFT length
+        if ( isempty(Nfft) )
+          Nfft = seg_len;
+        end
+        if ( is_sloppy )
+          Nfft = 2 ^ ceil( log(Nfft) * nearly_one / log_two );
+        end
+      end
+      // end of compatibility options
+      //
+      // minimum FFT length is seg_len
+      Nfft = max( Nfft, seg_len );
+      // Mean square of window is required for normalizing PSD amplitude.
+      win_meansq = (window.' * window) / seg_len;
+      //
+      // Set default or check overlap.
+      if ( isempty(overlap) )
+        overlap = fix(seg_len /2);
+      elseif ( overlap >= seg_len )
+        error( 'pwelch: arg (overlap=%d) too big. Must be <max(size(window)=%d',...
+               overlap, seg_len );
+      end
+      //
+      // Pad data with zeros if shorter than segment. This should not happen.
+      if ( x_len < seg_len )
+        x = [x; zeros(seg_len-x_len,1)];
+        if ( arg2_is_y )
+          y = [y; zeros(seg_len-x_len,1)];
+        end
+        x_len = seg_len;
+      end
+      // end SETUP REMAINING PARAMETERS
+      //
+      //
+      // MAIN CALCULATIONS
+      // Remove mean from the data
+      if ( rm_mean == 3 )
+        n_ffts = max( 0, fix( (x_len-seg_len)/(seg_len-overlap) ) ) + 1;
+        x_len  = min( x_len, (seg_len-overlap)*(n_ffts-1)+seg_len );
+        if ( need_Pxx || need_Pxy )
+          x = x - sum( x(1:x_len) ) / x_len;
+        end
+        if ( arg2_is_y || need_Pxy)
+          y = y - sum( y(1:x_len) ) / x_len;
+        end
+      end
+      //
+      // Calculate and accumulate periodograms
+      //   xx and yy are padded data segments
+      //   Pxx, Pyy, Pyy are periodogram sums, Vxx is for confidence interval
+      xx = zeros(Nfft,1);
+      yy = xx;
+      Pxx = xx;
+      Pxy = xx;
+      Pyy = xx;
+      if ( conf>0 )
+        Vxx = xx;
+      else
+        Vxx = [];
+      end
+      n_ffts = 0;
+      for start_seg = [1:seg_len-overlap:x_len-seg_len+1]
+        end_seg = start_seg+seg_len-1;
+        // Don't truncate/remove the zero padding in xx and yy
+        if ( need_Pxx || need_Pxy )
+          if ( rm_mean==1 ) // remove mean from segment
+            xx(1:seg_len) = window .* ( ...
+              x(start_seg:end_seg) - sum(x(start_seg:end_seg)) / seg_len);
+          elseif ( rm_mean == 2 ) // remove linear trend from segment
+            xx(1:seg_len) = window .* detrend( x(start_seg:end_seg) );
+          else // rm_mean==0 or 3
+            xx(1:seg_len) = window .* x(start_seg:end_seg);
+          end
+          fft_x = fft1(xx);
+        end
+        if ( need_Pxy || need_Pyy )
+          if ( rm_mean==1 ) // remove mean from segment
+            yy(1:seg_len) = window .* ( ...
+              y(start_seg:end_seg) - sum(y(start_seg:end_seg)) / seg_len);
+          elseif ( rm_mean == 2 ) // remove linear trend from segment
+            yy(1:seg_len) = window .* detrend( y(start_seg:end_seg) );
+          else // rm_mean==0 or 3
+            yy(1:seg_len) = window .* y(start_seg:end_seg);
+          end
+          fft_y = fft1(yy);
+        end
+        if ( need_Pxx )
+          // force Pxx to be real; pgram = periodogram
+          pgram = real(fft_x .* conj(fft_x));
+          Pxx = Pxx + pgram;
+          // sum of squared periodograms is required for confidence interval
+          if ( conf>0 )
+            Vxx = Vxx + pgram .^2;
+          end
+        end
+        if ( need_Pxy )
+          // Pxy (cross power spectrum) is complex. Do not force to be real.
+          Pxy = Pxy + fft_y .* conj(fft_x);
+        end
+        if ( need_Pyy )
+          // force Pyy to be real
+          Pyy = Pyy + real(fft_y .* conj(fft_y));
+        end
+        n_ffts = n_ffts +1;
+      end
+      //
+      // Calculate confidence interval
+      //    -- incorrectly assumes that the periodogram has Gaussian probability
+      //       distribution (actually, it has a single-sided (e.g. exponential)
+      //       distribution.
+      // Sample variance of periodograms is (Vxx-Pxx.^2/n_ffts)/(n_ffts-1).
+      //    This method of calculating variance is more susceptible to round-off
+      //  error, but is quicker, and for double-precision arithmetic and the
+      //  inherently noisy periodogram (variance==mean^2), it should be OK.
+      if ( conf>0 && need_Pxx )
+        if ( n_ffts<2 )
+          Vxx = zeros(Nfft,1);
+        else
+          // Should use student distribution here (for unknown variance), but tinv
+          // is not a core Matlab function (is in statistics toolbox. Grrr)
+          Vxx = (erfinv(conf)*sqrt(2*n_ffts/(n_ffts-1))) * sqrt(Vxx-Pxx.^2/n_ffts);
+        end
+      end
+      //
+      // Convert two-sided spectra to one-sided spectra (if range == 0).
+      // For one-sided spectra, contributions from negative frequencies are added
+      // to the positive side of the spectrum -- but not at zero or Nyquist
+      // (half sampling) frequencies.  This keeps power equal in time and spectral
+      // domains, as required by Parseval theorem.
+      //
+      if (  ~ range  )
+        if (modulo(Nfft,2) == 0 )    // one-sided, Nfft is even
+          psd_len = Nfft/2+1;
+          if ( need_Pxx )
+            Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1); 0];
+            if ( conf>0 )
+              Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1); 0];
+            end
+          end
+          if ( need_Pxy )
+            Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1); 0]);
+          end
+          if ( need_Pyy )
+            Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1); 0];
+          end
+        else                    // one-sided, Nfft is odd
+          psd_len = (Nfft+1)/2;
+          if ( need_Pxx )
+            Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1)];
+            if ( conf>0 )
+              Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1)];
+            end
+          end
+          if ( need_Pxy )
+            Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1)]);
+          end
+          if ( need_Pyy )
+            Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1)];
+          end
+        end
+      else                      // two-sided (and shifted)
+        psd_len = Nfft;
+      end
+      // end MAIN CALCULATIONS
+      //
+      // SCALING AND OUTPUT
+      // Put all results in matrix, one row per spectrum
+      //   Pxx, Pxy, Pyy are sums of periodograms, so "n_ffts"
+      //   in the scale factor converts them into averages
+      spectra    = zeros(psd_len,n_results);
+      spect_type = zeros(n_results,1);
+      scale = n_ffts * seg_len * Fs * win_meansq;
+      if ( do_power )
+        spectra(:,do_power) = Pxx / scale;
+        spect_type(do_power) = 1;
+        if ( conf>0 )
+          Vxx = [Pxx-Vxx Pxx+Vxx]/scale;
+        end
+      end
+      if ( do_cross )
+        spectra(:,do_cross) = Pxy / scale;
+        spect_type(do_cross) = 2;
+      end
+      if ( do_trans )
+        spectra(:,do_trans) = Pxy ./ Pxx;
+        spect_type(do_trans) = 3;
+      end
+      if ( do_coher )
+        // force coherence to be real
+        spectra(:,do_coher) = real(Pxy .* conj(Pxy)) ./ Pxx ./ Pyy;
+        spect_type(do_coher) = 4;
+      end
+      if ( do_ypower )
+        spectra(:,do_ypower) = Pyy / scale;
+        spect_type(do_ypower) = 5;
+      end
+      freq = [0:psd_len-1].' * ( Fs / Nfft );
+      //
+      // range='shift': Shift zero-frequency to the middle
+      if ( range == 2 )
+        len2 = fix((Nfft+1)/2);
+        spectra = [ spectra(len2+1:Nfft,:); spectra(1:len2,:)];
+        freq    = [ freq(len2+1:Nfft)-Fs; freq(1:len2)];
+        if ( conf>0 )
+          Vxx = [ Vxx(len2+1:Nfft,:); Vxx(1:len2,:)];
+        end
+      end
+      //
+      //  RETURN RESULTS or PLOT
+      // droping non postive values for plots
+      function res= postives(x)
+            j = 1 ;res = [];
+            for i=1:size(x,2)
+                if or( x(:,i) < 0 )
+                    j = j - 1
+                    warning("warning: axis: omitting non-positive data in log plot")
+                else
+                    res(:,j) = x (:,i)
+                end
+                j = j + 1
+            end
+      endfunction
+      // --------------------
+      if ( nargout>=2 && conf>0 )
+        varargout(2) = Vxx;
+      end
+      if ( nargout>=(2+(conf>0)) )
+        // frequency is 2nd or 3rd return value,
+        // depends on if 2nd is confidence interval
+        varargout(2+(conf>0)) = freq;
+      end
+      if ( nargout>=1 )
+        varargout(1) = spectra;
+      else
+        //
+        // Plot the spectra if there are no return variables.
+        plot_title=['power spectrum x ';
+                    'cross spectrum   ';
+                    'transfer function';
+                    'coherence        ';
+                    'power spectrum y ' ];
+        for ii = 1: n_results
+          if ( conf>0 && spect_type(ii)==1 )
+            Vxxxx = Vxx;
+          else
+            Vxxxx = [];
+          end
+          if ( n_results > 1 )
+            figure();
+          end
+          if ( plot_type == 1 )
+            plot(freq,[abs(spectra(:,ii)) Vxxxx]);
+          elseif ( plot_type == 2 )
+            semilogx(freq,[abs(spectra(:,ii)) Vxxxx]);
+          elseif ( plot_type == 3 )
+            semilogy(freq,[abs(spectra(:,ii)) postives(Vxxxx)]);
+          elseif ( plot_type == 4 )
+            loglog(freq,[abs(spectra(:,ii)) Vxxxx]);
+          elseif ( plot_type == 5 )  // db
+            ylabel( 'amplitude (dB)' );
+            plot(freq,[10*log10(abs(spectra(:,ii))) 10*log10(abs(Vxxxx))]);
+          end
+          title( char(plot_title(spect_type(ii),:)) );
+          ylabel( 'amplitude' );
+          // Plot phase of cross spectrum and transfer function
+          if ( spect_type(ii)==2 || spect_type(ii)==3 )
+            figure();
+            if ( plot_type==2 || plot_type==4 )
+              semilogx(freq,180/%pi*angle(spectra(:,ii)));
+            else
+              plot(freq,180/%pi*angle(spectra(:,ii)));
+            end
+            title( char(plot_title(spect_type(ii),:)) );
+            ylabel( 'phase' );
+          end
+        end
+      end
     end
-end
-
-endfunction
+  endfunction 
+  /*
+  
+// demo 1  // use "hold on" and "hold off" for octave...
+   pi = %pi ; i = %i;
+   Fs = 25;
+   a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
+   white = rand(1,16384);
+   signal = detrend(filter(0.70181,a,white));
+   skewed = signal.*exp(2*pi*i*2/25*[1:16384]);
+   compat = pwelch ([]);
+   hold on;
+   pwelch(skewed,[],[],[],Fs,'shift','semilogy');
+   pwelch(skewed,[],[],[],Fs,0.95,'shift','semilogy');
+   pwelch('R12+');
+   pwelch(signal,'squared');
+   pwelch (compat);
+   hold off;
+// demo 2 // use "hold on" and "hold off" for octave...
+ a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
+ white = rand(1,16384);
+ signal = detrend(filter(0.70181,a,white));
+ skewed = signal.*exp(2*pi*i*2/25*[1:16384]);
+ compat = pwelch ([]);
+ hold on;
+ pwelch(signal);
+ pwelch(skewed);
+ pwelch(signal,'shift','semilogy');
+ pwelch (compat);
+ hold off
+// demo 3 
+ a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
+ white = rand(1,16384);
+ signal = detrend(filter(0.70181,a,white));
+ compat = pwelch ([]);
+ pwelch(signal,3640,[],4096,2*pi,[],'no-strip');
+ pwelch (compat);  
+// demo 4 // use "hold on" and "hold off" for octave...
+ a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
+ white = rand(1,16384);
+ signal = detrend(filter(0.70181,a,white));
+ compat = pwelch ([]);
+ hold on;
+ pwelch(signal,[],[],[],2*pi,0.95,'no-strip');
+ pwelch(signal,64,[],[],2*pi,'no-strip');
+ pwelch(signal,64,[],256,2*pi,'no-strip');
+ pwelch (compat);
+ hold off;
+ //demo 5
+ a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
+ white = rand(1,16384);
+ signal = detrend(filter(0.70181,a,white));
+ compat = pwelch ('psd');
+ pwelch(signal,'squared');
+ pwelch({});
+ pwelch(white,signal,'trans','coher','short')
+ pwelch (compat);
+ */
diff --git a/macros/rceps.sci b/macros/rceps.sci
index f0433a3..831611f 100644
--- a/macros/rceps.sci
+++ b/macros/rceps.sci
@@ -1,48 +1,48 @@
+ /*Description
+        Real cepstrum and minimum-phase reconstruction
+        If  called with two output arguments,
+        the minimum phase reconstruction of the signal x is returned in ym.Calling Sequence
+        [y,ym] = rceps(x)
+    Parameters
+        x : A vector or a Matirx
+        y : Real cepstrum
+        ym : Minimum-phase reconstruction
+    Dependencies fft1 , ifft1
+    Example:
+        // create a 45 Hz sine wave sampled at 100 Hz. 
+        t = 0:0.01:1.27;
+        s1 = sin(2*%pi*45*t);
+        s2 = s1 + 0.5*[zeros(1,20) s1(1:108)]; //Add an echo of the signal, with half the amplitude, 0.2 seconds after the beginning of the signal.
+        c = rceps(s2); // real cepstrum of signal 
+        plot(t,c)
+*/
 function [y, xm]= rceps(x)
-//Produce the cepstrum of the signal x, and if desired, the minimum phase reconstruction of the signal x.
-//Calling Sequence
-//[y, xm] = rceps(x)
-//Parameters 
-//x: real or complex vector input
-//Produce the cepstrum of the signal x, and if desired, the minimum phase reconstruction of the signal x. If x is a matrix, do so for each column of the matrix.
-//Examples
-// f0 = 70; Fs = 10000;                   # 100 Hz fundamental, 10kHz sampling rate
-// a = poly (0.985 * exp (1i*pi*[0.1, -0.1, 0.3, -0.3])); # two formants
-// s = 0.005 * randn (1024, 1);           # Noise excitation signal
-// 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
-
-if(argn(2)~= 1 )
-error("Wrong number of Input Arguments");
-end
-
-if(argn(1)>2)
-error("Wrong number of Output Arguments")
-end
-  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)];
+    if(argn(2)~= 1 )
+    error("Wrong number of Input Arguments");
+    end
+    if(argn(1)>2)
+    error("Wrong number of Output Arguments")
+    end
+      f = abs(fft1(x));
+      if (or(f == 0))
+        error ("The spectrum of x contains zeros, unable to compute real cepstrum");
       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))];
+      y = real(ifft1(log(f)));
+      if argn(1) == 2 then
+        n=max(size(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
-    end
-    xm = real(ifft1(exp(fft1(xm)')));
-  end
-  
 endfunction
diff --git a/macros/shanwavf.sci b/macros/shanwavf.sci
index a597214..e088f60 100644
--- a/macros/shanwavf.sci
+++ b/macros/shanwavf.sci
@@ -1,30 +1,39 @@
+/*Description
+    Compute the Complex Shannon wavelet.
+    The complex Shannon wavelet is defined by a bandwidth parameter fb, a wavelet center frequency fc, and the expression
+    psi(x) = f * b^{1/2}sinc(fb . x) e^{2 pi i f c x}
+    on an n-point regular grid in the interval of lb to ub.
+Calling Sequence
+    [psi, x]= shanwavf(lb, ub, n, fb, fc)
+Input  Parameters
+    lb, ub (Real valued scalers) : Interval endpoints lb ≤ ub, specified as a pair of real-valued scalars.
+    n (Real valued integer strictly positive)` : Number of regularly spaced points in the interval [lb,ub], specified as a positive integer.
+    fb	: Time-decay parameter of the wavelet (bandwidth in the frequency domain). Must be a positive scalar.
+    fc	: Center frequency of the complex Shannon wavelet, specified as a positive scalar.
+Output Parameters
+     psi : Complex Shannon wavelet evaluated on the n point regular grid x in the interval [lb,ub], returned as a 1-by-n vector.
+     x   : Grid where the complex Shannon wavelet is evaluated, returned as a 1-by-n vector. The sample points are evenly distributed between lb and ub.
+Examples
+      1.[a,b]=shanwavf (2,8,3,1,6)
+           a =   [-3.8982e-17 + 1.1457e-31i   3.8982e-17 - 8.4040e-31i  -3.8982e-17 + 4.5829e-31i]
+           b =   [2   5   8]   */
 function [psi,x]=shanwavf(lb,ub,n,fb,fc)
-
-// Complex Shannon Wavelet
-// Calling Sequence
-//	[psi,x]=shanwavf(lb,ub,n,fb,fc)
-// Parameters
-//	lb: Real or complex valued vector or matrix
-//	ub: Real or complex valued vector or matrix
-//	n: Real valued integer strictly positive
-//	fb: Real or complex valued vector or matrix, strictly positive value for scalar input
-//	fc: Real or complex valued vector or matrix, strictly positive value for scalar input
-// Description
-//	This is an Octave function
-//	This function implements the complex Shannon wavelet function and returns the value obtained. The complex Shannon wavelet is defined by a bandwidth parameter FB, a wavelet center frequency FC on an N point regular grid in the interval [LB,UB].
-// Examples
-// 1.	[a,b]=shanwavf (2,8,3,1,6)
-//	a =   [-3.8982e-17 + 1.1457e-31i   3.8982e-17 - 8.4040e-31i  -3.8982e-17 + 4.5829e-31i]
-//	b =   [2   5   8]
-// 2.	[a,b]=shanwavf(1,2,1,[2,2;i,2],[-1,2;-i,i])
-//	a =   [-5.5128e-17 - 2.7005e-32i  -5.5128e-17 + 5.4010e-32i;
-//	       8.6404e+06 + 8.6404e+06i  -1.9225e-22 - 0.0000e+00i]
-//	b =  2
-
-funcprot(0);
-rhs=argn(2);
-if (rhs~=5) then
-	error ("Wrong number of input arguments.")
-else [psi,x]=callOctave("shanwavf",lb,ub,n,fb,fc)
-end
+    funcprot(0);
+    rhs=argn(2);
+    if (rhs~=5) then
+        error ("Wrong number of input arguments.")
+    else 
+        if (n <= 0 || floor(n) ~= n)
+            error("n must be an integer strictly positive");
+        elseif (fc <= 0 || fb <= 0)
+            error("fc and fb must be strictly positive");
+        end
+        x = linspace(lb,ub,n);
+        sincx=x;
+        for i=1:n
+            sincx(i)=sin(fb*x(i)*%pi)/(fb*x(i)*%pi);
+        end    
+        psi = (fb.^0.5).*(sincx.*exp(2.*%i.*%pi.*fc.*x));
+    end
 endfunction
+
diff --git a/macros/spectral_adf.sci b/macros/spectral_adf.sci
index 38b3990..e1cc5a2 100644
--- a/macros/spectral_adf.sci
+++ b/macros/spectral_adf.sci
@@ -1,33 +1,67 @@
-function y= spectral_adf(x, varargin)
-// Return the spectral density estimator given a vector of autocovariances C, window name WIN, and bandwidth, B.
-//Calling Sequence
-//spectral_adf(C)
-//spectral_adf(C, WIN)
-//spectral_adf(C, WIN, B)
-//Parameters 
-//C: Autocovariances
-//WIN: Window names
-//B: Bandwidth 
-//Description
-//Return the spectral density estimator given a vector ofautocovariances C, window name WIN, and bandwidth, B.
-//The window name, e.g., "triangle" or "rectangle" is used to search for a function called 'WIN_lw'.
-//If WIN is omitted, the triangle window is used.
-//If B is omitted, '1 / sqrt (length (C))' is used.
-
-   
-funcprot(0);
-rhs= argn(2);
-if(rhs<1 | rhs>3)
-error("Wrong number of input arguments")
-end
-
-select(rhs)
-	case 1 then
-		y= callOctave("spectral_adf", x);
-	case 2 then
-		y= callOctave("spectral_adf", x , varargin(1));
-	case 3 then
-		y= callOctave("spectral_adf", x , varargin(1), varargin(2));
-
-end
-endfunction
+/* spectral_adf  
+  Calling Sequence
+      spectral_adf (c)
+      spectral_adf (c, win)
+      spectral_adf (c, win, b)
+  Parameters
+      c : Vector of autocovariances
+      win : The window name . Default window is "triangle"
+      b : Bandwidth. Default is 1/sqrt(length(c))
+  Description
+      Return the spectral density estimator given a vector of autocovariances c, window name win, and bandwidth, b.
+      The window name, e.g., "triangle" or "rectangle" is used to search for a function called win_lw.
+      If win is omitted, the triangle window is used.
+      If b is omitted, 1 / sqrt (length (c)) is used.
+  Dependencies: fft1  */
+function sde = spectral_adf (c, win, b)
+    //c should be a vector
+    if ~isvector(c)
+      error("spectral_adf: input c should be a vector")
+    end
+    // maximum along the columns
+    function max_value = col_max (A)
+      max_value = zeros ( 1, size(A,2))
+      for i = 1:length(max_value)
+          max_value(i) = max(A(:,i))
+      end
+    endfunction
+    //window functions
+    function c = triangle_lw (n, b)
+      c = 1 - (0 : n-1)' * b;
+      c = [c' ; zeros(1, n)]
+      c = col_max (c)'
+    endfunction
+    function c = rectangle_lw (n, b)
+      c = zeros (n, 1);
+      t = floor (1 / b);
+      c(1:t) = 1;
+    endfunction
+    // main part
+    nargin = argn(2)
+    if (nargin < 1)
+      error("wrong number of inputs.");
+    end
+    cr = max(size (c) );
+    if (size (c , 2) > 1)
+      c = c';
+    end
+    if (nargin < 3)
+      b = 1 / ceil (sqrt (cr));
+    end
+    if (nargin == 1)
+      w = triangle_lw (cr, b);
+    elseif (~ (type (win) == 10) )
+      error ("spectral_adf: WIN must be a string");
+    elseif (~strcmp (win , "rectangle" ) )
+      w = rectangle_lw(cr , b) ;
+    elseif (~strcmp (win , "triangle" ) ) 
+        w = triangle_lw(cr , b) ;
+    else 
+        error("Invalid window or this window is not supported yet")
+    end 
+    c=c .* w;
+    sde = 2 * real (fft1 (c)) - c(1);
+    zer= zeros(cr, 1);
+    sde = [zer sde];
+    sde(:, 1) = (0 : cr-1)' / cr;
+  endfunction
diff --git a/macros/spectral_xdf.sci b/macros/spectral_xdf.sci
index f0d457b..6081c71 100644
--- a/macros/spectral_xdf.sci
+++ b/macros/spectral_xdf.sci
@@ -1,31 +1,66 @@
-function y= spectral_xdf(x, varargin)
-// Return the spectral density estimator given a data vector X, window name WIN, and bandwidth, B.
-//Calling Sequence
-//spectral_xdf(X)
-//spectral_xdf(X, WIN)
-//spectral_xdf(X, WIN, B)
-//Parameters 
-//X: Data Vector
-//WIN: Window names
-//B: Bandwidth 
-//Description
-//Return the spectral density estimator given a data vector X, window name WIN, and bandwidth, B.
-//The window name, e.g., "triangle" or "rectangle" is used to search for a function called 'WIN_lw'.
-//If WIN is omitted, the triangle window is used.
-//If B is omitted, '1 / sqrt (length (X))' is used.
-funcprot(0);
-rhs= argn(2);
-if(rhs<1 | rhs>3)
-error("Wrong number of input arguments")
-end
-
-select(rhs)
-	case 1 then
-		y= callOctave("spectral_xdf", x);
-	case 2 then
-		y= callOctave("spectral_xdf", x , varargin(1));
-	case 3 then
-		y= callOctave("spectral_xdf", x , varargin(1), varargin(2));
-
-end
+/* Description
+        Return the spectral density estimator given a data vector x, window name win, and bandwidth, b.
+        The window name, e.g., "triangle" or "rectangle" is used to search for a function called win_sw.
+        If win is omitted, the triangle window is used.
+        If b is omitted, 1 / sqrt (length (x)) is used.
+    Calling Sequence
+        spectral_xdf (x)
+        spectral_xdf (x, win)
+        spectral_xdf (x, win, b)
+    Parameters
+        x : Data vector
+        win : the window name . Default "triangle" is used .
+        b : Bandwidth . Default value 1/sqrt(length(x))
+Dependencies:fft1 ifft1 */
+function sde = spectral_xdf (x, win, b)
+    // check x is a vector or not
+    if ~isvector(x) 
+        error("spectral_xdf :  x must a data vector")
+    end
+    // available windows - rectangle
+    function c = rectangle_sw (n, b)
+          c = zeros (n, 1);
+          c(1) = 2 / b + 1;
+          l = (2:n)' - 1;
+          l = 2 * %pi * l / n;
+          c(2:n) = sin ((2/b + 1) * l / 2) ./ sin (l / 2);
+    endfunction
+    // triangle
+    function retval = triangle_sw (n, b)
+          retval = zeros (n,1);
+          retval(1) = 1 / b;
+          l = (2:n)' - 1;
+          l = 2 * %pi * l / n;
+          retval(2:n) = b * (sin (l / (2*b)) ./ sin (l / 2)).^2;
+    endfunction
+    // main function
+    nargin = argn(2)
+    if (nargin < 1)
+      error("invalid inputs");
+    end
+    xr = max(size (x) );
+    if (size (x, 2) > 1)
+      x = x';
+    end
+    if (nargin < 3)
+      b = 1 / ceil (sqrt (xr));
+    end
+    if (nargin == 1)
+      w = triangle_sw (xr, b);
+    elseif (~ (type(win)== 10))
+      error ("spectral_xdf: WIN must be a string");
+    elseif (~strcmp (win , "triangle" ) ) 
+        w = triangle_sw (xr , b);
+    elseif (~strcmp (win , "rectangle" ) )     
+        w = rectangle_sw ( xr , b);
+    else
+      error("Invalid window or this window is not supported yet");
+    end
+    x = x - sum (x) / xr;
+    sde = (abs (fft1 (x)) / xr).^2;
+    sde = real (ifft1 (fft1 (sde) .* fft1 (w)));
+    zer = zeros(xr,1);
+    sde = [zer sde];
+    sde(:, 1) = (0 : xr-1)' / xr;
+    sde = clean(sde);
 endfunction
diff --git a/macros/tfe.sci b/macros/tfe.sci
index 5acc391..5a6cf5a 100644
--- a/macros/tfe.sci
+++ b/macros/tfe.sci
@@ -1,27 +1,41 @@
-function [Pxx,freqs] = tfe(x,y,Nfft,Fs,win,overlap,ran,plot_type,detrends)
-//Estimate transfer function of system with input "x" and output "y". Use the Welch (1967) periodogram/FFT method.
-//Calling Sequence
-// [Pxx,freq] = tfe(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
-//Parameters 
-//x: [non-empty vector] system-input time-series data
-//y: [non-empty vector] system-output time-series data
-//win:[real vector] of window-function values between 0 and 1; the data segment has the same length as the window. Default window shape is Hamming. [integer scalar] length of each data segment.  The default value is window=sqrt(length(x)) rounded up to the nearest integer power of 2; see 'sloppy' argument.
-//overlap:[real scalar] segment overlap expressed as a multiple of window or segment length.   0 <= overlap < 1, The default is overlap=0.5 .
-//Nfft:[integer scalar] Length of FFT.  The default is the length of the "window" vector or has the same value as the scalar "window" argument.  If Nfft is larger than the segment length, "seg_len", the data segment is padded with "Nfft-seg_len" zeros.  The default is no padding. Nfft values smaller than the length of the data segment (or window) are ignored silently.
-//Fs:[real scalar] sampling frequency (Hertz); default=1.0 
-//range:'half',  'onesided' : frequency range of the spectrum is zero up to but not including Fs/2.  Power from negative frequencies is added to the positive side of the spectrum, but not at zero or Nyquist (Fs/2) frequencies.  This keeps power equal in time and spectral domains.  See reference [2]. 'whole', 'twosided' : frequency range of the spectrum is-Fs/2 to Fs/2, with negative frequenciesstored in "wrap around" order after the positivefrequencies; e.g. frequencies for a 10-point 'twosided'spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1 'shift', 'centerdc' : same as 'whole' but with the first half of the spectrum swapped with second half to put the zero-frequency value in the middle. (See "help fftshift". If data (x and y) are real, the default range is 'half', otherwise default range is 'whole'.
-//plot_type: 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db': specifies the type of plot.  The default is 'plot', which means linear-linear axes. 'squared' is the same as 'plot'. 'dB' plots "10*log10(psd)".  This argument is ignored and a spectrum is not plotted if the caller requires a returned value.
-//detrends:'no-strip', 'none' -- do NOT remove mean value from the data'short', 'mean' -- remove the mean value of each segment from each segment of the data. 'linear',-- remove linear trend from each segment of the data.'long-mean'-- remove the mean value from the data before splitting it into segments. This is the default.
-//Description
-//Estimate transfer function of system with input "x" and output "y". Use the Welch (1967) periodogram/FFT method.
-	funcprot(0);
-	rhs= argn(2);
-	lhs= argn(1);
-	if(rhs < 10 | rhs > 10)
-		error("Wrong number of input arguments");
-	end
-	select(rhs)
-	case 10 then
-		[Pxx,freqs] = callOctave("tfe",x,y,Nfft,Fs,win,overlap,ran,plot_type,detrends);
-	end
-endfunction
-\ No newline at end of file
+/*
+Dependency : pwelch
+Calling Sequence:
+    [Pxx,freq] = tfe(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
+        Estimate transfer function of system with input "x" and output "y".
+        Use the Welch (1967) periodogram/FFT method.
+        Compatible with Matlab R11 tfe and earlier.
+        See "help pwelch" for description of arguments, hints and references — especially hint (7) for Matlab R11 defaults.*/
+function varargout = tfe(varargin)
+    nargout = argn (1)
+    nargin = argn(2)
+    // Check fixed argument
+    if ( nargin<2 )
+      error( 'tfe: Need at least 2 args. Use help tfe.' );
+    end
+    nvarargin = max(size(varargin));
+    // remove any pwelch RESULT args and add 'trans'
+    for iarg=1:nvarargin
+      arg = varargin(iarg);
+      if ( ~isempty(arg) && type(arg) == [10 ] && ( ~strcmp(arg,'power') || ...
+             ~strcmp(arg,'cross') || ~strcmp(arg,'trans') || ...
+             ~strcmp(arg,'coher') || ~strcmp(arg,'ypower') ))
+        varargin(iarg) = [];
+      end
+    end
+    varargin(nvarargin+1) = 'trans';
+    saved_compatib = pwelch('R11-');
+    if ( nargout==0 )
+      pwelch(varargin(:));
+    elseif ( nargout==1 )
+      Pxx = pwelch(varargin(:));
+      varargout(1) = Pxx;
+    elseif ( nargout>=2 )
+      [Pxx,f] = pwelch(varargin(:));
+      varargout(1) = Pxx;
+      varargout(2) = f;
+    end
+    pwelch(saved_compatib);
+  endfunction
+
+
diff --git a/macros/unwrap2.sci b/macros/unwrap2.sci
index b8ea9de..7844e38 100644
--- a/macros/unwrap2.sci
+++ b/macros/unwrap2.sci
@@ -1,35 +1,211 @@
-function Y = unwrap2 (X, TOL, DIM)
-//Unwrap radian phases by adding or subtracting multiples of 2*pi.
-//Calling Sequence
-//B = unwrap(X)
-//B = unwrap(X, TOL)
-//B = unwrap(X, TOL, DIM)
-//Parameters
-//Description
-//This function unwraps radian phases by adding or subtracting multiples of 2*pi as appropriate to remove jumps greater than TOL.
-//
-//    TOL defaults to pi.
-//
-//Unwrap will work along the dimension DIM.  If DIM is unspecified it defaults to the first non-singleton dimension.
-//Examples
-//unwrap2([1,2,3])
-//ans = 
-//        1.    2.    3.  
-funcprot(0);
-lhs = argn(1);
-rhs = argn(2);
-if (rhs < 1 | rhs > 3)
-error("Wrong number of input arguments.");
-end
-
-select(rhs)
-	
-	case 1 then
-        Y = callOctave("unwrap",X);
-	case 2 then
-	Y = callOctave("unwrap",X,TOL);
-        case 3 then
-        Y = callOctave("unwrap",X,TOL,DIM);  
-        end
-
+/*Description:
+        The unwrap function adjusts radian phases in the input array x by adding or subtracting multiples of
+        2π as necessary to remove phase jumps that exceed the specified tolerance tol. If tol is not provided, it defaults to 𝜋
+        Radian Phases: These are typically angles or phases expressed in radians, commonly encountered in signal processing and communication systems.
+        Tolerance (tol): Determines the maximum allowable discontinuity in the phases.
+        If the difference between consecutive phases exceeds tol, unwrap adjusts the phase by adding or subtracting 2π.
+        Dimension (dim): Specifies the dimension along which the unwrapping operation is applied.
+        By default, unwrap operates along the first non-singleton dimension of the input array x.
+Calling Sequence:
+        b = unwrap(x)
+        b = unwrap(x, tol)
+        b = unwrap(x, tol, dim)
+Parameters:
+        x: Input array containing radian phases to be unwrapped.
+        tol (optional): Tolerance parameter specifying the maximum jump allowed between consecutive phases before adding or subtracting 2π. Defaults to 𝜋
+        dim (optional): Dimension along which to unwrap the phases. If unspecified, dim defaults to the first non-singleton dimension of the array x.
+Dependencies : ipermute*/
+function retval = unwrap2 (x, tol, dim)
+  nargin = argn(2)
+  if (nargin < 1)
+    error("invalid number of inputs");
+  end
+  if (~ (type(x) == [ 1 5 8]) || or(type(x)==[4,6]))
+    error ("unwrap2: X must be numeric");
+  end
+  if (nargin < 2 || isempty (tol))
+    tol = %pi;
+  end
+  // Don't let anyone use a negative value for TOL.
+  tol = abs (tol);
+  nd = ndims (x);
+  sz = size (x);
+  if (nargin == 3)
+    if (~(or(type(dim)==[1 5 8])&& isscalar (dim) && ...
+            dim == fix (dim)) || ~(1 <= dim))
+      error ("unwrap2: DIM must be an integer and a valid dimension");
+    end
+  else
+    // Find the first non-singleton dimension.
+    dim = find (sz > 1, 1)
+    if isempty(dim)
+      dim = 1;
+    end   
+  end
+  rng = 2*%pi;
+  // Handle case where we are trying to unwrap a scalar, or only have
+  // one sample in the specified dimension (a given when dim > nd).
+  if ((dim > nd) || ( sz(dim) == 1))
+    retval = x;
+    return;
+  end
+  if (and(abs(x(:))<%inf ) )
+    // Take first order difference so that wraps will show up as large values
+    // and the sign will show direction.
+    if length(sz) < 3 
+      sz(3) = 1 ;
+    end
+    sz(dim) = 1;
+    zero_padding = zeros (sz(1),sz(2),sz(3));
+    d = cat (dim, zero_padding, -diff (x, 1, dim));
+    // Find only the peaks and multiply them by the appropriate amount
+    // of ranges so that there are kronecker deltas at each wrap point
+    // multiplied by the appropriate amount of range values.
+    p = round (abs (d)./rng) .* rng .* (double(double(d > tol) > 0) - double(double(d < -tol) > 0));
+    // Integrate this so that the deltas become cumulative steps to shift
+    // the original data.
+    retval = cumsum (p, dim) + x;
+  else
+    // Unwrap needs to skip over NaN, NA, Inf in wrapping calculations.
+    if (isvector (x))
+      // Simlpified path for vector inputs.
+      retval = x;
+      xfin_idx = abs(x)<%inf ;
+      xfin = x(xfin_idx);
+      d = cat (dim, 0, -diff(xfin, 1, dim));
+      p = round (abs (d)./rng) .* rng .* (double(double(d > tol) > 0) - double(double(d < -tol) > 0));
+      retval(xfin_idx) = xfin + cumsum (p, dim);
+    else
+      // For n-dimensional arrays with a possibly unequal number of non-finite
+      // values, mask entries with values that do not impact calcualation.
+            // Locate nonfinite values.
+      nf_idx = ~ abs(x)<%inf;
+      if (and(nf_idx(:)))
+        // Trivial case, all non-finite values
+        retval = x;
+        return;
+      end
+      // Permute all operations to occur along dim 1.  Inverse permute at end.
+      permuteflag = dim ~= 1;
+      if (permuteflag)
+        perm_idx = [1 : nd];
+        perm_idx([1, dim]) = [dim, 1];
+        x = permute (x, perm_idx);
+        nf_idx = permute (nf_idx, perm_idx);
+        sz([1, dim]) = sz([dim, 1]);
+        dim = 1;
+      end
+      // Substitute next value in dim direction for nonfinite values(ignoring
+      // any at trailing end) to prevent calculation impact.
+      x_nf = x(nf_idx); // Store nonfinite values.
+      zero_padding = zeros ([1, sz(2:$)]);
+      x = __fill_nonfinite_columnwise__ (x, nf_idx, zero_padding, sz, nd);
+      d = [zero_padding; -diff(x, 1, 1)];
+      p = round (abs (d)./rng) .* rng .* ...
+          (((d > tol) > 0) - ((d < -tol) > 0));
+      retval = x + cumsum (p, 1);
+      // Restore nonfinite values.
+      retval(nf_idx) = x_nf;
+      // Invert permutation.
+      if (permuteflag)
+        retval = ipermute (retval, perm_idx);
+      end
+    end
+  end
 endfunction
+function y = repelems(x,r)
+  y = [];
+  for i = 1:size(r,2)
+      y = [y, x(r(1,i)*ones(1, r(2,i)))];
+  end
+endfunction
+function x = __fill_nonfinite_columnwise__ (x, nonfinite_loc, zero_padding, szx, ndx)
+  // Replace non-finite values of x, as indicated by logical index
+  // nonfinite_loc, with next values.
+  flip_idx(1:ndx) = {':'};
+  flip_idx(1) = {szx(1):-1:1};
+  // Isolate nf values by location:
+  nf_front = cumprod (nonfinite_loc, 1);
+  nf_back = cumprod (nonfinite_loc(flip_idx{:}), 1)(flip_idx{:});
+  nf_middle = nonfinite_loc & ~ (nf_back | nf_front);
+  // Process bound/middle elements
+  locs_before = [diff(nf_middle, 1, 1); zero_padding] == 1;
+  locs_after = diff ([zero_padding; nf_middle], 1, 1) == -1;
+  mid_gap_sizes = find (locs_after) - find (locs_before) - 1;
+  x(nf_middle) = repelems (x(locs_after), ...
+                          [1 : numel(mid_gap_sizes); mid_gap_sizes'])';
+  // Process front edge elements
+  nf_front = nf_front & ~ and (nonfinite_loc, 1); // Remove all nf columns.
+  locs_after = diff ([zero_padding; nf_front], 1, 1) == -1;
+  front_gap_sizes = (sum (nf_front, 1))(any (nf_front, 1))(:);
+  x(nf_front) = repelems (x(locs_after), ...
+                             [1:length(front_gap_sizes); front_gap_sizes'])';
+endfunction
+/*
+//Test cases
+i = 0;
+t = [];
+r = [0:100];                         // original vector
+w = r - 2*%pi*floor ((r+%pi)/(2*%pi));  // wrapped into [-pi,pi]
+tol = 1e3*%eps;
+assert_checkalmostequal (r,  unwrap2 (w),  tol)
+assert_checkalmostequal (r', unwrap2 (w'), tol)
+assert_checkalmostequal ([r',r'], unwrap2 ([w',w']), tol)
+assert_checkalmostequal ([r; r ], unwrap2 ([w; w ], [], 2), tol)
+assert_checkalmostequal(r + 10, unwrap2 (10 + w), tol)
+assert_checkequal (w', unwrap2 (w', [], 2))
+assert_checkequal(w,  unwrap2 (w,  [], 1))
+assert_checkequal([w; w], unwrap2 ([w; w]))
+// Test that small values of tol have the same effect as tol = pi
+assert_checkalmostequal(r, unwrap2 (w, 0.1), tol)
+assert_checkalmostequal(r, unwrap2 (w, %eps), tol)
+// Test that phase changes larger than 2*pi unwrap properly
+assert_checkalmostequal([0;  1],        unwrap2([0;  1]))
+assert_checkalmostequal([0;  4 - 2*%pi], unwrap2 ([0;  4]))
+assert_checkalmostequal([0;  7 - 2*%pi], unwrap2 ([0;  7]))
+assert_checkalmostequal([0; 10 - 4*%pi], unwrap2 ([0; 10]))
+assert_checkalmostequal([0; 13 - 4*%pi], unwrap2 ([0; 13]))
+assert_checkalmostequal([0; 16 - 6*%pi], unwrap2 ([0; 16]))
+assert_checkalmostequal([0; 19 - 6*%pi], unwrap2 ([0; 19]))
+//test
+A = [%pi*(-4), %pi*(-2+1/6), %pi/4, %pi*(2+1/3), %pi*(4+1/2), %pi*(8+2/3), %pi*(16+1), %pi*(32+3/2), %pi*64];
+assert_checkalmostequal (unwrap2 (A), unwrap2 (A, %pi));
+assert_checkalmostequal (unwrap2 (A, %pi), unwrap2 (A, %pi, 2));
+assert_checkalmostequal (unwrap2 (A', %pi), unwrap2 (A', %pi, 1));
+//test
+A = [%pi*(-4); %pi*(2+1/3); %pi*(16+1)];
+B = [%pi*(-2+1/6); %pi*(4+1/2); %pi*(32+3/2)];
+C = [%pi/4; %pi*(8+2/3); %pi*64];
+D = [%pi*(-2+1/6); %pi*(2+1/3); %pi*(8+2/3)];
+E(:, :, 1) = [A, B, C, D];
+E(:, :, 2) = [A+B, B+C, C+D, D+A];
+F(:, :, 1) = [unwrap2(A), unwrap2(B), unwrap2(C), unwrap2(D)];
+F(:, :, 2) = [unwrap2(A+B), unwrap2(B+C), unwrap2(C+D), unwrap2(D+A)];
+assert_checkalmostequal (unwrap2 (E), F);
+// Test trivial return for m = 1 and dim > nd
+assert_checkalmostequal (unwrap2 (ones(4,1), [], 1), ones(4,1))
+assert_checkalmostequal (unwrap2 (ones(4,1), [], 2), ones(4,1))
+assert_checkalmostequal (unwrap2 (ones(4,1), [], 3), ones(4,1))
+assert_checkalmostequal (unwrap2 (ones(4,3,2), [], 99), ones(4,3,2))
+// Test empty input return
+assert_checkalmostequal (unwrap2 ([]), [])
+assert_checkalmostequal (unwrap2 (ones (1,0)), ones (1,0))
+assert_checkalmostequal (unwrap2 (ones (1,0), [], 1), ones (1,0))
+assert_checkalmostequal (unwrap2 (ones (1,0), [], 2), ones (1,0))
+assert_checkalmostequal (unwrap2 (ones (1,0), [], 3), ones (1,0))
+// Test trivial return for m = 1 and dim > nd
+assert_checkalmostequal (unwrap2 (ones(4,1), [], 1), ones(4,1))
+assert_checkalmostequal (unwrap2 (ones(4,1), [], 2), ones(4,1))
+assert_checkalmostequal (unwrap2 (ones(4,1), [], 3), ones(4,1))
+assert_checkalmostequal (unwrap2 (ones(4,3,2), [], 99), ones(4,3,2))
+// Test empty input return
+assert_checkalmostequal (unwrap2 ([]), [])
+assert_checkalmostequal (unwrap2 (ones (1,0)), ones (1,0))
+assert_checkalmostequal (unwrap2 (ones (1,0), [], 1), ones (1,0))
+assert_checkalmostequal (unwrap2 (ones (1,0), [], 2), ones (1,0))
+assert_checkalmostequal (unwrap2 (ones (1,0), [], 3), ones (1,0))
+// Test handling of non-finite values
+x = %pi * [-%inf, 0.5, -1, %nan, %inf, -0.5, 1];
+assert_checkalmostequal (unwrap2 (x), %pi * [-%inf, 0.5, 1, %nan, %inf, 1.5, 1], %eps)
+assert_checkalmostequal (unwrap2 (x.'), %pi * [-%inf, 0.5, 1, %nan, %inf, 1.5, 1].', %eps)
+*/
+\ No newline at end of file
diff --git a/macros/xcorr2.sci b/macros/xcorr2.sci
index 688dc3a..40d3b57 100644
--- a/macros/xcorr2.sci
+++ b/macros/xcorr2.sci
@@ -1,34 +1,70 @@
-function c = xcorr2 (a, b, biasflag)
-//
-//Calling Sequence
-//c = xcorr2 (a)
-//c = xcorr2 (a, b)
-//c = xcorr2 (a, b, biasflag)
-//Parameters 
-//a:
-//b:
-//biasflag:
-//Description
-//This is an Octave function.
-
-//Examples
-//xcorr2(5,0.8,'coeff')
-//ans =  1
-
-funcprot(0);
-
-rhs = argn(2)
-if(rhs<1 | rhs>3)
-error("Wrong number of input arguments.");
-end
-
-	select(rhs)
-	case 1 then
-	c = callOctave("xcorr2",a);
-	case 2 then
-	c = callOctave("xcorr2",a,b);
-	case 3 then
-	c = callOctave("xcorr2",a,b,biasflag);
-	end;
+/*Calling Sequence
+        c = xcorr2 (a)
+        c = xcorr2 (a, b)
+        c = xcorr2 (a, b, scale)
+Description:
+        Compute the 2D cross-correlation of matrices a and b.
+        If b is not specified, computes autocorrelation of a, i.e., same as xcorr (a, a).
+        The optional argument scale, defines the type of scaling applied to the cross-correlation matrix. Possible values are:
+        "none" (default)
+            No scaling.
+        "biased"
+            Scales the raw cross-correlation by the maximum number of elements of a and b involved in the generation of any element of c.
+        "unbiased"
+            Scales the raw correlation by dividing each element in the cross-correlation matrix by the number of products a and b used to generate that element.
+        "coeff"
+            Scales the normalized cross-correlation on the range of [0 1] so that a value of 1 corresponds to a correlation coefficient of 1.
+  Examples
+        xcorr2(5,0.8,'coeff')
+            ans =  1 */
+function c = xcorr2 (a, b, scale)
+  funcprot(0);
+  nargin=argn(2);
+  if nargin < 3 then
+    scale = "none"
+  end
+  if (nargin < 1 || nargin > 3)
+    error("Wrong number of inputs")
+  elseif (nargin == 2 && type (b) == 10 )
+    scale = b;
+    b        = a;
+  elseif (nargin == 1)
+    // we have to set this case here instead of the function line, because if
+    // someone calls the function with zero argument, since a is never set, we
+    // will fail with "`a' undefined" error rather that print_usage
+    b = a;
+  end
+  if (ndims (a) ~= 2 || ndims (b) ~= 2)
+    error ("xcorr2: input matrices must must have only 2 dimensions");
+  end
+  // compute correlation
+  [ma,na] = size(a);
+  [mb,nb] = size(b);
+  c = conv2 (a, conj (b (mb:-1:1, nb:-1:1)));
+  // bias routines by Dave Cogdell (cogdelld@asme.org)
+  // optimized by Paul Kienzle (pkienzle@users.sf.net)
+  // coeff routine by Carnë Draug (carandraug+dev@gmail.com)
+  switch  (scale)
+    case {"none"}
+      // do nothing, it's all done
+    case {"biased"}
+      c = c / ( min ([ma, mb]) * min ([na, nb]) );
+    case {"unbiased"}
+      lo   = min ([na,nb]);
+      hi   = max ([na, nb]);
+      row  = [ 1:(lo-1), lo*ones(1,hi-lo+1), (lo-1):-1:1 ];
+      lo   = min ([ma,mb]);
+      hi   = max ([ma, mb]);
+      col  = [ 1:(lo-1), lo*ones(1,hi-lo+1), (lo-1):-1:1 ]';
+      bias = col*row;
+      c    = c./bias;
+    case {"coeff"}
+      a = double (a);
+      b = double (b);
+      a = conv2 (a.^2, ones (size (b,1) , size( b ,2)));
+      b = sum(b(:).* conj(b(:)));
+      c(:,:) = c(:,:) ./ sqrt (a(:,:) * b);
+    else
+      error ("xcorr2: invalid type of scale %s", scale);
+  end
 endfunction
-