function [t,u,N,RMSu,Au,Bu,Cu,Du] = inputs( INPUT, r, dt, N, tau )
% [t,u,N,RMSu,Au,Bu,Cu,Du] = inputs( INPUT, r, dt, N, tau  )
%
%    INPUT           DESCRIPTION
% -----------       ------------------------------------------------------
%  INPUT             'IMPULSE', 'CHIRP', 'POISSON', 'FILTER' 'FLICKER'
%      r              number of independent inputs
%     dt              time step increment
%      N              number of time steps
%    tau              characteristic time scale of the random process
%                     refer to comments within each case
%
%   OUTPUT           DESCRIPTION
% -----------       ------------------------------------------------------
%      t             time axis data                           [ 1 x N ]
%      u             input time series                        [ r x N ]
%      N             actual number of time steps
%     RMSu           root mean square of u               
%  Au,Bu,Cu,Du       continuous time state space model for noise process
%
% H.P. Gavin, Duke Univ, 2019-04-09

  Au = 0; Bu = 0; Cu = 0; Du = 0;
  t = [0:N-1]*dt; 

  switch INPUT

  case 'IMPULSE'  % ===== "tau" is the duration of the impuse
     u = zeros(r,N);
     n = round(tau/dt); 
     tau = n*dt;
     u(:,1:n+1) = ones(r,1) * (1 - cos(2*pi*[0:n]*dt/tau)) * sqrt(N*dt)/tau;
     RMSu = 0.5;

   case 'CHIRP'    % ===== "tau" is the longest period in the chirp 
     Ao = 10;  Af = 2; Fo = 1/tau; Ff = 10/tau; CYCLES = 40, P = 1.5; N = 0.5;
     [aa,vv,dd,ns,t] = chirp ( Ao, Af, Fo, Ff, CYCLES, t, P, N);
     u = zeros(r,ns+11*r);
     for ii=1:r
       u(ii,10*ii+1:10*ii+ns) = vv; 
     end 
     RMSu = sqrt(var(u'));
     N = size(u,2);
     t = [0:N-1]*dt; 
   
   case 'POISSON'  % ===== "tau" is the return period of the Poison process
%    u = Poisson_rnd(tau/dt,r,N).*(2*(rand(r,N)>0.5)-1)*sqrt(tau)/dt;
     u = poissrnd(dt/tau,r,N).*(2*(rand(r,N)>0.5)-1)*sqrt(tau)/dt;
     RMSu = sqrt(1/dt);                % ... for Bi-Polar Poisson noise.   

   case 'FILTER'   % ===== "tau" is the central period of the process

     wn = 2*pi/tau*(1+0.1*randn(1,r));    % modify to fit your purpose
     zf = 0.5*(1+0.1*randn(1,r));       % modify to fit your purpose

     Au = zeros(2*r);
     Bu = zeros(2*r,r);
     Cu = zeros(r,2*r);
     Du = zeros(r,r);

     for ii=1:r
       Ai = [ 0 1 ; -wn(ii)^2   -2*wn(ii)*zf(ii) ];
       Au( 2*ii-1:2*ii , 2*ii-1:2*ii ) = Ai;
       Bu(2*ii-1:2*ii,ii) = [ 0 ; 1 ];
       Cu(ii,2*ii-1:2*ii) = [ 0  2*wn(ii)*zf(ii) ];
     end

     u = lsym(Au,Bu,Cu,Du,randn(r,N)/sqrt(dt),t);  % filtered white noise

     RMSu = sqrt(diag(real(Cu*liap(Au,Bu*Bu')*Cu')));   % RMS of each input 

   otherwise % 'FLICKER' % ===== "tau" is the cut-off period of the "1/f" noise
     Au = -(1/tau)*eye(r);  Bu = sqrt(2/tau)*eye(r); Cu = eye(r); Du = zeros(r);
     u = lsym(Au,Bu,Cu,Du,randn(r,N)/sqrt(dt),t);       % filtered noise
     RMSu = sqrt(diag(real(Cu*liap(Au,Bu*Bu')*Cu')));   % RMS of each input 

  end % =====================================================================
