24 files changed, 580 insertions, 47 deletions

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/autoreg_matrix.sci b/macros/autoreg_matrix.sci
index bc36969..8778a84 100644
--- a/macros/autoreg_matrix.sci
+++ b/macros/autoreg_matrix.sci
@@ -1,5 +1,22 @@
 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)
diff --git a/macros/cceps.sci b/macros/cceps.sci
index 8360085..3756a5d 100644
--- a/macros/cceps.sci
+++ b/macros/cceps.sci
@@ -3,10 +3,19 @@ function y = cceps (x,correct)
 //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/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/dwt.sci b/macros/dwt.sci
index 60dc7f6..b564c7c 100644
--- a/macros/dwt.sci
+++ b/macros/dwt.sci
@@ -1,5 +1,22 @@
 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)
diff --git a/macros/fftconv.sci b/macros/fftconv.sci
index 4ec2815..d39441d 100644
--- a/macros/fftconv.sci
+++ b/macros/fftconv.sci
@@ -1,5 +1,16 @@
 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)
@@ -10,6 +21,6 @@ end
 	case 2 then
 	y = callOctave("fftconv", X, Y);
 	case 3 then
-	y = callOctave("ifftn",X, Y, varargin(1));
+	y = callOctave("fftconv",X, Y, varargin(1));
 	end
 endfunction
diff --git a/macros/fftn.sci b/macros/fftn.sci
index ba7f352..b8b4170 100644
--- a/macros/fftn.sci
+++ b/macros/fftn.sci
@@ -1,5 +1,16 @@
 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)
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
index 2a17374..3c32928 100644
--- a/macros/filter1.sci
+++ b/macros/filter1.sci
@@ -1,5 +1,38 @@
 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)
diff --git a/macros/filter2.sci b/macros/filter2.sci
index 82f3f0b..b8b48e4 100644
--- a/macros/filter2.sci
+++ b/macros/filter2.sci
@@ -1,5 +1,20 @@
 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)
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
index e3ff152..70c95b7 100644
--- a/macros/fir1.sci
+++ b/macros/fir1.sci
@@ -1,7 +1,27 @@
 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)
+rhs = argn(2);
 if(rhs<2 | rhs>5)
 error("Wrong number of input arguments.");
 end
@@ -9,11 +29,11 @@ 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));
+        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
index 418a9e4..990fb3e 100644
--- a/macros/fir2.sci
+++ b/macros/fir2.sci
@@ -1,7 +1,29 @@
 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)
+rhs = argn(2);
 if(rhs<3 | rhs>6)
 error("Wrong number of input arguments.");
 end
@@ -9,11 +31,11 @@ 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));
+        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/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/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/ifftn.sci b/macros/ifftn.sci
index 269a77d..3d26c04 100644
--- a/macros/ifftn.sci
+++ b/macros/ifftn.sci
@@ -1,5 +1,17 @@
 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)
diff --git a/macros/lib b/macros/lib
index d265058..4261b2e 100644
--- a/macros/lib
+++ b/macros/lib
diff --git a/macros/names b/macros/names
index c6a8d88..fb6e9a7 100644
--- a/macros/names
+++ b/macros/names
@@ -2,6 +2,7 @@ ac2poly
 ac2rc
 arParEst
 arburg
+arch_fit
 arcov
 armcov
 aryule
@@ -38,6 +39,7 @@ cmorwavf
 convmtx
 corrmtx
 cplxreal
+cpsd
 cummax
 cummin
 czt
@@ -65,6 +67,7 @@ fft2
 fftconv
 fftfilt
 fftn
+fftshift1
 fht
 filter1
 filter2
@@ -72,6 +75,7 @@ filternorm
 filtfilt
 filtic
 filtord
+findpeaks
 fir1
 fir2
 firpmord
@@ -80,6 +84,7 @@ flattopwin
 fracshift
 fractdiff
 freqs
+freqz
 fwhm
 fwhmjlt
 fwht
@@ -102,6 +107,7 @@ idct2
 idst1
 ifft
 ifft2
+ifftn
 ifftshift1
 ifht
 ifwht
@@ -170,6 +176,7 @@ pulseperiod
 pulsesep
 pulsewidth
 pulstran
+pwelch
 rc2ac
 rc2is
 rc2lar
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/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/unwrap2.sci b/macros/unwrap2.sci
index 5b19283..b8ea9de 100644
--- a/macros/unwrap2.sci
+++ b/macros/unwrap2.sci
@@ -1,20 +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)
+lhs = argn(1);
+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
-        Y = callOctave("unwrap",X)
+        Y = callOctave("unwrap",X);
 	case 2 then
-		Y = callOctave("unwrap",X,TOL)
-    case 3 then
-        Y = callOctave("unwrap",X,TOL,DIM)  
-    end
+	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