function [mag,ph,w,FocusInfo] = freqresp(D,Grade,wspec,isPlotted)
% Generates frequency response data (magnitude+phase) for
% MIMO LTI models. Used by BODE, NICHOLS, and NYQUIST.

%   Copyright 2020 The MathWorks, Inc.
Ts = abs(D.Ts);
[w,Focus] = parseFreqSpec(D,wspec,isPlotted && isreal(D)); 
FocusInfo = struct('Focus',Focus,'DynRange',Focus,'Soft',false);
         
% Gather external delays
[ny,nu] = iosize(D);
linDelay = D.Delay.Input(:,ones(1,ny)).' + D.Delay.Output(:,ones(1,nu));
if Ts>0
   linDelay = linDelay * Ts;  % equivalent continuous-time delays
end

% Compute response without external delays
D.Delay.Input(:) = 0;   D.Delay.Output(:) = 0;
[w,is] = sort(w);  % sort frequencies (required for phase unwrapping)
h = permute(fresp(D,w),[3 1 2]);
[mag,ph] = ltipack.getMagPhase(h,1,isPlotted);
   
% Add phase shift due to external delays
if norm(linDelay,1)>0
   ph = ph - reshape(w*reshape(linDelay,[1 numel(linDelay)]),size(ph));
end

% Undo frequency sorting for [mag,ph] = xxx(sys,w) syntax
if ~isPlotted && any(diff(is)<0)
   w(is) = w;    mag(is,:,:) = mag;   ph(is,:,:) = ph;
end

