function y = utSSFastFsim(A,B,C,D,Ts,w,cw,u,X0,InpDel)
% Frequency output calculation, typically performed by fsim.
% Optimized for speed.
% 
% Normal route: y = idpack.ltidata.fsim(sys, w, u, X0);

% Note: InpDel exists only when delays are estimated in (CT, CT) case. cw
% is complex frequency (iw or exp(jwT)) scaled by 1/sqrt(nw); w is
% frequency vector X0 must be a column vector.
%
% X0 are time domain initial conditions. 

% Copyright 1986-2011 The MathWorks, Inc.

[ny,nu] = size(D);
if any(X0)
   B = [B, X0]; D = [D, zeros(ny,1)];
   if Ts==0
      u = [u, ones(size(cw))]; 
   else
      u = [u, cw];
   end
end
[A,B,C] = ltipack.xscale(A,B,C,D,[],Ts,'Warn',false);
h = ssfresp(A,B,C,D,[],cw);
[nr, nc, N] = size(h);

if any(InpDel)
   h(:,1:nu,:) = idpack.utAddDelayFresp(h(:,1:nu,:),InpDel,[],w,Ts);
end

if nc>1
   y = reshape((sum((reshape(shiftdim(h,1),nc,nr*N)).*(repmat(u.',1,nr)))).',N,nr);
else
   y = reshape(((reshape(shiftdim(h,1),nc,nr*N)).*(repmat(u.',1,nr))).',N,nr);
end
