17 files changed, 323 insertions, 12 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_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/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/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/detrend1.sci b/macros/detrend1.sci
index b9929f4..52ac26f 100644
--- a/macros/detrend1.sci
+++ b/macros/detrend1.sci
@@ -1,4 +1,15 @@
 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");
diff --git a/macros/diffpara.sci b/macros/diffpara.sci
index 3f9e2d9..c36ae7c 100644
--- a/macros/diffpara.sci
+++ b/macros/diffpara.sci
@@ -1,4 +1,21 @@
 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)
diff --git a/macros/durbinlevinson.sci b/macros/durbinlevinson.sci
index cd56e48..74dbb46 100644
--- a/macros/durbinlevinson.sci
+++ b/macros/durbinlevinson.sci
@@ -1,4 +1,17 @@
 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");
diff --git a/macros/fftshift1.sci b/macros/fftshift1.sci
index 2928850..d8721cd 100644
--- a/macros/fftshift1.sci
+++ b/macros/fftshift1.sci
@@ -1,4 +1,25 @@
 function y= fftshift1(X,DIM)
+//Perform a shift of the vector X, for use with the 'fft' and 'ifft' 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 'fft' and 'ifft' 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 (fft (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');
diff --git a/macros/fractdiff.sci b/macros/fractdiff.sci
index 1731dc1..979a079 100644
--- a/macros/fractdiff.sci
+++ b/macros/fractdiff.sci
@@ -1,4 +1,10 @@
 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)
diff --git a/macros/hamming.sci b/macros/hamming.sci
index 6fa9a12..49a970f 100644
--- a/macros/hamming.sci
+++ b/macros/hamming.sci
@@ -1,4 +1,14 @@
 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)
diff --git a/macros/hanning.sci b/macros/hanning.sci
index d90faf2..60ca783 100644
--- a/macros/hanning.sci
+++ b/macros/hanning.sci
@@ -1,5 +1,14 @@
 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);
diff --git a/macros/hilbert1.sci b/macros/hilbert1.sci
index 1c0fe0e..fbcb136 100644
--- a/macros/hilbert1.sci
+++ b/macros/hilbert1.sci
@@ -9,7 +9,7 @@ function h= hilbert1(f, varargin)
 //N: The result will have length N
 //dim : It analyses the input in this dimension
 //Description
-//h = hilbert (f) computes the extension of the real valued signal f to an analytic signal. If f is a matrix, the transformation is applied to each column. For N-D arrays, the transformation is applied to the first non-singleton dimension.
+//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.
 //
diff --git a/macros/ifftshift1.sci b/macros/ifftshift1.sci
index 6c5ae8e..7426130 100644
--- a/macros/ifftshift1.sci
+++ b/macros/ifftshift1.sci
@@ -1,4 +1,16 @@
 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');
diff --git a/macros/lib b/macros/lib
index 4261b2e..4da3137 100644
--- a/macros/lib
+++ b/macros/lib
diff --git a/macros/names b/macros/names
index fb6e9a7..001dba2 100644
--- a/macros/names
+++ b/macros/names
@@ -1,7 +1,9 @@
 ac2poly
 ac2rc
 arParEst
+ar_psd
 arburg
+arch_test
 arch_fit
 arcov
 armcov
@@ -11,6 +13,7 @@ barthannwin
 bartlett
 besselap
 besself
+bilinear
 bitrevorder
 blackman
 blackmanharris
@@ -36,6 +39,7 @@ check
 chirp
 clustersegment
 cmorwavf
+cohere
 convmtx
 corrmtx
 cplxreal
@@ -228,6 +232,7 @@ synthesis
 tf2sos
 tf2zp
 tf2zpk
+tfe
 transpose
 trial_iirlp2mb
 triang
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/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