--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ssm/noiseSpectrum.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,238 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% 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
+