Documentation of etfe

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

g = etfe(z,M,N,T)

Help text

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 also 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 BODEPLOT, FFPLOT and SPA.

Cross-Reference Information

Listing of function etfe

function g = etfe(z,M,N,T)

%   L. Ljung 7-7-87
%   Revised 3-8-89, 4-20-91
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 2.3 $  $Date: 1997/12/02 03:44:16 $

if nargin < 1
   disp('Usage: G = ETFE(Z)')
   disp('       G = ETFE(Z,FREQ_RESOL,No_OF_FREQUENCIES,T)')
   return
end
%       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;