function ffplot(G,sd,C,mode)
%FFPLOT Plots frequency functions
%
%	ffplot(G)  or  ffplot(G,SD)
%
%	The syntax is identical to that of BODEPLOT
%	FFPLOT plots frequency functions with linear frequency scales
%	and uses Hz as the frequency unit.
%	BODEPLOT gives logarithmic frequency scales and uses rad/sec as the
%	frequency unit. 

%	L. Ljung 10-4-90,2-10-92
%	Copyright (c) 1990-92 by the MathWorks, Inc.
%	All Rights Reserved.

if nargin<4 , mode=[];end
if nargin<3, C=[];end
if nargin<2, sd=[];end
if isempty(mode),mode='sep';end
if isempty(C),C='B';end
if isempty(sd),sd=0;end

abp=['a';'A';'b';'B';'p';'P'];
if sum(sd(1)==abp)>0, if length(C)>1,mode=C;end,C=sd;sd=0;end
if length(sd)>1,mode=sd;sd=0;end

mode=mode(1:3);
if mode=='SEP', mode='sep';end
if mode=='SAM', mode='sam';end
if mode~='sam', mode='sep';end

if C == 'a', C='A'; end
if C == 'p', C='P'; end
if C == 'b', C='B'; end

[m,n]=size(G);
inputs1=G(1,:);mo=fix(inputs1/1000);mmo=max(mo);
kko=[]; % The output indices
for k=[0:mmo],if length(find(mo==k))>0,kko=[kko k];end,end
G=G(2:m,:);
for ko=kko+1

inputs=inputs1-(ko-1)*1000;
mi=max(inputs(find(inputs<19)));

kki=[]; % The input indices
for k=[1:mi 0], if length(find(inputs==k))>0,kki=[kki k];end,end
if sum(kki)==0 & C=='B', C='A';end
if mode=='sep'
  for k=kki
  if k==0 & C=='B',C='A';end
  if C=='B',subplot(211),else subplot(111),end
  if C=='A'| C=='B'

	semilogy(G(:,find(inputs==100+k))/2/pi,G(:,find(inputs==k)))
	if sd>0,hold on,
	ind=find(k+50==inputs);
	semilogy(G(:,ind-2)/2/pi,G(:,ind-1)+sd*G(:,ind),'-.')
	semilogy(G(:,ind-2)/2/pi,G(:,ind-1)-sd*G(:,ind),'-.') %Change here 2.04
	hold off,end
	xlabel('frequency (Hz)')

	if k~=0,title(['AMPLITUDE PLOT, input # ',int2str(k),' output # ',int2str(ko)])
	else title(['SPECTRUM output # ',int2str(ko)]),end
  end
  if ((C=='B')|(C=='P'))&k~=0
        if (C=='B'), subplot(212), end
        if (C=='P'),subplot(111),end

	plot(G(:,find(inputs==100+k))/2/pi,G(:,find(inputs==20+k)))
	if sd>0,hold on,
	ind=find(k+70==inputs);
	plot(G(:,ind-4)/2/pi,G(:,ind-1)+G(:,ind)*sd,'-.')
	plot(G(:,ind-4)/2/pi,G(:,ind-1)-sd*G(:,ind),'-.') %Change here 2.04
	hold off,end
	title(['PHASE PLOT, input # ',int2str(k),' output # ',int2str(ko)])
	xlabel('frequency (Hz)')
	ylabel('phase')
  end


  if k~=kki(length(kki)),pause,end
end
end
if mode=='sam' 

	if C=='B',subplot(211),else subplot(111),end
      if C=='B' | C=='A'	
	semilogy(G(:,find(inputs>99&inputs<120))/2/pi,G(:,find(inputs<20&inputs>-1)))
	if sd>0,hold on,
	ind=find(inputs<69 & inputs>49);
	semilogy(G(:,ind-2)/2/pi,G(:,ind-1)+G(:,ind)*sd,'-.')
	semilogy(G(:,ind-2)/2/pi,G(:,ind-1)-sd*G(:,ind),'-.'),hold off,end
	xlabel('frequency (Hz)')
	title(['AMPLITUDE PLOT output #',int2str(ko)])
      end
      if C=='B'| C=='P'
         if C=='B',subplot(212),else subplot(111),end
        plot(G(:,find(inputs>100&inputs<120))/2/pi,G(:,find(inputs<39 & inputs>19)))
	if sd>0,hold on
	ind=find(inputs<89 & inputs>69);
	plot(G(:,ind-4)/2/pi,G(:,ind-1)+G(:,ind)*sd,'-.')
	plot(G(:,ind-4)/2/pi,G(:,ind-1)-sd*G(:,ind),'-.'),hold off,end
	xlabel('frequency (Hz)')
	ylabel('phase')
	title(['PHASE PLOT output #',int2str(ko)])
      end
	
end,if ko~=kko(length(kko))+1,pause,end

end

hold off
set(gcf,'NextPlot','replace');

