Merge and solve ifft2 conflict

60 files changed, 1594 insertions, 38 deletions

diff --git a/macros/ar_psd.sci b/macros/ar_psd.sci
new file mode 100644
index 0000000..540dfec
--- /dev/null
+++ b/macros/ar_psd.sci
@@ -0,0 +1,40 @@
+function [P, F]= ar_psd(A, varargin)
+//Calculate the power spectrum of the autoregressive model
+//Calling Sequence
+// [PSD,F_OUT]=ar_psd (A, V)
+// [PSD,F_OUT]=ar_psd (A, V, FREQ)
+// [PSD,F_OUT]=ar_psd (A, V, FREQ, FS)
+// [PSD,F_OUT]=ar_psd (..., RANGE)
+// [PSD,F_OUT]=ar_psd (..., METHOD)
+// [PSD,F_OUT]=ar_psd (..., PLOTTYPE)
+//Parameters 
+//A:List of M=(order+1) autoregressive model coefficients. The first element of "ar_coeffs" is the zero-lag coefficient, which always has a value of 1.
+//V:Square of the moving-average coefficient of the AR model.
+//FREQ:Frequencies at which power spectral density is calculated, or a scalar indicating the number of uniformly distributed frequency values at which spectral density is calculated.  (default = 256)
+//FS:Sampling frequency (Hertz) (default=1)
+//Range: 'half', 'onesided' : frequency range of the spectrum is from zero up to but not including sample_f/2. Power from negative frequencies is added to the positive side of the spectrum.'whole', 'twosided' : frequency range of the spectrum is-sample_f/2 to sample_f/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 "freq" is vector, 'shift' is ignored.  If model coefficients "ar_coeffs" are real, the default range is 'half', otherwise default range is 'whole'.
+// Method:'fft': use FFT to calculate power spectrum.  'poly': calculate power spectrum as a polynomial of 1/z N.B. this argument is ignored if the "freq" argument is a vector.  The default is 'poly' unless the "freq" argument is an integer power of 2.
+// 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.
+//PSD: estimate of power-spectral density
+//F_OUT: frequency values
+//Description
+//If the FREQ argument is a vector (of frequencies) the spectrum is calculated using the polynomial method and the METHOD argument is ignored.  For scalar FREQ, an integer power of 2, or METHOD = "FFT", causes the spectrum to be calculated by FFT. Otherwise, the spectrum is calculated as a polynomial.  It may be computationally more efficient to use the FFT method if length of the model is not much smaller than the number of frequency values.  The spectrum is scaled so that spectral energy (area under spectrum) is the same as the time-domain energy (mean square of the signal).	
+//Examples
+//[a,b]= ar_psd([1,2,3],2)
+	
+	funcprot(0);
+	rhs= argn(2);
+	if(rhs <2 | rhs>5)
+		error("Wrong number of input arguments");
+	end
+	select(rhs)
+	case 2 then
+		[P,F]= callOctave("ar_psd", A, varargin(1));
+	case 3 then
+		[P,F]= callOctave("ar_psd", A, varargin(1), varargin(2));
+	case 4 then
+		[P,F]= callOctave("ar_psd", A, varargin(1), varargin(2), varargin(3));
+	case 5 then
+		[P,F]= callOctave("ar_psd", A, varargin(1), varargin(2), varargin(3), varargin(4));
+	end
+endfunction
diff --git a/macros/arch_fit.sci b/macros/arch_fit.sci
new file mode 100644
index 0000000..3e431a6
--- /dev/null
+++ b/macros/arch_fit.sci
@@ -0,0 +1,35 @@
+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
diff --git a/macros/arch_test.sci b/macros/arch_test.sci
new file mode 100644
index 0000000..3c53fc5
--- /dev/null
+++ b/macros/arch_test.sci
@@ -0,0 +1,46 @@
+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
diff --git a/macros/autoreg_matrix.sci b/macros/autoreg_matrix.sci
new file mode 100644
index 0000000..8778a84
--- /dev/null
+++ b/macros/autoreg_matrix.sci
@@ -0,0 +1,30 @@
+function y = autoreg_matrix(Y, varargin)
+// Given a time series (vector) Y, return a matrix with ones in the first column and the first K lagged values of Y in the other columns.
+//Calling Sequence
+//autoreg_matrix(Y, K)
+//Parameters
+//Y: Vector
+//K: Scalar or Vector
+//Description
+// Given a time series (vector) Y, return a matrix with ones in the first column and the first K lagged values of Y in the other columns.
+//
+//In other words, for T > K, '[1, Y(T-1), ..., Y(T-K)]' is the t-th row of the result.
+//
+//The resulting matrix may be used as a regressor matrix in autoregressions.
+//Examples
+//autoreg_matrix([1,2,3],2)
+//ans = 
+//      1.    0.    0.  
+//      1.    1.    0.  
+//      1.    2.    1.   
+funcprot(0);
+rhs = argn(2)
+if(rhs<2 | rhs>2)
+error("Wrong number of input arguments.");
+end
+
+	select(rhs)
+	case 2 then
+	y = callOctave("autoreg_matrix", Y, varargin(1));
+	end
+endfunction
diff --git a/macros/bilinear.sci b/macros/bilinear.sci
new file mode 100644
index 0000000..387b8d0
--- /dev/null
+++ b/macros/bilinear.sci
@@ -0,0 +1,40 @@
+function [Zb, Za, Zg]= bilinear(Sb,varargin)
+// Transform a s-plane filter specification into a z-plane specification
+//Calling Sequence
+// [ZB, ZA] = bilinear (SB, SA, T)
+// [ZB, ZA] = bilinear (SZ, SP, SG, T)
+// [ZZ, ZP, ZG] = bilinear (...)
+//Description
+//Transform a s-plane filter specification into a z-plane specification.  Filters can be specified in either zero-pole-gain or transfer function form.  The input form does not have to match the output form.  1/T is the sampling frequency represented in the z plane.
+//
+//Note: this differs from the bilinear function in the signal processing toolbox, which uses 1/T rather than T.
+//
+//Theory: Given a piecewise flat filter design, you can transform it from the s-plane to the z-plane while maintaining the band edges by means of the bilinear transform. This maps the left hand side of the s-plane into the interior of the unit circle. The mapping is highly non-linear, so you must design your filter with band edges in the s-plane positioned at 2/T tan(w*T/2) so that they will be positioned at w after the bilinear transform is complete.
+//Examples
+//[ZB,ZA]=bilinear([1],[2,3],3)
+	funcprot(0);
+	lhs= argn(1);
+	rhs= argn(2);
+	if(rhs < 3 | rhs > 4)
+		error("Wrong number of input arguments");
+	end
+	if(lhs < 2 | lhs > 3)
+		error("Wrong number of output arguments");
+	end
+	select(rhs)
+	case 3 then
+		select(lhs)
+		case 2 then
+			[Zb, Za]= callOctave("bilinear", Sb, varargin(1), varargin(2));
+		case 3 then
+			[Zb, Za, Zg]= callOctave("bilinear", Sb, varargin(1), varargin(2));
+		end
+	case 4 then
+		select(lhs)
+		case 2 then
+			[Zb, Za]= callOctave("bilinear", Sb, varargin(1), varargin(2), varargin(3));
+		case 3 then
+			[Zb, Za, Zg]= callOctave("bilinear", Sb, varargin(1), varargin(2), varargin(3));
+		end
+	end
+endfunction		
\ No newline at end of file
diff --git a/macros/cceps.sci b/macros/cceps.sci
index 99e4a7d..3756a5d 100644
--- a/macros/cceps.sci
+++ b/macros/cceps.sci
@@ -1,7 +1,21 @@
 function y = cceps (x,correct)
+//Return the complex cepstrum of the vector x
+//Calling Sequence
+//cceps (x)
+//cceps(x, correct)
+//Parameters
+//x: vector.
+//correct: if 1, a correction method  is applied.
+//Description
+//This function return the complex cepstrum of the vector x. If the optional argument correct has the value 1, a correction method is applied. The default is not to do this.
+//Examples
+//cceps([1,2,3],1)
+//ans = 
+//     1.92565
+//     0.96346
+//    -1.09735
 
 funcprot(0);
-//
 rhs = argn(2)
 if(rhs<1 | rhs>2)
 error("Wrong number of input arguments.")
diff --git a/macros/cohere.sci b/macros/cohere.sci
new file mode 100644
index 0000000..ba013ab
--- /dev/null
+++ b/macros/cohere.sci
@@ -0,0 +1,28 @@
+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
diff --git a/macros/cpsd.sci b/macros/cpsd.sci
new file mode 100644
index 0000000..5e56d7a
--- /dev/null
+++ b/macros/cpsd.sci
@@ -0,0 +1,50 @@
+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
+
+endfunction
diff --git a/macros/czt.sci b/macros/czt.sci
index b027489..84a0253 100644
--- a/macros/czt.sci
+++ b/macros/czt.sci
@@ -1,4 +1,23 @@
 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);
diff --git a/macros/detrend1.sci b/macros/detrend1.sci
new file mode 100644
index 0000000..52ac26f
--- /dev/null
+++ b/macros/detrend1.sci
@@ -0,0 +1,23 @@
+function y = detrend1(x, varargin)
+//This function removes the best fit of a polynomial of order P from the data X
+//Calling Sequence
+//detrend1(X,P)
+//Parameters 
+//X: Input vecor or matrix.
+//P: The order of polnomial
+//Description
+//If X is a vector, 'detrend1(X, P)' removes the best fit of apolynomial of order P from the data X.If X is a matrix, 'detrend1(X, P)' does the same for each column in X.
+//
+//The second argument P is optional.  If it is not specified, a value of 1 is assumed.  This corresponds to removing a linear trend.
+//The order of the polynomial can also be given as a string, in which case P must be either "constant" (corresponds to 'P=0') or "linear"(corresponds to 'P=1')
+	rhs= argn(2);
+	if(rhs<1 | rhs> 2)
+		error("Wrong number of input arguments");
+	end
+	select(rhs)
+	case 1 then
+		y= callOctave("detrend", x);
+	case 2 then
+		y= callOctave("detrend", x , varargin(1));
+	end
+endfunction
diff --git a/macros/diffpara.sci b/macros/diffpara.sci
new file mode 100644
index 0000000..c36ae7c
--- /dev/null
+++ b/macros/diffpara.sci
@@ -0,0 +1,50 @@
+function [D,DD] = diffpara(X,varargin)
+//Return the estimator D for the differencing parameter of an integrated time series
+//Calling Sequence
+// [D, DD] = diffpara (X)
+// [D, DD] = diffpara (X, A)
+// [D, DD] = diffpara (X, A, B)
+//Parameters
+//X: Input scalar or vector. 
+//DD:The estimators for all frequencies in the intervals described above.
+//D:The mean of DD 
+//Description
+//Return the estimator D for the differencing parameter of an integrated time series.
+//
+//The frequencies from [2*pi*a/t, 2*pi*b/T] are used for the estimation. If B is omitted, the interval [2*pi/T, 2*pi*a/T] is used.  If both B and A are omitted then a = 0.5 * sqrt (T) and b = 1.5 * sqrt (T) is used, where T is the sample size.  If X is a matrix, the differencing parameter of each column is estimated.
+//
+//The estimators for all frequencies in the intervals described above is returned in DD.
+//
+//The value of D is simply the mean of DD.
+	lhs= argn(1);
+	rhs= argn(2);
+	if(rhs <1 | rhs> 3)
+		error("Wrong number of input parameters");
+	end
+	if(lhs<1 | lhs>2)
+		error("Wrong number of output parameters");
+	end
+	select(rhs)
+	case 1 then
+		select(lhs)
+		case 1 then
+			D= diffpara(X);
+		case 2 then 
+			[D, DD]= diffpara(X);
+		end
+	case 2 then
+		select(lhs)
+		case 1 then
+			D= diffpara(X, varargin(1));
+		case 2 then 
+			[D, DD]= diffpara(X, varargin(1));
+		end
+	case 3 then
+		select(lhs)
+		case 1 then
+			D= diffpara(X, varargin(1), varargin(2));
+		case 2 then 
+			[D, DD]= diffpara(X, varargin(1), varargin(2));
+		end
+	end
+endfunction
diff --git a/macros/dst1.sci b/macros/dst1.sci
index 3b3b4b6..7f2165f 100644
--- a/macros/dst1.sci
+++ b/macros/dst1.sci
@@ -1,5 +1,14 @@
 function y = dst1(x, varargin)
-
+//Computes the type I discrete sine transform of x
+//Calling Sequence
+//y= dst1(x)
+//y= dst1(x,n)
+//Parameters 
+//x: real or complex valued vector
+//n= Length to which x is trimmed before transform 
+//Description
+//This is an Octave function.
+//Computes the 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.
 funcprot(0);
 
 lhs= argn(1);
diff --git a/macros/durbinlevinson.sci b/macros/durbinlevinson.sci
new file mode 100644
index 0000000..74dbb46
--- /dev/null
+++ b/macros/durbinlevinson.sci
@@ -0,0 +1,27 @@
+function y= durbinlevinson(C, varargin)
+// Perform one step of the Durbin-Levinson algorithm..
+//Calling Sequence
+// durbinlevinson (C);
+// durbinlevinson (C, OLDPHI);
+// durbinlevinson (C, OLDPHI, OLDV);
+//Parameters 
+//C: The vector C specifies the autocovariances '[gamma_0, ..., gamma_t]' from lag 0 to T.
+//OLDPHI: It specifies the coefficients based on C(T-1).
+//OLDV: It specifies the corresponding error.
+//Description
+//This is an Octave function.
+//Perform one step of the Durbin-Levinson.
+//If OLDPHI and OLDV are omitted, all steps from 1 to T of the algorithm are performed.
+	rhs=argn(2);
+	if(rhs<1 | rhs>3)
+		error("Wrong number of input arguments");
+	end
+	select(rhs)
+	case 1 then
+		y=callOctave("durbinlevinson",C);
+	case 2 then
+		y=callOctave("durbinlevinson",C, varargin(1));
+	case 3 then
+		y=callOctave("durbinlevinson",C, varargin(1), varargin(2));
+	end
+endfunction
\ No newline at end of file
diff --git a/macros/dwt.sci b/macros/dwt.sci
new file mode 100644
index 0000000..b564c7c
--- /dev/null
+++ b/macros/dwt.sci
@@ -0,0 +1,34 @@
+function [U, V] = dwt(X, varargin)
+//Discrete wavelet transform (1D)
+//Calling Sequence
+//[U, V] = dwt(X, WNAME)
+//[U, V] = dwt(X, HP, GP)
+//[U, V] = dwt(X, HP, GP,...)
+//Parameters
+//Inputs:
+//X: Signal Vector.
+//WNAME: Wavelet name.
+//HP: Coefficients of low-pass decomposition FIR filter.
+//GP: Coefficients of high-pass decomposition FIR filter.
+//Outputs:
+//U: Signal vector of average, approximation.
+//V: Signal vector of difference, detail.
+//Description
+//This function calculates the discrete wavelet transform (1D).
+//Examples
+//
+funcprot(0);
+rhs = argn(2)
+if(rhs<2 | rhs>4)
+error("Wrong number of input arguments.");
+end
+
+	select(rhs)
+	case 2 then
+	[U, V] = callOctave("dwt", X, varargin(1));
+    case 3 then
+    [U, V] = callOctave("dwt", X, varargin(1), varargin(2));
+    case 4 then
+    [U, V] = callOctave("dwt", X, varargin(1), varargin(2), varargin(3));
+	end
+endfunction
diff --git a/macros/fftconv.sci b/macros/fftconv.sci
new file mode 100644
index 0000000..d39441d
--- /dev/null
+++ b/macros/fftconv.sci
@@ -0,0 +1,26 @@
+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
diff --git a/macros/fftn.sci b/macros/fftn.sci
new file mode 100644
index 0000000..b8b4170
--- /dev/null
+++ b/macros/fftn.sci
@@ -0,0 +1,26 @@
+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
+endfunction
diff --git a/macros/fftshift1.sci b/macros/fftshift1.sci
new file mode 100644
index 0000000..9026025
--- /dev/null
+++ b/macros/fftshift1.sci
@@ -0,0 +1,33 @@
+function y= fftshift1(X,DIM)
+//Perform a shift of the vector X, for use with the 'fft1' and 'ifft1' functions, in order the move the frequency 0 to the center of the vector or matrix.
+//Calling Sequence
+// fftshift1 (X)
+// fftshift1 (X, DIM)
+//Parameters 
+//X:It is a vector of N elements corresponding to time samples
+//DIM: The optional DIM argument can be used to limit the dimension along which the permutation occurs
+//Description
+//This is an Octave function.
+//Perform a shift of the vector X, for use with the 'fft1' and 'ifft1' functions, in order the move the frequency 0 to the center of the vector or matrix.
+//
+//If X is a vector of N elements corresponding to N time samples spaced by dt, then 'fftshift1 (fft1 (X))' corresponds to frequencies
+//
+//f = [ -(ceil((N-1)/2):-1:1)*df 0 (1:floor((N-1)/2))*df ]
+//
+//where df = 1 / dt.
+//
+//If X is a matrix, the same holds for rows and columns.  If X is an array, then the same holds along each dimension.
+//
+//The optional DIM argument can be used to limit the dimension along
+     which the permutation occurs.
+	rhs= argn(2);
+	if(rhs <1 | rhs >2)
+		error('Wrong number of Input arguments');
+	end
+	select(rhs)
+	case 1 then
+		y=callOctave("fftshift",X);
+	case 2 then
+		y=callOctave("fftshift",X,DIM);
+	end
+endfunction
diff --git a/macros/fht.sci b/macros/fht.sci
index c6fd1bd..f7e5c8e 100644
--- a/macros/fht.sci
+++ b/macros/fht.sci
@@ -1,4 +1,17 @@
 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)
diff --git a/macros/filter1.sci b/macros/filter1.sci
new file mode 100644
index 0000000..3c32928
--- /dev/null
+++ b/macros/filter1.sci
@@ -0,0 +1,67 @@
+function [Y, SF] = filter1 (B, A, X, SI, DIM)
+//Apply a 1-D digital filter to the data X.
+//Calling Sequence
+//Y = filter1(B, A, X)
+//[Y, SF] = filter1(B, A, X, SI)
+//[Y, SF] = filter1(B, A, X, [], DIM)
+//[Y, SF] = filter1(B, A, X, SI, DIM)
+//Parameters
+//B: Matrix or Integer
+//A: Matrix or Integer
+//X: Matrix or Integer 
+//Description
+//'filter' returns the solution to the following linear, time-invariant difference equation:
+//
+//          N                   M
+//
+//         SUM a(k+1) y(n-k) = SUM b(k+1) x(n-k)    for 1<=n<=length(x)
+//
+//         k=0                 k=0
+//
+//where N=length(a)-1 and M=length(b)-1.  The result is calculated over the first non-singleton dimension of X or over DIM if supplied.
+//
+//An equivalent form of the equation is:
+//
+//                    N                   M
+//
+//          y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k)  for 1<=n<=length(x)
+//
+//                   k=1                 k=0
+//
+//    where c = a/a(1) and d = b/a(1).
+//Examples
+//filter([1,2,3], [3,4,5], [5,6,7])
+//ans = 
+//    1.6666667    3.1111111    4.4074074  
+funcprot(0);
+lhs = argn(1)
+rhs = argn(2)
+if (rhs < 3 | rhs > 5)
+error("Wrong number of input arguments.")
+end
+
+select(rhs)
+	
+	case 3 then
+		if(lhs==1)
+		Y=callOctave("filter",B,A,X)
+		elseif(lhs==2)
+		[Y, SF] = callOctave("filter",B,A,X)
+		else
+		error("Wrong number of output arguments.")
+		end
+	case 4 then
+		if(lhs==2)
+		[Y, SF] = callOctave("filter",B,A,X,SI)
+		else
+		error("Wrong number of output arguments.")
+	    end
+	case 5 then
+		if(lhs==2)
+		[Y, SF] = callOctave("filter",B,A,X,SI,DIM)
+		else
+		error("Wrong number of output arguments.")
+	    end
+	
+	end
+endfunction
diff --git a/macros/filter2.sci b/macros/filter2.sci
new file mode 100644
index 0000000..b8b48e4
--- /dev/null
+++ b/macros/filter2.sci
@@ -0,0 +1,33 @@
+function Y = filter2 (B, X, SHAPE)
+//Apply the 2-D FIR filter B to X.
+//Calling Sequence
+//Y = filter2(B, X)
+//Y = filter2(B, X, SHAPE)
+//Parameters
+//B, X: Matrix  
+// SHAPE: 
+//        'full':  pad X with zeros on all sides before filtering.
+//        'same': unpadded X (default)
+//        'valid': trim X after filtering so edge effects are no included. 
+//Description
+//This function applies the 2-D FIR filter B to X. If the argument SHAPE is specified, return an array of the desired shape.
+//Examples
+//filter2([1,3], [4,5])
+//ans = 
+//     19.    5.  
+funcprot(0);
+lhs = argn(1)
+rhs = argn(2)
+if (rhs < 2 | rhs > 3)
+error("Wrong number of input arguments.")
+end
+
+select(rhs)
+	
+	case 2 then
+        Y=callOctave("filter2",B,X)
+	case 3 then
+		Y = callOctave("filter2",B,X,SHAPE)
+    end
+
+endfunction
diff --git a/macros/findpeaks.sci b/macros/findpeaks.sci
new file mode 100644
index 0000000..3e7ce47
--- /dev/null
+++ b/macros/findpeaks.sci
@@ -0,0 +1,58 @@
+function [PKS, LOC, EXTRA] = findpeaks(DATA, varargin)
+//This function find peaks on DATA.
+//Calling Sequence
+//[PKS, LOC, EXTRA] = findpeaks(DATA)
+//[PKS, LOC, EXTRA] = findpeaks(..., PROPERTY, VALUE)
+//[PKS, LOC, EXTRA] = findpeaks(..., "DoubleSided")
+//Description
+//Peaks of a positive array of data are defined as local maxima. For double-sided data, they are maxima of the positive part and minima of the negative part. DATA is expected to be a single column vector.
+//
+//The function returns the value of DATA at the peaks in PKS. The index indicating their position is returned in LOC.
+//
+//The third output argument is a structure with additional information:
+//
+//"parabol":
+//          A structure containing the parabola fitted to each returned peak. The structure has two fields, "x" and "pp". The field "pp" contains the coefficients of the 2nd degree polynomial and "x" the extrema of the intercal here it was fitted.
+//
+//"height":
+//          The estimated height of the returned peaks (in units of DATA).
+//
+//"baseline":
+//          The height at which the roots of the returned peaks were calculated (in units of DATA).
+//
+//"roots":
+//        The abscissa values (in index units) at which the parabola fitted to each of the returned peaks crosses the "baseline" value. The width of the peak is calculated by 'diff(roots)'.
+//
+//This function accepts property-value pair given in the list below:
+//
+//"MinPeakHeight":
+//                Minimum peak height (positive scalar). Only peaks that exceed this value will be returned. For data taking positive and negative values use the option "DoubleSided". Default value '2*std (abs (detrend (data,0)))'.
+//
+//"MinPeakDistance":
+//                  Minimum separation between (positive integer). Peaks separated by less than this distance are considered a single peak.  This distance is also used to fit a second order polynomial to the peaks to estimate their width, therefore it acts as a smoothing parameter.  Default value 4.
+//
+//"MinPeakWidth":
+//               Minimum width of peaks (positive integer). The width of the peaks is estimated using a parabola fitted to the neighborhood of each peak.  The neighborhood size is equal to the value of "MinPeakDistance". The width is evaluated at the half height of the peak with baseline at "MinPeakHeight". Default value 2.
+//
+//"DoubleSided":
+//               Tells the function that data takes positive and negative values. The base-line for the peaks is taken as the mean value of the function.  This is equivalent as passing the absolute value of the data after removing the mean.
+funcprot(0);
+rhs=argn(2);
+lhs=argn(1)
+if(rhs<1 | rhs>2) then
+    error("Wrong number of input arguments.");
+end
+if(lhs<3 | lhs>3) then
+    error("Wrong number of output arguments.");
+end
+
+select(rhs)
+case 1 then
+    [PKS, LOC, EXTRA] = callOctave("findpeaks", DATA);
+case 2 then
+    [PKS, LOC, EXTRA] = callOctave("findpeaks", DATA, varargin(1));
+case 3 then
+    [PKS, LOC, EXTRA] = callOctave("findpeaks", DATA, varargin(1), varargin(2));
+end
+
+endfunction
diff --git a/macros/fir1.sci b/macros/fir1.sci
new file mode 100644
index 0000000..70c95b7
--- /dev/null
+++ b/macros/fir1.sci
@@ -0,0 +1,39 @@
+function B = fir1(N, W, varargin)
+//Produce an order N FIR filter with the given frequency cutoff, returning the N+1 filter coefficients in B.
+//Calling Sequence
+//B = fir1(N, W)
+//B = fir1(N, W, TYPE)
+//B = fir1(N, W, TYPE, WINDOW)
+//B = fir1(N, W, TYPE, WINDOW, NOSCALE)
+//Parameters
+//N: Integer
+//W: Integer or Vector
+//Description
+// Produce an order N FIR filter with the given frequency cutoff W, returning the N+1 filter coefficients in B. If W is a scalar, it specifies the frequency cutoff for a lowpass or highpass filter. If W is a two-element vector, the two values specify the edges of a bandpass or bandstop filter. If W is an N-element vector, each value specifies a band edge of a multiband pass/stop filter.
+//
+//The filter TYPE can be specified with one of the following strings: "low", "high", "stop", "pass", "bandpass", "DC-0", or "DC-1". The default is "low" is W is a scalar, "pass" if W is a pair, or "DC-0" if W is a vector with more than 2 elements.
+//
+//An optional shaping WINDOW can be given as a vector with length N+1. If not specified, a Hamming window of length N+1 is used.
+//
+//With the option "noscale", the filter coefficients are not normalized. The default is to normalize the filter such that the magnitude response of the center of the first passband is 1.
+//Examples
+// fir1 (5, 0.4)
+//ans = 
+//        9.2762e-05   9.5482e-02   4.0443e-01   4.0443e-01   9.5482e-02   9.2762e-05 
+funcprot(0);
+rhs = argn(2);
+if(rhs<2 | rhs>5)
+error("Wrong number of input arguments.");
+end
+
+	select(rhs)
+	case 2 then
+	B = callOctave("fir1", N, W);
+        case 3 then
+        B = callOctave("fir1", N, W, varargin(1));
+        case 4 then
+        B = callOctave("fir1", N, W, varargin(1), varargin(2));
+        case 5 then
+        B = callOctave("fir1", N, W, varargin(1), varargin(2), varargin(3));
+	end
+endfunction
diff --git a/macros/fir2.sci b/macros/fir2.sci
new file mode 100644
index 0000000..990fb3e
--- /dev/null
+++ b/macros/fir2.sci
@@ -0,0 +1,41 @@
+function B = fir2(N, F, M, varargin)
+//Produce an order N FIR filter with arbitrary frequency response M over frequency bands F, returning the N+1 filter coefficients in B.
+//Calling Sequence
+//B = fir2(N, F, M)
+//B = fir2(N, F, M, GRID_N)
+//B = fir1(N, F, M, GRID_N, RAMP_N)
+//B = fir1(N, F, M, GRID_N, RAMP_N, WINDOW)
+//Parameters
+//N: Integer
+//F, M: Vector
+//Description
+//Produce an order N FIR filter with arbitrary frequency response M over frequency bands F, returning the N+1 filter coefficients in B. The vector F specifies the frequency band edges of the filter response and M specifies the magnitude response at each frequency.
+//
+//The vector F must be nondecreasing over the range [0,1], and the first and last elements must be 0 and 1, respectively. A discontinuous jump in the frequency response can be specified by duplicating a band edge in F with different values in M.
+//
+//The resolution over which the frequency response is evaluated can be controlled with the GRID_N argument. The default is 512 or the next larger power of 2 greater than the filter length.
+//
+//The band transition width for discontinuities can be controlled with the RAMP_N argument. The default is GRID_N/25. Larger values will result in wider band transitions but better stopband rejection.
+//
+//An optional shaping WINDOW can be given as a vector with length N+1. If not specified, a Hamming window of length N+1 is used.
+//Examples
+// fir2 (10, [0, 0.5, 1], [1, 2, 3])
+//ans = 
+//     -0.00130   0.00000  -0.01792   0.00000  -0.36968   2.00000  -0.36968   0.00000  -0.01792   0.00000  -0.00130 
+funcprot(0);
+rhs = argn(2);
+if(rhs<3 | rhs>6)
+error("Wrong number of input arguments.");
+end
+
+	select(rhs)
+	case 3 then
+	B = callOctave("fir2", N, F, M);
+        case 4 then
+        B = callOctave("fir2", N, F, M, varargin(1));
+        case 5 then
+        B = callOctave("fir2", N, F, M, varargin(1), varargin(2));
+        case 6 then
+        B = callOctave("fir2", N, F, M, varargin(1), varargin(2), varargin(3));
+	end
+endfunction
diff --git a/macros/fractdiff.sci b/macros/fractdiff.sci
new file mode 100644
index 0000000..979a079
--- /dev/null
+++ b/macros/fractdiff.sci
@@ -0,0 +1,17 @@
+function y= fractdiff(x,d)
+//Compute the fractional differences (1-L)^d x where L denotes the lag-operator and d is greater than -1.
+//Calling Sequence
+// fractdiff (X, D)
+//Description
+//This is an Octave function.
+//Compute the fractional differences (1-L)^d x where L denotes the lag-operator and d is greater than -1.
+	funcprot(0);
+	rhs= argn(2);
+	if(rhs < 2 | rhs >2)
+		error("Wrong number of input arguments");
+	end
+	select(rhs)
+	case 2 then
+		y= callOctave("fractdiff",x,d);
+	end
+endfunction
\ No newline at end of file
diff --git a/macros/freqs.sci b/macros/freqs.sci
index 88773fe..05dd294 100644
--- a/macros/freqs.sci
+++ b/macros/freqs.sci
@@ -24,6 +24,6 @@ end
 
 select (rhs)
 	case 3 then
-		y = callOctave("freqs",b, a, w)
+		h = callOctave("freqs",b, a, w)
 	end
 endfunction
diff --git a/macros/freqz.sci b/macros/freqz.sci
new file mode 100644
index 0000000..dc71b75
--- /dev/null
+++ b/macros/freqz.sci
@@ -0,0 +1,71 @@
+function [H, W] = freqz(B, varargin) 
+//This function returns the complex frequency response H of the rational IIR filter whose numerator and denominator coefficients are B and A, respectively.
+//Calling Sequence
+//[H, W] = freqz(B, A, N, "whole")
+//[H, W] = freqz(B)
+//[H, W] = freqz(B, A)
+//[H, W] = freqz(B, A, N)
+//H = freqz(B, A, W)
+//[H, W] = freqz(..., FS)
+//freqz(...)
+//Parameters
+//B, A, N: Integer or Vector
+//Description
+// Return the complex frequency response H of the rational IIR filter whose numerator and denominator coefficients are B and A, respectively.
+//
+//The response is evaluated at N angular frequencies between 0 and 2*pi.
+//
+//The output value W is a vector of the frequencies.
+// 
+//If A is omitted, the denominator is assumed to be 1 (this corresponds to a simple FIR filter).
+//
+//If N is omitted, a value of 512 is assumed. For fastest computation, N should factor into a small number of small primes.
+//
+//If the fourth argument, "whole", is omitted the response is evaluated at frequencies between 0 and pi.
+//
+//     'freqz (B, A, W)'
+//
+//Evaluate the response at the specific frequencies in the vector W. The values for W are measured in radians.
+//
+//     '[...] = freqz (..., FS)'
+//
+//Return frequencies in Hz instead of radians assuming a sampling rate FS. If you are evaluating the response at specific frequencies W, those frequencies should be requested in Hz rather than radians.
+//
+//     'freqz (...)'
+//
+//Plot the magnitude and phase response of H rather than returning them.
+//Examples
+//H = freqz([1,2,3], [4,3], [1,2,5])
+//ans = 
+//       0.4164716 - 0.5976772i  - 0.4107690 - 0.2430335i    0.1761948 + 0.6273032i  
+funcprot(0);
+rhs=argn(2);
+lhs=argn(1);
+if(rhs<2 | rhs>4) then
+    error("Wrong number of input arguments.");
+end
+if (lhs<1 | lhs>2)
+    error("Wrong number of output arguments.");
+end
+if (lhs==1) then
+select(rhs)
+case 1 then
+    H = callOctave("freqz",B);
+case 2 then
+    H = callOctave("freqz",B, varargin(1));
+case 3 then 
+    H = callOctave("freqz",B, varargin(1), varargin(2));
+end
+elseif (lhs==2) then 
+    select(rhs)
+case 1 then
+    [H, W] = callOctave("freqz",B);
+case 2 then
+    [H, W] = callOctave("freqz",B, varargin(1));
+case 3 then 
+    [H, W] = callOctave("freqz",B, varargin(1), varargin(2));
+case 4 then 
+    [H, W] = callOctave("freqz", B, varargin(1), varargin(2), varargin(3));
+end    
+end
+endfunction
diff --git a/macros/fwht.sci b/macros/fwht.sci
index 280bbee..e0a15ff 100644
--- a/macros/fwht.sci
+++ b/macros/fwht.sci
@@ -1,4 +1,28 @@
 function y = fwht(x, varargin)
+//Compute the Walsh-Hadamard transform of x using the Fast Walsh-Hadamard Transform (FWHT) algorithm
+//Calling Sequence
+//fwht (x)
+//fwht (x, n)
+//fwht (x, n, order)
+//Parameters 
+//x: real or complex valued scalar or vector
+//n: x is truncated or extended to have length n 
+//order: Specification of order in which coefficients should be arranged 
+//Description
+//Compute the Walsh-Hadamard transform of x using the Fast Walsh-Hadamard Transform (FWHT) algorithm. If the input is a matrix, the FWHT is calculated along the columns of x.
+//
+//The number of elements of x must be a power of 2; if not, the input will be extended and filled with zeros. If a second argument is given, the input is truncated or extended to have length n.
+//
+//The third argument specifies the order in which the returned Walsh-Hadamard transform coefficients should be arranged. The order may be any of the following strings:
+//
+//"sequency"
+//The coefficients are returned in sequency order. This is the default if order is not given.
+//
+//"hadamard"
+//The coefficients are returned in Hadamard order.
+//
+//"dyadic"
+//The coefficients are returned in Gray code order.
 
 funcprot(0);
 rhs= argn(2);
diff --git a/macros/hamming.sci b/macros/hamming.sci
new file mode 100644
index 0000000..49a970f
--- /dev/null
+++ b/macros/hamming.sci
@@ -0,0 +1,24 @@
+function y = hamming(m, varargin)
+//Return the filter coefficients of a Hamming window of length M
+//Calling Sequence
+//hamming (M)
+//hamming (M, "periodic")
+//hamming (M, "symmetric")
+//Parameters 
+//M: real scalar, which will be the length of hamming window
+//Description
+//Return the filter coefficients of a Hamming window of length M.
+//If the optional argument "periodic" is given, the periodic form of the window is returned.  This is equivalent to the window of length M+1 with the last coefficient removed.  The optional argument "symmetric" is equivalent to not specifying a second argument.
+funcprot(0);
+rhs= argn(2);
+if(rhs <1 | rhs>2)
+error("Wrong number of Input parameters");
+end
+
+select(rhs)
+	case 1 then
+		y= callOctave("hamming", m);
+	case 2 then
+		y= callOctave("hamming", m , varargin(1));
+end
+endfunction
diff --git a/macros/hanning.sci b/macros/hanning.sci
new file mode 100644
index 0000000..60ca783
--- /dev/null
+++ b/macros/hanning.sci
@@ -0,0 +1,25 @@
+function y = hanning(m, varargin)
+//Return the filter coefficients of a Hanning window of length M
+//Calling Sequence
+//hanning (M)
+//hanning (M, "periodic")
+//hanning (M, "symmetric")
+//Parameters 
+//M: real scalar, which will be the length of hanning window
+//Description
+//Return the filter coefficients of a Hanning window of length M.
+//If the optional argument "periodic" is given, the periodic form of the window is returned.  This is equivalent to the window of length M+1 with the last coefficient removed.  The optional argument "symmetric" is equivalent to not specifying a second argument.	
+
+funcprot(0);
+rhs= argn(2);
+if(rhs <1 | rhs>2)
+error("Wrong number of Input parameters");
+end
+
+select(rhs)
+	case 1 then
+		y= callOctave("hanning", m);
+	case 2 then
+		y= callOctave("hanning", m , varargin(1));
+end
+endfunction
diff --git a/macros/hilbert1.sci b/macros/hilbert1.sci
index 864067f..fbcb136 100644
--- a/macros/hilbert1.sci
+++ b/macros/hilbert1.sci
@@ -1,4 +1,26 @@
 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);
diff --git a/macros/hurst.sci b/macros/hurst.sci
index 1a99f8c..27507fb 100644
--- a/macros/hurst.sci
+++ b/macros/hurst.sci
@@ -1,5 +1,13 @@
 function y = hurst(x)
-
+//  Estimate the Hurst parameter of sample X via the rescaled r statistic.
+//Calling Sequence
+//hurst(X)
+//variable=hurst(X)
+//Parameters 
+//X: X is a matrix, the parameter of sample X via the rescaled r statistic
+//Description
+//This is an Octave function.
+//This function estimates the Hurst parameter of sample X via the rescaled rstatistic.
 funcprot(0);
 rhs= argn(2);
 if(rhs<1 | rhs>1)
diff --git a/macros/idct1.sci b/macros/idct1.sci
index d398e56..5015187 100644
--- a/macros/idct1.sci
+++ b/macros/idct1.sci
@@ -1,14 +1,27 @@
 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)
+rhs=argn(2);
 if (rhs<1 | rhs>2) then
-    error("Wrong number of input arguments.")
+    error("Wrong number of input arguments.");
 end
 select(rhs)
 case 1 then
-    y=callOctave("idct",x)
+    y=callOctave("idct",x);
 case 2 then 
-    y=callOctave("idct",x,n)
+    y=callOctave("idct",x,n);
 end
 
 endfunction
diff --git a/macros/idct2.sci b/macros/idct2.sci
index 1cfacab..c48980a 100644
--- a/macros/idct2.sci
+++ b/macros/idct2.sci
@@ -1,8 +1,22 @@
 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)
+rhs=argn(2);
 if (rhs<1 | rhs>3) then
-    error("Wrong number of input arguments.")
+    error("Wrong number of input arguments.");
 end
 select(rhs)
 case 1 then
diff --git a/macros/idst1.sci b/macros/idst1.sci
index 5276a70..98fc874 100644
--- a/macros/idst1.sci
+++ b/macros/idst1.sci
@@ -1,4 +1,17 @@
 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
diff --git a/macros/ifft21.sci b/macros/ifft2.sci
index 25ed82f..8ecb1a5 100644
--- a/macros/ifft21.sci
+++ b/macros/ifft2.sci
@@ -1,4 +1,4 @@
-function res = ifft21 (A, m, n)
+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)
@@ -15,7 +15,7 @@ function res = ifft21 (A, m, n)
 //x = [1 2 3; 4 5 6; 7 8 9]
 //m = 4
 //n = 4
-//ifft21 (A, m, n)
+//ifft2 (A, m, n)
 //ans =
 //
 //   2.81250 + 0.00000i  -0.37500 + 0.93750i   0.93750 + 0.00000i  -0.37500 - 0.93750i
diff --git a/macros/ifftn.sci b/macros/ifftn.sci
new file mode 100644
index 0000000..3d26c04
--- /dev/null
+++ b/macros/ifftn.sci
@@ -0,0 +1,27 @@
+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
diff --git a/macros/ifftshift1.sci b/macros/ifftshift1.sci
new file mode 100644
index 0000000..7426130
--- /dev/null
+++ b/macros/ifftshift1.sci
@@ -0,0 +1,24 @@
+function y= ifftshift1(X,DIM)
+//Undo the action of the 'fftshift1' function.
+//Calling Sequence
+// ifftshift1 (X)
+// ifftshift1 (X, DIM)
+//Parameters 
+//X:It is a vector of N elements corresponding to time samples
+//DIM: The optional DIM argument can be used to limit the dimension along which the permutation occurs
+//Description
+//This is an Octave function.
+//Undo the action of the 'fftshift1' function.
+//
+//For even length X, 'fftshift1' is its own inverse, but odd lengths differ slightly.
+	rhs= argn(2);
+	if(rhs <1 | rhs >2)
+		error('Wrong number of Input arguments');
+	end
+	select(rhs)
+	case 1 then
+		y=callOctave("ifftshift",X);
+	case 2 then
+		y=callOctave("ifftshift",X,DIM);
+	end
+endfunction
\ No newline at end of file
diff --git a/macros/ifht.sci b/macros/ifht.sci
index aacbd7a..7f865cb 100644
--- a/macros/ifht.sci
+++ b/macros/ifht.sci
@@ -1,4 +1,22 @@
 function m = ifht(d, varargin)
+//Calculate the inverse Fast Hartley Transform of real input D
+//Calling Sequence
+//m= ifht (d)
+//m= ifht (d,n)
+//m= ifht (d,n,dim)
+//Parameters 
+//d: real or complex valued scalar or vector
+//n: Similar to the options of FFT function
+//dim: Similar to the options of FFT function 
+//Description
+//Calculate the inverse Fast Hartley Transform of real input d. If d is a matrix, the inverse Hartley transform is calculated along the columns by default. The options n and dim are similar to the options of FFT function.
+//
+//The forward and inverse Hartley transforms are the same (except for a scale factor of 1/N for the inverse hartley transform), but implemented using different functions.
+//
+//The definition of the forward hartley transform for vector d, m[K] = 1/N \sum_{i=0}^{N-1} d[i]*(cos[K*2*pi*i/N] + sin[K*2*pi*i/N]), for 0 <= K < N. m[K] = 1/N \sum_{i=0}^{N-1} d[i]*CAS[K*i], for 0 <= K < N.
+//Examples
+//ifht(1 : 4)
+//ifht(1:4, 2)
 funcprot(0);
 rhs= argn(2);
 if(rhs<1 | rhs>3)
diff --git a/macros/ifwht.sci b/macros/ifwht.sci
index 163c032..c1517ac 100644
--- a/macros/ifwht.sci
+++ b/macros/ifwht.sci
@@ -1,4 +1,27 @@
 function y= ifwht(x, varargin)
+// Compute the inverse Walsh-Hadamard transform of x using the Fast Walsh-Hadamard Transform (FWHT) algorithm
+// Calling Sequence
+// ifwht (x)
+// ifwht (x, n)
+// ifwht (x, n, order)
+//Parameters 
+//x: real or complex valued scalar or vector
+//n: Input is truncated or extended to have a length of n
+//order: Specifies the order in which the returned inverse Walsh-Hadamard transform
+//Description
+//Compute the inverse Walsh-Hadamard transform of x using the Fast Walsh-Hadamard Transform (FWHT) algorithm. If the input is a matrix, the inverse FWHT is calculated along the columns of x.
+//The number of elements of x must be a power of 2; if not, the input will be extended and filled with zeros. If a second argument is given, the input is truncated or extended to have length n.
+//The third argument specifies the order in which the returned inverse Walsh-Hadamard transform coefficients should be arranged. The order may be any of the following strings:
+//
+//"sequency"
+//The coefficients are returned in sequency order. This is the default if order is not given.
+//
+//"hadamard"
+//The coefficients are returned in Hadamard order.
+//
+//"dyadic"
+//The coefficients are returned in Gray code order.
+
 
 funcprot(0);
 rhs= argn(2);
diff --git a/macros/invfreq.sci b/macros/invfreq.sci
index dde23ef..a720115 100644
--- a/macros/invfreq.sci
+++ b/macros/invfreq.sci
@@ -1,6 +1,28 @@
 function [B,A] = invfreq(H,F,nB,nA,W,iter,tol, plane)
-
-
+// Calculates inverse frequency vectors
+//
+// Calling Sequence
+//[B,A] = invfreq(H,F,nB,nA)
+//[B,A] = invfreq(H,F,nB,nA,W)
+//[B,A] = invfreq(H,F,nB,nA,W,[],[],plane)
+//[B,A] = invfreq(H,F,nB,nA,W,iter,tol,plane)
+// 
+// Parameters
+// H: desired complex frequency response,It is assumed that A and B are real polynomials, hence H is one-sided.
+// F: vector of frequency samples in radians
+// nA: order of denominator polynomial A
+// nB: order of numerator polynomial B
+//
+// Description
+//Fit filter B(z)/A(z) or B(s)/A(s) to complex frequency response at frequency points F. A and B are real polynomial coefficients of order nA and nB respectively.  Optionally, the fit-errors can be weighted vs frequency according to the weights W. Also, the transform plane can be specified as either 's' for continuous time or 'z' for discrete time. 'z' is chosen by default.  Eventually, Steiglitz-McBride iterations will be specified by iter and tol.
+//
+// Examples
+//  [B,A] = butter(12,1/4);
+//  [H,w] = freqz(B,A,128);
+//  [Bh,Ah] = invfreq(H,F,4,4);
+//  Hh = freqz(Bh,Ah);
+//  disp(sprintf('||frequency response error|| = %f',norm(H-Hh)));
+// 
 funcprot(0);
 lhs= argn(1);
 rhs= argn(2);
diff --git a/macros/lib b/macros/lib
index 47dc348..54f79c5 100644
--- a/macros/lib
+++ b/macros/lib
diff --git a/macros/mscohere.sci b/macros/mscohere.sci
new file mode 100644
index 0000000..357da58
--- /dev/null
+++ b/macros/mscohere.sci
@@ -0,0 +1,82 @@
+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
+
diff --git a/macros/names b/macros/names
index b3d1310..800df1a 100644
--- a/macros/names
+++ b/macros/names
@@ -1,16 +1,21 @@
 ac2poly
 ac2rc
 arParEst
+ar_psd
 arburg
 arch_rnd
+arch_test
+arch_fit
 arcov
 arma_rnd
 armcov
 aryule
+autoreg_matrix
 barthannwin
 bartlett
 besselap
 besself
+bilinear
 bitrevorder
 blackman
 blackmanharris
@@ -37,9 +42,11 @@ chirp
 cl2bp
 clustersegment
 cmorwavf
+cohere
 convmtx
 corrmtx
 cplxreal
+cpsd
 cummax
 cummin
 czt
@@ -47,11 +54,15 @@ db
 db2pow
 dctmtx
 decimate
+detrend1
 dftmtx
+diffpara
 diric
 downsample
 dst1
+durbinlevinson
 dutycycle
+dwt
 ellip
 ellipap
 ellipord
@@ -60,18 +71,28 @@ eqtflength
 falltime
 fft1
 fft21
+fftconv
 fftfilt
 fftw1
+fftn
+fftshift1
 fht
+filter1
+filter2
 filternorm
 filtfilt
 filtic
 filtord
+findpeaks
+fir1
+fir2
 firpmord
 firtype
 flattopwin
 fracshift
+fractdiff
 freqs
+freqz
 fwhm
 fwhmjlt
 fwht
@@ -82,7 +103,9 @@ gausswin
 gmonopuls
 goertzel
 grpdelay
+hamming
 hann
+hanning
 helperHarmonicDistortionAmplifier
 hilbert1
 hurst
@@ -92,6 +115,8 @@ idct2
 idst1
 ifft1
 ifft21
+ifftn
+ifftshift1
 ifht
 ifwht
 iirlp2mb
@@ -131,6 +156,7 @@ midcross
 modulate
 morlet
 movingrms
+mscohere
 musicBase
 ncauer
 nnls
@@ -163,6 +189,7 @@ pulsesep
 pulsewidth
 pulstran
 pyulear
+pwelch
 qp_kaiser
 rc2ac
 rc2is
@@ -217,6 +244,7 @@ tf2sos
 tf2zp
 tf2zpk
 tfestimate
+tfe
 transpose
 trial_iirlp2mb
 triang
@@ -227,6 +255,7 @@ udecode
 uencode
 ultrwin
 unshiftdata
+unwrap2
 upfirdn
 upsample
 upsamplefill
diff --git a/macros/periodogram.sci b/macros/periodogram.sci
index 33e15af..ca8995e 100755
--- a/macros/periodogram.sci
+++ b/macros/periodogram.sci
@@ -1,20 +1,46 @@
 function [d,n]=periodogram(a,b,c,d,e)
+//Return the periodogram (Power Spectral Density) of X
+//Calling Sequence
+// [PXX, W] = periodogram (X)
+// [PXX, W] = periodogram (X, WIN)
+// [PXX, W] = periodogram (X, WIN, NFFT)
+// [PXX, W] = periodogram (X, WIN, NFFT, FS)
+// [PXX, W] = periodogram (..., "RANGE")
+//Parameters 
+// X:data vector.  If X is real-valued a one-sided spectrum is estimated. If X is complex-valued, or "RANGE" specifies "twosided", the full spectrum is estimated.
+//WIN: window weight data.  If window is empty or unspecified a default rectangular window is used.  Otherwise, the window is applied to the signal ('X .* WIN') before computing th periodogram.  The window data must be a vector of the same length as X.
+//NFFT:number of frequency bins.  The default is 256 or the next higher power of 2 greater than the length of X ('max (256,2.^nextpow2 (length (x)))').  If NFFT is greater than the length of the input then X will be zero-padded to the length of NFFT.
+//FS:sampling rate.  The default is 1.
+//RANGE:range of spectrum.  "onesided" computes spectrum from [0..nfft/2+1]."twosided" computes spectrum from [0..nfft-1].
+//Description
+//The optional second output W are the normalized angular frequencies.  For a one-sided calculation W is in the range [0, pi]. If NFFT is even and [0, pi) if NFFT is odd.  Similarly, for a two-sided calculation W is in the range [0, 2*pi] or [0, 2*pi)depending on NFFT.
+//
+//If a sampling frequency is specified, FS, then the output frequencies F will be in the range [0, FS/2] or [0, FS/2) for one-sided calculations. For two-sided calculations the range will be [0, FS).
+//
+//When called with no outputs the periodogram is immediately plotted in the current figure window.
     funcprot(0);
-    [nargout,nargin]=argn();
-    select nargin
+    lhs= argn(1);
+    rhs= argn(2);
+    if(rhs<1 | rhs>5)
+        error("Wrong number of input arguments");
+    end
+    if(lhs>2 | lhs< 2)
+        error("Wrong number of output arguments");
+    end
+    select(rhs)
     case 1 then
-        [d,n]=callOctave('periodogram',a);
+            [d,n]= callOctave('periodogram',a);
     case 2 then
-        [d,n]=callOctave('periodogram',a,b);
+            [d,n]= callOctave('periodogram',a,b);
+
     case 3 then
-        [d,n]=callOctave('periodogram',a,b,c);
+            [d,n]= callOctave('periodogram',a,b,c);
+
     case 4 then
-        [d,n]=callOctave('periodogram',a,b,c,d);  
+            [d,n]= callOctave('periodogram',a,b,c,d);
+  
     case 5 then
-        [d,n]=callOctave('periodogram',a,b,c,d,e);    
-    else
-        error("Incorrect no. of Input Arguments");
+            [d,n]= callOctave('periodogram',a,b,c,d,e);
     end
-
 endfunction
 
diff --git a/macros/pwelch.sci b/macros/pwelch.sci
new file mode 100644
index 0000000..3e4aa8d
--- /dev/null
+++ b/macros/pwelch.sci
@@ -0,0 +1,67 @@
+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));
+    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));
+    end
+end
+
+endfunction
diff --git a/macros/rceps.sci b/macros/rceps.sci
index 065e482..8b4d89a 100644
--- a/macros/rceps.sci
+++ b/macros/rceps.sci
@@ -1,5 +1,17 @@
 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
 funcprot(0)
 lhs= argn(1)
 rhs= argn(2)
diff --git a/macros/remez1.sci b/macros/remez1.sci
index aabb25e..d9cfd29 100644
--- a/macros/remez1.sci
+++ b/macros/remez1.sci
@@ -1,4 +1,19 @@
 function b = remez1(n,f,a, varargin)
+//Parks-McClellan optimal FIR filter design
+//Calling Sequence
+//b = remez1 (n, f, a)
+//b = remez1 (n, f, a, w)
+//b = remez1 (n, f, a, w, ftype)
+//b = remez1 (n, f, a, w, ftype, griddensity)
+//Parameters 
+//n: gives the number of taps in the returned filter
+//f:gives frequency at the band edges [b1 e1 b2 e2 b3 e3 …]
+//a:gives amplitude at the band edges [a(b1) a(e1) a(b2) a(e2) …]
+//w:gives weighting applied to each band
+//ftype:is "bandpass", "hilbert" or "differentiator"
+//griddensity:determines how accurately the filter will be constructed. The minimum value is 16, but higher numbers are slower to compute.
+//Description
+// Frequency is in the range (0, 1), with 1 being the Nyquist frequency.
 
 funcprot(0);
 rhs= argn(2);
diff --git a/macros/sigmoid_train.sci b/macros/sigmoid_train.sci
index bba751d..181584a 100644
--- a/macros/sigmoid_train.sci
+++ b/macros/sigmoid_train.sci
@@ -1,4 +1,16 @@
 function y =sigmoid_train(t, ranges, rc)
+// Evaluate a train of sigmoid functions at T.
+//Calling Sequence
+//y = sigmoid_train(t, ranges, rc)
+//Parameters
+//t: integer
+//ranges: matrix
+//Description
+//The number and duration of each sigmoid is determined from RANGES. Each row of RANGES represents a real interval, e.g.  if sigmoid 'i' starts at 't=0.1' and ends at 't=0.5', then 'RANGES(i,:) = [0.1 0.5]'.  The input RC is an array that defines the rising and falling time constants of each sigmoid.  Its size must equal the size of RANGES.
+//Examples
+//sigmoid_train(0.1,[1:3],4)
+//ans = 
+//     0.27375
 funcprot(0);
 rhs=argn(2);
 if (rhs<3 | rhs>3) then
diff --git a/macros/sinetone.sci b/macros/sinetone.sci
index 006a2e0..407fa45 100644
--- a/macros/sinetone.sci
+++ b/macros/sinetone.sci
@@ -1,5 +1,17 @@
 function y= sinetone(x, varargin)
-
+//Return a sinetone of the input
+//Calling Sequence
+//y= sinetone(FREQ)
+//y= sinetone(FREQ, RATE)
+//y= sinetone(FREQ, RATE, SEC)
+//y= sinetone(FREQ, RATE, SEC, AMPL)
+//Parameters 
+//FREQ: frequency of sinetone
+//RATE: Sampling rate
+//SEC: Length in seconds
+//AMPL: Amplitude
+//Description
+//Return a sinetone of frequency FREQ with a length of SEC seconds atsampling rate RATE and with amplitude AMPL.The arguments FREQ and AMPL may be vectors of common size.The defaults are RATE = 8000, SEC = 1, and AMPL = 64.
 funcprot(0);
 rhs= argn(2);
 if(rhs<1 | rhs>4)
diff --git a/macros/sinewave.sci b/macros/sinewave.sci
index d5bdfa9..5263f35 100644
--- a/macros/sinewave.sci
+++ b/macros/sinewave.sci
@@ -1,5 +1,16 @@
 function y= sinewave(x, varargin)
-
+//Return an M-element vector with I-th element given by 'sin(2* pi *(I+D-1)/N).'
+//Calling Sequence
+//y= sinewave(M)
+//y= sinewave(M,N)
+//y= sinewave(M,N,D)
+//Parameters 
+//M: Input vector
+//N: The default value for N is M
+//D: The default value for D is 0 
+//AMPL: Amplitude
+//Description
+//Return an M-element vector with I-th element given by 'sin(2* pi *(I+D-1)/N).'
 funcprot(0);
 rhs= argn(2);
 if(rhs<1 | rhs>3)
diff --git a/macros/spectral_adf.sci b/macros/spectral_adf.sci
index ea81f42..38b3990 100644
--- a/macros/spectral_adf.sci
+++ b/macros/spectral_adf.sci
@@ -1,5 +1,20 @@
 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)
diff --git a/macros/spectral_xdf.sci b/macros/spectral_xdf.sci
index fe93327..f0d457b 100644
--- a/macros/spectral_xdf.sci
+++ b/macros/spectral_xdf.sci
@@ -1,5 +1,18 @@
 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)
diff --git a/macros/spencer.sci b/macros/spencer.sci
index 537ca8a..63a1b83 100644
--- a/macros/spencer.sci
+++ b/macros/spencer.sci
@@ -1,5 +1,11 @@
 function y= spencer(x)
-
+//Return Spencer's 15 point moving average of each column of X.
+//Calling Sequence
+//spencer(X)
+//Parameters 
+//X: Real scalar or vector
+//Description
+//Return Spencer's 15 point moving average of each column of X.
 funcprot(0);
 
 rhs= argn(2);
diff --git a/macros/stft.sci b/macros/stft.sci
index 01d8f9b..7c30360 100644
--- a/macros/stft.sci
+++ b/macros/stft.sci
@@ -1,4 +1,39 @@
 function [y,c]= stft(x, varargin)
+//Compute the short-time Fourier transform of the vector X
+//Calling Sequence
+//Y = stft (X)
+//Y = stft (X, WIN_SIZE)
+//Y = stft (X, WIN_SIZE, INC)
+//Y = stft (X, WIN_SIZE, INC, NUM_COEF)
+//Y = stft (X, WIN_SIZE, INC, NUM_COEF, WIN_TYPE)
+//[Y,C] = stft (X)
+//[Y,C] = stft (X, WIN_SIZE)
+//[Y,C] = stft (X, WIN_SIZE, INC)
+//[Y,C] = stft (X, WIN_SIZE, INC, NUM_COEF)
+//[Y,C] = stft (X, WIN_SIZE, INC, NUM_COEF, WIN_TYPE)
+//Parameters 
+//X: Real scalar or vector
+//WIN_SIZE: Size of the window used
+//INC: Increment
+//WIN_TYPE: Type of window
+//Description
+//Compute the short-time Fourier transform of the vector X with NUM_COEF coefficients by applying a window of WIN_SIZE data points and an increment of INC points.
+//
+//Before computing the Fourier transform, one of the following windows is applied:
+//
+//"hanning" -> win_type = 1
+//
+//"hamming" -> win_type = 2
+//
+//"rectangle" -> win_type = 3
+//
+//The window names can be passed as strings or by the WIN_TYPE number.
+//
+//The following defaults are used for unspecified arguments:WIN_SIZE= 80, INC = 24, NUM_COEF = 64, and WIN_TYPE = 1.
+//
+//Y = stft (X, ...)' returns the absolute values of the Fourier coefficients according to the NUM_COEF positive frequencies.
+//
+//'[Y, C] = stft (x, ...)' returns the entire STFT-matrix Y and a 3-element vector C containing the window size, increment, and window type, which is needed by the 'synthesis' function.
 
 funcprot(0);
 lhs= argn(1);
diff --git a/macros/synthesis.sci b/macros/synthesis.sci
index e787c86..7224686 100644
--- a/macros/synthesis.sci
+++ b/macros/synthesis.sci
@@ -1,5 +1,14 @@
 function x= synthesis(Y,C)
-
+//Compute a signal from its short-time Fourier transform
+//Calling Sequence
+//X= synthesis(Y,C)
+//Parameters 
+//Y: Shirt-time fourier transform
+//C: 3-element vector C specifying window size, increment, window type.
+//Description
+//Compute a signal from its short-time Fourier transform Y and a 3-element vector C specifying window size, increment, and window type.
+//The values Y and C can be derived by
+//[Y, C] = stft (X , ...)
 funcprot(0);
 lhs= argn(1);
 rhs= argn(2);
diff --git a/macros/tfe.sci b/macros/tfe.sci
new file mode 100644
index 0000000..5acc391
--- /dev/null
+++ b/macros/tfe.sci
@@ -0,0 +1,27 @@
+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
diff --git a/macros/unwrap2.sci b/macros/unwrap2.sci
new file mode 100644
index 0000000..b8ea9de
--- /dev/null
+++ b/macros/unwrap2.sci
@@ -0,0 +1,35 @@
+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
+
+endfunction
diff --git a/macros/xcorr2.sci b/macros/xcorr2.sci
index e73f00f..688dc3a 100644
--- a/macros/xcorr2.sci
+++ b/macros/xcorr2.sci
@@ -4,12 +4,10 @@ function c = xcorr2 (a, b, biasflag)
 //c = xcorr2 (a)
 //c = xcorr2 (a, b)
 //c = xcorr2 (a, b, biasflag)
-
 //Parameters 
 //a:
 //b:
 //biasflag:
-
 //Description
 //This is an Octave function.
 
@@ -21,16 +19,16 @@ funcprot(0);
 
 rhs = argn(2)
 if(rhs<1 | rhs>3)
-error("Wrong number of input arguments.")
+error("Wrong number of input arguments.");
 end
 
 	select(rhs)
 	case 1 then
-	c = callOctave("xcorr2",a)
+	c = callOctave("xcorr2",a);
 	case 2 then
-	c = callOctave("xcorr2",a,b)
+	c = callOctave("xcorr2",a,b);
 	case 3 then
-	c = callOctave("xcorr2",a,b,biasflag)
-	end
+	c = callOctave("xcorr2",a,b,biasflag);
+	end;
 endfunction
 
diff --git a/macros/xcov1.sci b/macros/xcov1.sci
index 4c1f37e..b0bff93 100644
--- a/macros/xcov1.sci
+++ b/macros/xcov1.sci
@@ -1,4 +1,22 @@
 function [R,lag] = xcov1(X, Y, biasflag)
+// Compute covariance at various lags [=correlation(x-mean(x),y-mean(y))]. 
+//Calling Sequence
+//[R, lag] = xcov (X)
+//... = xcov (X, Y)
+//... = xcov (..., maxlag)
+//... = xcov (..., scale)
+//Parameters 
+//X: Input vector
+//Y: if specified, compute cross-covariance between X and Y, otherwise compute autocovariance of X.
+//maxlag: is specified, use lag range [-maxlag:maxlag], otherwise use range [-n+1:n-1].
+//scale:
+//     'biased': for covariance=raw/N,
+//     'unbiased': for covariance=raw/(N-|lag|),
+//     'coeff': for covariance=raw/(covariance at lag 0),
+//     'none': for covariance=raw
+//     'none': is the default. 
+//Description
+//Compute covariance at various lags [=correlation(x-mean(x),y-mean(y))]. Returns the covariance for each lag in the range, plus an optional vector of lags.
 
 funcprot(0);
 rhs = argn(2)
diff --git a/macros/yulewalker.sci b/macros/yulewalker.sci
index ebcb8aa..40fcadb 100644
--- a/macros/yulewalker.sci
+++ b/macros/yulewalker.sci
@@ -1,5 +1,13 @@
 function [A,V]= yulewalker(C)
-
+// Fit an AR (p)-model with Yule-Walker estimates given a vector C of autocovariances '[gamma_0, ..., gamma_p]'.
+//Calling Sequence
+//A    = yulewalker(C)
+//[A,V]= yulewalker(C)
+//Parameters 
+//C: Autocovariances
+//Description
+//Fit an AR (p)-model with Yule-Walker estimates given a vector C of autocovariances '[gamma_0, ..., gamma_p]'.
+//Returns the AR coefficients, A, and the variance of white noise, V.
 funcprot(0);
 lhs=argn(1);
 rhs= argn(2);