%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: noiseSpectrum makes spectra of ssm outputs from the given noises spectra.
%
% CALL: [ao_out] = noiseSpectrum(sys, plist_inputs)
%
% INPUTS:
%         - sys, (array of) ssm object
%         - plist_inputs contains names of noise inputs, spectra of
%         specified noise inputs,  returned
%         outputs and frequency's vector of returned spectrum
%
% OUTPUTS:
%          _ ao_out contains spectra of specified outputs and states
%
% <a href="matlab:utils.helper.displayMethodInfo('ssm', 'noiseSpectrum')">Parameters Description</a>
%
% VERSION: $Id: noiseSpectrum.m,v 1.8 2011/04/08 08:56:23 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function varargout = noiseSpectrum(varargin)
  
  %% starting initial checks
  error('### this function is deprecated, use ssm/PSD or ssm/CPSD instead')
  
  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end
  
  utils.helper.msg(utils.const.msg.MNAME, ['running ', mfilename]);
  
  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end
  
  % Collect all SSMs and plists
  [sys, ssm_invars, rest] = utils.helper.collect_objects(varargin(:), 'ssm', in_names);
  [pl, invars2, rest]  = utils.helper.collect_objects(rest(:), 'plist');
  pl = combine(pl, plist(rest{:}) );
  plist_inputs = combine(pl, getDefaultPlist());
  
  
  %% retrieve system infos
  if numel(sys)~=1
    error('noisespectrum needs exactly one ssm as an input')
  end
  if ~sys.isnumerical
    error(['error because system ',sys.name,' is not numerical']);
  end
  
  %% collecting data
  noise_varnames = find(plist_inputs, 'noise variable names');
  if ischar(noise_varnames)
    noise_varnames = {noise_varnames};
  end
  Nnoise = numel(noise_varnames);
  In_noise_spectrum = find(plist_inputs, 'spectrum');
  freq_vector = find(plist_inputs, 'f');
  return_outputs = find(plist_inputs, 'return outputs');
  if ischar(return_outputs)
    return_outputs = {return_outputs};
  end
  Noutputs=numel(return_outputs);
  
  %% verify frequency
  Nspectrum=length(In_noise_spectrum);
  ur=1;
  diff_vec=[];
  vec=[];
  for u=1: Nspectrum
    noise=In_noise_spectrum{u};
    if isa(noise,'ao')==1
      vec(ur,:)=noise.x;
      ur=ur+1;
    end
  end
  
  if ~isempty(vec) && ur>=3
    for ui=1: ur-2
      diff_vec(ui,:)= vec(ui+1,:)-vec(ui,:);
    end
  end
  
  if ~isempty(diff_vec)
    if  any(diff_vec~=0)
      freq_status=0;
      error('### Please specify inputs in the same frequency range .')
    else
      freq_status=1;
    end
  end
  
  if isempty(vec) && isempty(freq_vector)
    f1 = find(plist_inputs,'f1');
    f2 = find(plist_inputs,'f2');
    nf = find(plist_inputs,'nf');
    if isempty(f1) || isempty(f2)|| isempty(nf)
      error('### Please specify frequency vector a start and stop frequency .');
    else
      freq_vector = linspace(f1, f2, nf);
    end
    
  end
  
  if ~isempty(freq_vector)
    if ur==2 || (ur>=3 && freq_status==1)
      if length(freq_vector) ~=length(vec(1,:)) || any((freq_vector - vec(1,:))~=0)
        error('### Please specify inputs in asked frequency vector .');
      else
        freq_vector=vec(1,:);
      end
    end
  else
    if ur==2 || (ur>=3 && freq_status==1)
      freq_vector=vec(1,:);
    end
  end
  
  %% spectra inputs
  if Nnoise>=1
    for u=1: Nnoise
      noise=In_noise_spectrum{u};
      if isa(noise,'ao')==1
        noise_uncor{u}= In_noise_spectrum{u}.y';
      elseif isa(noise,'double')==1
        noise_uncor{u}=In_noise_spectrum{u}.*ones(1,length(freq_vector));
      elseif isa(noise,'miir')==1
        rep=resp(In_noise_spectrum{u},plist('f',freq_vector));
        noise_uncor{u}= rep.y';
      elseif isa(noise,'pzmodel')==1
        rep=resp(In_noise_spectrum{u},plist('f',freq_vector));
        noise_uncor{u}=rep.y';
      else error('spectrum should be an Ao, a double, a miir or a pzmodel')
      end
    end
  else
    error('Specify at least one input noise spectrum')
  end
  
  %% SSM Transfer function
  plist_OutputTF=plist('inputs', noise_varnames,'outputs', return_outputs, 'f',freq_vector);
  if ~isempty(ver('control')) % check that control system toolbox is installed
    rep_TF= sys.bodecst(plist_OutputTF);
  else
    rep_TF= sys.bode(plist_OutputTF);
  end
  
  [Ninp Nout]=size(rep_TF);
  if (Ninp~=Nnoise)||(Nout~=Noutputs)
    error('Asked input /output does not exist. Please verify names of return outputs and noise variable names')
  end
  
  %% Spectral density
  dd=[noise_uncor{1,:}];
  noise=reshape(dd,length(freq_vector),Nnoise);
  for u=1:length(freq_vector)
    mat_TF=reshape(rep_TF.y(u),Noutputs,Nnoise);
    for ui=1:Noutputs
      mat_TFN(ui,:)=  mat_TF(ui,:).* abs(noise(u,:));
    end
    resu(:,:,u)=mat_TFN*conj(transpose(mat_TF));
  end
  
  %% saving in aos
  
  ao_out = ao.initObjectWithSize(Noutputs,Noutputs);
  for i=1:Noutputs
    
    for j=1:Noutputs
      if j==i
        nameao=[sprintf('PSD(%s)',  return_outputs{i})];
      else
        nameao=[sprintf('CPD(%s)',  [return_outputs{i} ',' return_outputs{j}])];
      end
      for u=1:length(freq_vector)
        ao_out_y(u)= resu(j,i,u);
      end
      ao_out(i,j) = ao(fsdata(freq_vector,ao_out_y));
      ao_out(i,j).setName(nameao);
      ao_out(i,j).setXunits('Hz');
      ao_out(i,j).setYunits('arb');
      ao_out(i,j).setDescription(...
        ['Cross Power Spectrum of ' ,  return_outputs{j} , return_outputs{i}]);
    end
    if nargout == numel(ao_out)
      for ii = 1:numel(ao_out)
        varargout{ii} = ao_out(ii);
      end
    else
      varargout{1} = ao_out;
    end
  end
  varargout  = {ao_out };
end


%--------------------------------------------------------------------------
% Get Info Object
%--------------------------------------------------------------------------
function ii = getInfo(varargin)
  
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pl   = [];
  else
    sets = {'Default'};
    pl   = getDefaultPlist;
  end
  % Build info object
  ii = minfo(mfilename, 'ssm', 'ltpda', utils.const.categories.op, '$Id: noiseSpectrum.m,v 1.8 2011/04/08 08:56:23 hewitson Exp $', sets, pl);
  
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------
function pl = getDefaultPlist()
  pl = plist();
  
  p = param({'noise variable names', 'A cell-array of strings specifying the desired input variable names.'}, {} );
  pl.append(p);
  
  p = param({'spectrum', 'The spectrum of this noise between input ports for the <i>time-continuous</i> noise model.'}, []);
  pl.append(p);
  
  p = param({'f', 'frequency vector of output spectrum .'}, []);
  pl.append(p);
  
  p = param({'return outputs', 'A cell-array of output ports to return.'}, {});
  pl.append(p);
  
end