function [G,PV]=trf(th,nnu,w,nny)
%TRF	Computes a model's frequency function.
%
%	G = trf(TH)  or  [G,NSP] = trf(TH)
%
%	TH: A matrix defining a model, as described in HELP THETA
%	G is returned as the transfer function estimate, and NSP (if specified)
%	as the noise spectrum, corresponding to the model TH. These matrices
%	contain also estimated standard deviations, calculated from the
%	covariance matrix in TH, and are  of the standard frequency function
%	format (see HELP FREQFUNC).If TH describes a time series, G is returned
%	as its spectrum.
%	
%	If the model TH has several inputs, G will be returned as the transfer
%	functions of selected inputs # j1 j2 .. jk by 
%	G = trf(TH,[j1 j2 ... jk])  [default is all inputs]. The functions 
%	are computed at 128 equally spaced frequency-values between 0(excluded)
%	and pi/T, where T is the sampling interval specified by TH. The func-
%	tions can be computed at arbitrary frequencies w RAD/S, (a row vector,
%	generated e.g. by LOGSPACE) by G = trf(TH,ku,w). The transfer function
%	can be plotted by BODEPLOT. bodeplot(trf(TH)) is a possible 
%	construction. If the model TH has several outputs, G will be returned
%	as the frequency function at selected outputs ky (a row vector) by
%	G=trf(TH,ku,w,ky); (Default is all outputs). See also TH2FF.

%	L. Ljung 7-7-87, 1-25-92
%	Copyright (c) 1987-92 by the MathWorks, Inc.
%	All Rights Reserved.
 
T=gett(th);

% *** Set up default values ***
if nargin<4, nny=[];end
if nargin<3, w=[];end
if nargin<2 , nnu=[];end
if T>0,wdef=pi*[1:128]/128/T;
else wdef=logspace(log10(pi/abs(T)/100),log10(10*pi/abs(T)),128);
end
if isempty(w),w=wdef;end
if length(w)==1, if w<0, w=wdef;end,end
[nrw,ncw]=size(w);if nrw>ncw, w=w.';end
if T>0 & any(w>pi/T),disp('WARNING: Frequencies in w exceed the Nyquist frequency!'),end

if isempty(nnu), nnu=1:th(1,3);end
if length(nnu)==1,if nnu==0,nnu=[];end,if nnu<0,nnu=1:th(1,3);end,end
if isthss(th), if nargout==1, eval('G=trfss(th,nnu,nny,w);')
	     else eval('[G,PV]=trfss(th,nnu,nny,w);'),end,return,end
% *** Compute the model polynomials and form the basic
%      matrix of frequencies ***
nu=th(1,3);
[a,b,c,d,f]=th2poly(th);
na=th(1,4);nb=th(1,5:4+nu);nc=th(1,5+nu);nd=th(1,6+nu);
nf=th(1,7+nu:6+2*nu);nk=th(1,7+2*nu:6+3*nu);nnd=na+sum(nb)+nc+nd+sum(nf);
nm=max([length(a)+length(f)-1 length(b) length(c) length(d)+length(a)-1]);
   i=sqrt(-1);
if T>0,OM=exp(-i*[0:nm-1]'*w*T);end
if T<=0,
OM=ones(1,length(w));
for kom=1:nm-1
OM=[OM;(i*w).^kom];
end
[nrth,ncth]=size(th);
if nrth>3+nnd,delays=th(nrth,1:length(nk));else delays=zeros(1,length(nk));end
end

% *** Compute the transfer function(s) GC=B/AF ***

sc=1;kf=0;ks=0;nrc=length(w)+1;
for k=nnu
    G(1,sc:sc+2)=[100+k,k,k+20];	
    G(2:nrc,sc)=w';
     gn=conv(a,f(k,:));
if T>0,indb=1:length(b(k,:));indg=1:length(gn);
else indb=length(b(k,:)):-1:1; indg=length(gn):-1:1;end
     GC=(b(k,:)*OM(indb,:))./(gn*OM(indg,:));
    G(2:nrc,sc+1)=abs(GC');
    G(2:nrc,sc+2)=phase(GC)'*180/pi;
    if T<0, G(2:nrc,sc+2)=G(2:nrc,sc+2)-w'*delays(k)*180/pi;end
sc=sc+3;
end
if nargout==1 & nu>0, return, end

% *** Compute the transfer function H=C/AD *** 

hn=conv(a,d);
if T>0,indc=1:length(c);indh=1:length(hn);
else indc=length(c):-1:1; indh=length(hn):-1:1;end
H=(c*OM(indc,:))./(hn*OM(indh,:));
PV(1,1:3)=[100,0,50];
PV(2:nrc,1)=w';
PV(2:nrc,2)=abs(T)*(abs(H).^2)'*th(1,1);

if isempty(nnu) G=PV;end


