Global Index (short | long) | Local contents | Local Index (short | long)
[G,PV]=th2ff(th,nnu,w,nny)
TH2FF Computes a model's frequency function,along with its standard deviation G = TH2FF(TH) or [G,NSP] = TH2FF(TH) TH: A matrix defining a model, as described (see also) 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 also FREQFUNC). If TH describes a time series, G is returned as its spectrum. Both discrete and continuous time models are handled. If the model TH has several inputs, G will be returned as the transfer functions of selected inputs # j1 j2 .. jk by G = TH2FF(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 (a row vector, gene- rated e.g. by LOGSPACE) by G = TH2FF(TH,ku,w). The transfer function can be plotted by BODEPLOT. If the model TH has several outputs, G will be returned as the frequency function at selected outputs ky (a row vector) by G=TH2FF(TH,ku,w,ky); (Default is all outputs). See also BODEPLOT, FFPLOT, FREQFUNC, NYPLOT and TRF.
| This function calls | This function is called by |
|---|---|
function [G,PV]=th2ff(th,nnu,w,nny)
% L. Ljung 7-7-87,1-25-92
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 2.3 $ $Date: 1997/12/02 03:40:39 $
if nargin<1
disp('Usage: G = TH2FF(TH)')
disp(' [G,NSP] = TH2FF(TH,INPUTS,FREQUENCIES,OUTPUTS)')
return
end
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
[dum,P,dum]=th2par(th);
if any(any(isinf(P))),eval('[G,PV]=trf(th,nnu,w,nny);'),return,end
if isempty(P)|norm(P)==0
eval('[G,PV]=trf(th,nnu,w,nny);')
return
end
if isthss(th)
if length(th2par(th))>10,
disp(str2mat('This may take a while.',...
'TRF will be faster, but does not evaluate standard deviations.'))
end
if nargout==1
eval('G=trfsssd(th,nnu,nny,w);')
else eval('[G,PV]=trfsssd(th,nnu,nny,w);')
end
return
end
if length(nnu)==1, if nnu==0,nnu=[];end,if nnu<0, nnu=1:th(1,3);end,end
% *** Compute the model polynomials and form the basic
% matrix of frequencies ***
nu=th(1,3);
[a,b,c,d,f]=th2poly(th);[par,P]=th2par(th);if isempty(P),isP=0;else isP=1;end
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+4)=[100+k,k,k+50,k+20,k+70];
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+3)=phase(GC)'*180/pi;
if T<0, G(2:nrc,sc+3)=G(2:nrc,sc+3)-w'*delays(k)*180/pi;end
ll=[1:na,na+ks+1:na+ks+nb(k),na+ks+nb(k)+nc+nd+kf+1:na+ks+nb(k)+nc+nd+kf+nf(k)];
ks=ks+nb(k);kf=kf+nf(k);
if isP,
P=th(3+ll,ll);
if nb(k)==0
dga=zeros(length(OM),1);dgp=dga;
else
[dga,dgp]=ffsdcal(a,b(k,1:nb(k)+nk(k)),f(k,1:nf(k)+1),na,...
nb(k),nf(k),nk(k),GC,OM,P,T);
end
G(2:nrc,[sc+2 sc+4])=[dga,dgp];
end
sc=sc+5;
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(H).^2)'*th(1,1)*abs(T);
%
% *** Compute the standard deviation of the spectrum ***
if isP,Ncap=getncap(th);
if ~isempty(Ncap)
ll=[1:na,na+sum(nb)+1:na+sum(nb)+nc+nd];
P=th(3+ll,ll);
dga=ffsdcal(a,c,d,na,nc,nd,0,H,OM,P,T);
PV(2:nrc,3)=sqrt(2/Ncap)*PV(2:nrc,2)+...
2*abs(T)*th(1,1)*(dga.*abs(H)');
end
end
if isempty(nnu), G=PV;end