function g = etfe(z,M,N,T)
%ETFE	Computes the Empirical Transfer Function Estimate and Periodogram
%
%	G = etfe(Z)   or   G = etfe(Z,M)
%
%	Z contains the input-output data [y u] or a time series y. Only
%	single (or no) input systems can be handled. If an input is present
%	G is returned as the ETFE (the ratio of the output Fourier transform
%	to the input Fourier transform) for the data. For a time series G
%	is returned as the periodogram (the normed absolute square of the
%	Fourier transform) of the data. G is returned in the standard frequency
%	function format. See HELP FREQFUNC.
%
%	With M specified,a smoothing operation is performed on the raw spectral
%	estimates. The effect of M is then analogous to the same argument in
%	SPA.
%
%	The transfer function is estimated at 128 equally spaced frequencies
%	between 0 (excluded) and pi. This number can be changed to N (a power
%	of two) and the sampling interval can be changed from 1 (default) to T
%	by	G = etfe(Z,M,N,T).
%	The transfer function can be plotted by BODEPLOT. bodeplot(etfe(z))
%	is a possible construction. See also SPA.

%	L. Ljung 7-7-87
%	Revised 3-8-89, 4-20-91
%	Copyright (c) 1987-90 The MathWorks, Inc.
%	All Rights Reserved

%	Some initial checks on the arguments

if nargin<4,T=1;end, if T<0,T=1;end,if isempty(T),T=1;end
if nargin<3,N=128;end, if N<0, N=128;end,if isempty(N),N=128;end
if nargin<2,M=[];end
[Ncap,nz]=size(z);
if nz>2, error('This routine works only for single-input systems!'),return,end

% ** Add zeros to data so that it can support N frequency points **

if Ncap<2*N,z=[z;zeros(2*N-Ncap,nz)];end 

% ** Compute the Fourier transforms by FFT **
nfft = 2^nextpow2(max(Ncap,2*N)); % Corr 891007
if nz==1,Y=abs(fft(z(:,1),nfft)).^2;else Y=fft(z(:,1),nfft);end,l=length(Y);
if M<0,M=l;end,if isempty(M),M=l;end
M1=fix(l/M);sc=l/(2*N);
if fix(sc)~=sc, error('N must be a power of two!'),return,end
if M1>1,
ha = .54 - .46*cos(2*pi*(0:M1)'/M1);
ha=ha/(norm(ha)^2);Y=[Y(l-M1+2:l);Y];
	if nz==1,Y=filter(ha,1,Y);end
end
Yd=Y(M1+fix(M1/2)+sc:sc:M1+fix(M1/2)+l/2);
if nz==1,
	g(1,1:2)=[100 0];
	g(2:N+1,1)=(1:N)'*pi/T/N;%%%%%%Correction made here
	g(2:N+1,2)=T*Yd/Ncap;%Corr 910420
	return
end
if M1==1,clear Y,end
U=fft(z(:,2),nfft);
if M1>1,
	U=[U(l-M1+2:l);U];
	Y=filter(ha,1,Y.*conj(U));
	U=filter(ha,1,abs(U).^2);
	Yd=Y(M1+fix(M1/2)+sc:sc:M1+fix(M1/2)+l/2);
end
Ud=U(M1+fix(M1/2)+sc:sc:l/2+M1+fix(M1/2));clear U
g(1,1:3)=[101 1 21];
g(2:N+1,1)=(1:N)'*pi/N/T; %%%%%Correction made here
g(2:N+1,2)=abs(Yd./Ud);
g(2:N+1,3)=-180*phase((Yd./Ud)')'/pi;
