function [y,t] = fimpulse(w,f,iu)
%FIMPULSE Impulse response from freqency response (MIMO).
%
%        [y,t] = fimpulse(w,f,iu)
%            y = fimpulse(w,f,iu)
%
% f is MVFR matrix, iu is the input on which a unit
% impulse is applied.
% w is the frequency list, which must have (N/2)+1
% evenly spaced points on it. The first point should
% be 0. N should be a power of 2. w is most conveniently 
% generated using function T2W.
%
% The output y is a matrix with as many columns as there are
% outputs, and N rows, each representing one time point.
%
% Optionally, t is the list of time points at which
% the impulse response is computed (it can also be
% generated by T2W).
%
% This function allows the impulse response of systems
% containing time delays to be generated.
%
% See also: FSTEP

% J.M.Maciejowski, 1 June 1988. Modified 3 and 29 December 1988 by JMM.
% Copyright (c) 1988, by Cambridge Control Ltd.
% MRN0037

[nout,nin] = fsize(w,f);
lw = length(w);
N = 2*(lw-1);

% Check whether frequency list looks OK:
wok = 1;
if w(1)~=0, wok=0; w1flag=1; end  % 1st element not 0
M = round(log(N)/log(2)); 
if abs(N-(2^M)) > 0.1, wok=0; lwflag=1; end  % Length of w incorrect
wdiff = diff(w);
if max(wdiff)-min(wdiff) > 1000*eps*lw,  % Frequencies not
  wok= 0;                                      % uniformly spaced
  wdflag = 1;
end 
if ~wok,     % Something is wrong
  disp('Warning: Frequency list is not as expected')
  disp('         for the following reasons:')
  if w1flag, disp('           First frequency in list is not 0'); end
  if lwflag, disp('           Length of frequency list incorrect'); end
  if wdflag, disp('           Frequencies not uniformly spaced'); end
  disp('         Use T2W or see Reference Manual.  (FIMPULSE)');
end  % of checking frequency list

%First get the right input:
f = fget(w,f,[(1:nout)',iu*ones(nout,1)]);

% Make the response at the highest frequency real:
f(lw) = abs(f(lw));  % This assumes that abs(f(lw)) is so small that
                     % this change will have negligible effect.

% Put conjugate response on the lower half of unit circle:
ff = [f;conj(f(lw-1:-1:2,:))];

% Now get the response:
y = real(ifft(ff));   % ifft(ff) should be real anyway.

% Now the time vector:
T = 2*pi/w(2);  % Simulation horizon
dt = T/(N-1);   % Interval between successive time points
if nargout > 1,
  t = 0:dt:T;     % Time vector
end

% Now scale the response correctly:
y = y/dt;