--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/spsd.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,221 @@
+% SPSD implements the smoothed (binned) PSD algorithm for analysis objects.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: SPSD implements the smoothed PSD algorithm for analysis objects.
+%
+% CALL:        bs = spsd(a1,a2,a3,...,pl)
+%              bs = spsd(as,pl)
+%              bs = as.spsd(pl)
+%
+% INPUTS:      aN   - input analysis objects
+%              as   - input analysis objects array
+%              pl   - input parameter list
+%
+% OUTPUTS:     bs   - array of analysis objects, one for each input
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'spsd')">Parameters Description</a>
+%
+% VERSION:     $Id: spsd.m,v 1.21 2011/07/11 10:43:35 adrien Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = spsd(varargin)
+
+  import utils.const.*
+
+  % Check if this is a call for parameters
+  if utils.helper.isinfocall(varargin{:})
+    varargout{1} = getInfo(varargin{3});
+    return
+  end
+
+  utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
+
+  % Collect input variable names
+  in_names = cell(size(varargin));
+  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+
+  % Collect all AOs and plists
+  [as, ao_invars, rest] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+  [pl, pl_invars, rest] = utils.helper.collect_objects(rest(:), 'plist', in_names);
+
+  % Decide on a deep copy or a modify
+  bs = copy(as, nargout);
+
+  % Combine plists
+  pl = combine(pl, plist(rest(:)),  getDefaultPlist); 
+
+  inhists = [];
+
+  %% Go through each input AO
+  for jj = 1 : numel(bs)
+    % gather the input history objects
+    inhists = [inhists bs(jj).hist]; %#ok<AGROW>
+
+    % check this is a time-series object
+    if ~isa(bs(jj).data, 'tsdata')
+      warning('!!! spsd requires tsdata (time-series) inputs. Skipping AO %s', ao_invars{jj}); %#ok<WNTAG>
+    else
+
+      % Check the time range.
+      time_range = find(pl, 'times');
+      if ~isempty(time_range)
+        bs(jj) = split(bs(jj), plist('method', 'times', 'times', time_range));
+      end
+      % Check the length of the object
+      if bs(jj).len <= 0
+        error('### The object is empty! Please revise your settings ...');
+      end
+
+      % pl = utils.helper.process_spectral_options(pl, 'log');
+      pl = pl.combine(getDefaultPlist());
+
+      % getting data
+      y = bs(jj).y;
+
+      % Window function
+      Win = find(pl, 'Win');
+      nfft = length(y);
+      Win = ao( combine(plist('win', Win , 'length', nfft), pl) );
+
+      % detrend
+      order = find(pl,'order');
+      if ~(order < 0)
+        y = ltpda_polyreg(y,  order).';
+      else 
+        y = reshape(y, 1, nfft);
+      end
+
+      % computing PSD
+      window = Win.data.y;
+      window = window/norm(window)*sqrt(nfft);
+      yASD = real(fft(y.*window, nfft)).^2 + imag(fft(y.*window, nfft)).^2;
+      pow = [yASD(1) yASD(2:floor(nfft/2))*2];
+      pow = pow /  ( bs(jj).data.fs * nfft);
+      Freqs = linspace(0, bs(jj).data.fs/2, nfft/2);
+
+      % smoothing PSD
+      if ~isempty(find(pl,'frequencies'))
+        error('the option "frequencies" is deprecated, frequencies are "removed" by default')
+      end
+      [Freqs, pow, nFreqs, nDofs] = ltpda_spsd(Freqs, pow, find(pl,'linCoef'), find(pl,'logCoef') );
+      % create new output fsdata
+      scale = find(pl, 'Scale');
+      switch lower(scale)
+        case 'asd'
+          fsd = fsdata(Freqs, sqrt(pow), bs(jj).data.fs);
+          fsd.setYunits(bs(jj).data.yunits / unit('Hz^0.5'));
+          %           stdDev = 0.5 * sqrt( pow ./ nDofs ); % linear approximation of the sqrt of a distribution
+          % approximation knowing the STD of the PSD
+          % STD assuming amplitude samples are independent, Chi^1_2 distibuted
+          % (with both variables of powe expectancy pow/2), and of different
+          % magnitude
+          stdDev = 2 * sqrt(pow./nDofs) .* ( nDofs - 2*exp( 2*(gammaln((nDofs+1)/2)-gammaln(nDofs/2)) ) ); % std of the chi_2N^1 
+        case 'psd'
+          fsd = fsdata(Freqs, pow, bs(jj).data.fs);
+          fsd.setYunits(bs(jj).data.yunits.^2/unit('Hz'));
+          % STD assuming power samples are independent, Chi^2_2 distibuted
+          % (with both variables of expectancy pow/2), and of different
+          % magnitude
+          stdDev = sqrt(2) * (pow./nDofs) .* sqrt(2*nDofs);  % std of the chi_2N^2 
+        otherwise
+          error(['### Unknown scaling:' scale]);
+      end
+
+      fsd.setXunits('Hz');
+      fsd.setDx(nFreqs*Freqs(2)/2);
+      fsd.setEnbw(1);% WARNING HERE!!!
+      fsd.setT0(bs(jj).data.t0);
+      % make output analysis object
+      bs(jj).data = fsd;
+      % set name
+      bs(jj).name = ['SPSD(', ao_invars{jj}, ') ' upper(scale)];
+      % Add standard deviation
+      bs(jj).data.dy = stdDev;
+      % Add history
+      bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), inhists(jj));
+
+    end % End tsdata if/else
+  end % End AO loop
+
+  %% Set output
+  if nargout == numel(bs)
+    % List of outputs
+    for ii = 1:numel(bs)
+      varargout{ii} = bs(ii);
+    end
+  else
+    % Single output
+    varargout{1} = bs;
+  end
+
+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, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: spsd.m,v 1.21 2011/07/11 10:43:35 adrien Exp $', sets, pl);
+end
+
+%--------------------------------------------------------------------------
+% Get Default Plist
+%--------------------------------------------------------------------------
+function pl = getDefaultPlist()
+
+  % Plist for Welch-based, log-scale spaced spectral estimators.
+  pl = plist;
+
+  % Win
+  p = param({'Win',['the window to be applied to the data to remove the ', ...
+    'discontinuities at edges of segments. [default: taken from user prefs] <br>', ...
+    'Only the design parameters of the window object are used. Enter either: <ul>', ...
+    '<li> a specwin window object OR</li>', ...
+    '<li> a string value containing the window name</li></ul>', ...
+    'e.g., <tt>plist(''Win'', ''Kaiser'', ''psll'', 200)</tt>']}, paramValue.WINDOW);
+  pl.append(p);
+
+  % Psll
+  p = param({'Psll',['the peak sidelobe level for Kaiser windows.<br>', ...
+    'Note: it is ignored for all other windows']}, paramValue.DOUBLE_VALUE(200));
+  pl.append(p);
+
+  % Psll
+  p = param({'levelOrder','the contracting order for levelledHanning window'}, paramValue.DOUBLE_VALUE(2));
+  pl.append(p);
+
+  % Order
+  p = param({'Order',['order of segment detrending:<ul>', ...
+    '<li>-1 - no detrending</li>', ...
+    '<li>0 - subtract mean</li>', ...
+    '<li>1 - subtract linear fit</li>', ...
+    '<li>N - subtract fit of polynomial, order N</li></ul>']}, paramValue.DETREND_ORDER);
+  p.val.setValIndex(2);
+  pl.append(p);
+
+  % Times
+  p = param({'Times','time range. If not empty, sets the restricted interval to analyze'}, paramValue.DOUBLE_VALUE([]));
+  pl.append(p);
+
+  % Scale
+  p = param({'Scale',['scaling of output. Choose from:<ul>', ...
+    '<li>PSD - Power Spectral Density</li>', ...
+    '<li>ASD - Amplitude (linear) Spectral Density</li>'...
+    ]}, {1, {'PSD', 'ASD', 'PS', 'AS'}, paramValue.SINGLE});
+  pl.append(p);
+
+  p = param( {'lincoef', 'Linear scale smoothing coefficent (freq. bins)'}, 1);
+  pl.append(p);
+
+  p = param( {'logcoef', ['Logarithmic scale smoothing coefficent<br>', 'Best compromise for both axes is 2/3']}, 2/3);
+  pl.append(p);
+end