Mercurial > hg > ltpda

changeset 43:bc767aaa99a8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
CVS Update
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Tue, 06 Dec 2011 11:09:25 +0100
parents f0afece42f48
children 409a22968d5e
files m-toolbox/classes/+utils/@math/math.m m-toolbox/classes/+utils/@math/rjsample.m m-toolbox/classes/+utils/@math/welchscale.m m-toolbox/classes/+utils/@models/built_in_model_template.m m-toolbox/classes/+utils/@models/mainFnc.m m-toolbox/classes/@ao/ao.m m-toolbox/classes/@ao/cohere.m m-toolbox/classes/@ao/cpsd.m m-toolbox/classes/@ao/tfe.m m-toolbox/classes/@ao/wosa.m m-toolbox/classes/@ao/xspec.m m-toolbox/classes/@ltpda_uo/fromModel.m m-toolbox/classes/@ltpda_uo/retrieve.m m-toolbox/classes/@ltpda_uoh/fromModel.m m-toolbox/classes/@ltpda_uoh/ltpda_uoh.m m-toolbox/classes/@matrix/fromInput.m m-toolbox/classes/@ssm/ssm.m
diffstat 17 files changed, 460 insertions(+), 52 deletions(-) [+]
line wrap: on
line diff
--- a/m-toolbox/classes/+utils/@math/math.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/+utils/@math/math.m	Tue Dec 06 11:09:25 2011 +0100
@@ -11,7 +11,7 @@
 % HISTORY: M Hewitson 26-05-08
 %              Creation
 %
-% VERSION: $Id: math.m,v 1.77 2011/10/07 08:19:06 miquel Exp $
+% VERSION: $Id: math.m,v 1.79 2011/12/01 09:41:22 hewitson Exp $
 %
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -27,7 +27,7 @@
     %-------------------------------------------------------------
     % List other methods
     %-------------------------------------------------------------
-    
+    varargout = welchscale(varargin);
     varargout = intfact(varargin); % Compute two integers P and Q
     varargout = cpf(varargin)
     varargout = lp2z(varargin)
@@ -80,7 +80,7 @@
     best = diffStepFish(i1,i2,S11,S12,S21,S22,N,meval,params,ngrid,ranges,freqs,inNames,outNames)
     best = diffStepFish_1x1(i1,S11,N,meval,params,values,ngrid,ranges,freqs,inNames,outNames)    
     loglk = loglikelihood(varargin)
-    loglk = loglikelihood_ssm(varargin)
+    [loglk snr] = loglikelihood_ssm(varargin)
     [loglk snr] = loglikelihood_matrix(varargin)
     snrexp = stnr(tmplt1,tmplt2,out1,out2,InvS11,InvS22,InvS12,InvS21)
     loglk = loglikelihood_ssm_td(xp,in,out,parnames,model,inNames,outNames,Noise,varargin)
@@ -101,7 +101,7 @@
     Covar = corr2cov(CorrC,SigC)
     R = Rcovmat(x)
     smpl = mhsample_td(model,in,out,cov,number,limit,parnames,Tc,xi,xo,search,jumps,parplot,dbg_info,inNames,outNames,inNoise,inNoiseNames,cutbefore,cutafter)
-    Bxy = rjsample(model,in,out,nse,cov,number,limit,param,Tc,xi,xo,search,jumps,parplot,dbg_info,inNames,outNames,inModel,outModel)
+    [Bxy LogLambda] = rjsample(mmdl,fin,fout,mnse,cvar,N,rang,param,Tc,xi,xo,search,jumps,parplot,debug,inNames,outNames,anneal,SNR0,DeltaL,inModel,outModel); 
     [Fout,x] = ecdf(y)
     cVal = SKcriticalvalues(n1,n2,alph)
     x = Finv(p,n1,n2)
--- a/m-toolbox/classes/+utils/@math/rjsample.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/+utils/@math/rjsample.m	Tue Dec 06 11:09:25 2011 +0100
@@ -287,6 +287,7 @@
     figure (2)
      for jj = 3:(nummodels+2)
        %if  (mnacc(nacc,jj-2) ~= 0)
+       smpl(smpl==0) = nan;
        plot(smpl(:,jj-2),'color',coloor(jj-2,:))
        legen = [legen ; sprintf('model%d',jj-2)]; % legend
        legend(legen)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/welchscale.m	Tue Dec 06 11:09:25 2011 +0100
@@ -0,0 +1,54 @@
+% WELCHSCALE scales the output of welch to be in the required units
+%
+% $Id: welchscale.m,v 1.1 2011/12/01 09:41:22 hewitson Exp $
+%
+
+function [yy, dyy, info] = welchscale(xx, dxx, win, fs, norm, inunits)
+  
+  nfft = length(win);
+  S1   = sum(win);
+  S2   = sum(win.^2);
+  enbw = fs * S2 / (S1*S1);
+  
+  if isempty(norm)
+    norm = 'None';
+  end
+  switch lower(norm)
+    case 'asd'
+      yy = sqrt(xx);
+      if isempty(dxx)
+        dyy = dxx;
+      else
+        dyy = 1./2./sqrt(xx) .* dxx;
+      end
+      info.units = inunits ./ unit('Hz^0.5');
+    case 'psd'
+      yy = xx;
+      dyy = dxx;
+      info.units = inunits.^2/unit('Hz');
+    case 'as'
+      yy = sqrt(xx * enbw);
+      if isempty(dxx)
+        dyy = dxx;
+      else
+        dyy = 1./2./sqrt(xx) .* dxx * enbw;
+      end
+      info.units = inunits;
+    case 'ps'
+      yy = xx * enbw;
+      dyy = dxx * enbw;
+      info.units = inunits.^2;
+    case 'none'
+      yy = xx;
+      dyy = dxx;
+      info.units = inunits;
+    otherwise
+      error('Unknown normalisation');
+  end
+  
+  info.nfft = nfft;
+  info.enbw = enbw;
+  info.norm = norm;
+  
+end
+
--- a/m-toolbox/classes/+utils/@models/built_in_model_template.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/+utils/@models/built_in_model_template.m	Tue Dec 06 11:09:25 2011 +0100
@@ -20,7 +20,7 @@
 % REFERENCES:
 %
 %
-% VERSION:     $Id: built_in_model_template.m,v 1.1 2011/04/29 07:05:54 hewitson Exp $
+% VERSION:     $Id: built_in_model_template.m,v 1.2 2011/12/02 09:15:11 hewitson Exp $
 %
 % HISTORY:
 %
@@ -99,11 +99,14 @@
     return;
   end
   
-  % build model
-  pl = varargin{1};
-
+  %---  Input plist
+  % Changes to this plist will be captured in the history
+  historyPlist = varargin{1}; 
+  % This can be locally modified without the changes being reflected in the history
+  localPlist   = copy(historyPlist, 1); 
+  
   % Get parameters
-  p = pl.find('My Parameter');
+  p = localPlist.find('My Parameter');
   
   % Build the model object the way you want
   
@@ -124,6 +127,6 @@
 %--------------------------------------------------------------------------
 function v = getVersion
   
-  v = '$Id: built_in_model_template.m,v 1.1 2011/04/29 07:05:54 hewitson Exp $';
+  v = '$Id: built_in_model_template.m,v 1.2 2011/12/02 09:15:11 hewitson Exp $';
   
 end
--- a/m-toolbox/classes/+utils/@models/mainFnc.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/+utils/@models/mainFnc.m	Tue Dec 06 11:09:25 2011 +0100
@@ -20,7 +20,7 @@
 %             getVersion - a function handle to the getVersion function
 %           versionTable - a function handle to the versionTable function
 % 
-% VERSION: $Id: mainFnc.m,v 1.3 2011/04/29 14:25:14 hewitson Exp $
+% VERSION: $Id: mainFnc.m,v 1.4 2011/12/02 09:03:20 hewitson Exp $
 % 
 % 
 function varargout = mainFnc(inputs, modelFilename, getModelDescription, getModelDocumentation, getVersion, versionTable)
@@ -39,13 +39,12 @@
   end
   
   % Build the object
-  hpl = copy(pl, 1);
   out = fcn(pl);
   
   % Set the method version string in the minfo object
   if ~isempty(constructorInfo) && utils.helper.isSubclassOf(class(out), 'ltpda_uoh')
     % If this is a user-call via a constructor, then we add history
-    out = addHistoryStep(out, constructorInfo, hpl);
+    out = addHistoryStep(out, constructorInfo, pl);
   end
   
   if nargout > 0
--- a/m-toolbox/classes/@ao/ao.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ao/ao.m	Tue Dec 06 11:09:25 2011 +0100
@@ -29,7 +29,7 @@
 %
 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'ao')">Parameters Description</a>
 %
-% VERSION: $Id: ao.m,v 1.361 2011/08/22 05:23:45 hewitson Exp $
+% VERSION: $Id: ao.m,v 1.362 2011/12/01 08:30:52 hewitson Exp $
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -110,7 +110,7 @@
 %                             given in the parameter list
 %
 %
-% VERSION:    $Id: ao.m,v 1.361 2011/08/22 05:23:45 hewitson Exp $
+% VERSION:    $Id: ao.m,v 1.362 2011/12/01 08:30:52 hewitson Exp $
 %
 % Parameter sets for plist constructor (in order of priority):
 %
@@ -775,7 +775,7 @@
     end
     
     function out = VEROUT()
-      out = '$Id: ao.m,v 1.361 2011/08/22 05:23:45 hewitson Exp $';
+      out = '$Id: ao.m,v 1.362 2011/12/01 08:30:52 hewitson Exp $';
     end
     
     function ii = getInfo(varargin)
@@ -1358,6 +1358,7 @@
   methods (Static = true, Access = private)
     % constructor functions
     % Spectral estimate function
+    varargout = wosa(varargin);
     [yy, dyy, info] = welchscale(xx, dxx, win, fs, norm, inunits)
     [x,M,isreal_x,y,Ly,win,winName,winParam,noverlap,k,L,options] = welchparse(x,esttype,varargin)
     [P,f] = computeperiodogram(x,win,nfft,esttype,varargin)
--- a/m-toolbox/classes/@ao/cohere.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ao/cohere.m	Tue Dec 06 11:09:25 2011 +0100
@@ -20,7 +20,7 @@
 %
 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'cohere')">Parameters Description</a>
 %
-% VERSION:    $Id: cohere.m,v 1.43 2011/04/08 08:56:16 hewitson Exp $
+% VERSION:    $Id: cohere.m,v 1.44 2011/12/01 09:35:42 hewitson Exp $
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -82,7 +82,7 @@
     pl   = getDefaultPlist();
   end
   % Build info object
-  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: cohere.m,v 1.43 2011/04/08 08:56:16 hewitson Exp $', sets, pl);
+  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: cohere.m,v 1.44 2011/12/01 09:35:42 hewitson Exp $', sets, pl);
   ii.setModifier(false);
   ii.setArgsmin(2);
 end
@@ -110,5 +110,8 @@
     '<li>C  - Complex Coherence:<br><tt>Sxy / sqrt(Sxx * Syy)</tt></li></ul>']}, {1, {'C', 'MS'}, paramValue.SINGLE});
   pl.append(p);
   
+  p = param('newMethod', paramValue.TRUE_FALSE);
+  pl.append(p);  
+  
 end
 
--- a/m-toolbox/classes/@ao/cpsd.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ao/cpsd.m	Tue Dec 06 11:09:25 2011 +0100
@@ -14,7 +14,7 @@
 %
 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'cpsd')">Parameters Description</a>
 %
-% VERSION:    $Id: cpsd.m,v 1.38 2011/04/08 08:56:13 hewitson Exp $
+% VERSION:    $Id: cpsd.m,v 1.39 2011/12/01 08:51:34 hewitson Exp $
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -68,7 +68,7 @@
     pl   = getDefaultPlist();
   end
   % Build info object
-  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: cpsd.m,v 1.38 2011/04/08 08:56:13 hewitson Exp $', sets, pl);
+  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: cpsd.m,v 1.39 2011/12/01 08:51:34 hewitson Exp $', sets, pl);
   ii.setModifier(false);
   ii.setArgsmin(2);
 end
@@ -90,5 +90,7 @@
   % General plist for Welch-based, linearly spaced spectral estimators
   pl = plist.WELCH_PLIST;
   
+  p = param('newMethod', paramValue.TRUE_FALSE);
+  pl.append(p);
 end
 
--- a/m-toolbox/classes/@ao/tfe.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ao/tfe.m	Tue Dec 06 11:09:25 2011 +0100
@@ -17,7 +17,7 @@
 %
 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'tfe')">Parameters Description</a>
 %
-% VERSION:    $Id: tfe.m,v 1.36 2011/04/08 08:56:13 hewitson Exp $
+% VERSION:    $Id: tfe.m,v 1.37 2011/12/01 09:01:48 hewitson Exp $
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -71,7 +71,7 @@
     pl   = getDefaultPlist();
   end
   % Build info object
-  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: tfe.m,v 1.36 2011/04/08 08:56:13 hewitson Exp $', sets, pl);
+  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: tfe.m,v 1.37 2011/12/01 09:01:48 hewitson Exp $', sets, pl);
   ii.setModifier(false);
   ii.setArgsmin(2);
 end
@@ -92,5 +92,7 @@
   % General plist for Welch-based, linearly spaced spectral estimators
   pl = plist.WELCH_PLIST;
   
+  p = param('newMethod', paramValue.TRUE_FALSE);
+  pl.append(p);
 end
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/wosa.m	Tue Dec 06 11:09:25 2011 +0100
@@ -0,0 +1,317 @@
+% WOSA implements Welch's overlaped segmented averaging algorithm with
+% segment detrending and variance estimation.
+% 
+% [pxx, f, info] = wosa(x,type,pl)
+% [pxx, f, info] = wosa(x,y,type,pl)
+%
+% INPUTS:      x    - input analysis objects
+%              y    - input analysis objects
+%              type - type of estimation:
+%                       'psd'      - compute Power Spectral Denstiy (PSD)
+%                       'cpsd'     - compute cross-spectral density
+%                       'tfe'      - estimate transfer function between inputs
+%                       'mscohere' - estimate magnitude-squared cross-coherence
+%                       'cohere'   - estimate complex cross-coherence
+%              pl   - input parameter list
+%
+% PARAMETERS: 'Win'   - a specwin window object [default: Kaiser -200dB psll]
+%             'Olap' - segment percent overlap [default: taken from window function]
+%             'Nfft'  - number of samples in each fft [default: length of input data]
+%             'Scale' - one of
+%                                'ASD' - amplitude spectral density
+%                                'PSD' - power spectral density [default]
+%                                'AS'  - amplitude spectrum
+%                                'PS'  - power spectrum
+%                       * applies only to spectrum 'Type' 'psd'
+%             'Order' - order of segment detrending
+%                        -1 - no detrending
+%                         0 - subtract mean [default]
+%                         1 - subtract linear fit
+%                         N - subtract fit of polynomial, order N
+%
+% Version: $Id: wosa.m,v 1.5 2011/12/02 07:08:11 hewitson Exp $
+%
+
+function varargout = wosa(varargin)
+  import utils.const.*
+  
+  % Process inputs
+  if nargin == 3
+    a  = varargin{1};
+    esttype = varargin{2};
+    pl = varargin{3};
+    inunits = a.data.yunits;
+    L = a.len;
+  else
+    a  = varargin{1};
+    b  = varargin{2};
+    esttype = varargin{3};
+    pl = varargin{4};
+    if a.data.fs ~= b.data.fs
+      error('The two time-series have different sample rates.');
+    end
+    inunits = b.data.yunits / a.data.yunits;
+    L = min(a.len, b.len);
+  end
+  
+  % Parse inputs
+  win          = find(pl, 'Win');
+  nfft         = find(pl, 'Nfft');
+  percentOlap  = find(pl, 'Olap')/100;
+  scale        = find(pl, 'scale');
+  xOlap        = round(percentOlap*nfft);
+  detrendOrder = find(pl, 'order');
+  fs           = a.fs;
+  winVals      = win.win.'; % because we always get a column from ao.y
+  
+  % Compute segment details
+  
+  nSegments = fix((L - xOlap)./(nfft - xOlap));
+  utils.helper.msg(msg.PROC3, 'N segment: %d', nfft);
+  utils.helper.msg(msg.PROC3, 'N overlap: %d', xOlap);
+  utils.helper.msg(msg.PROC3, 'N segments: %d', nSegments);
+  
+  % Compute start and end indices of each segment
+  segmentStep = nfft-xOlap;
+  segmentStarts = 1:segmentStep:nSegments*segmentStep;
+  segmentEnds   = segmentStarts+nfft-1;
+  
+  % Estimate the averaged periodogram for the desired quantity
+  switch lower(esttype)
+    case 'psd'
+      % Compute averaged periodogram
+      [Sxx, Svxx] = psdPeriodogram(a, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
+    case 'cpsd'
+      [Sxx, Svxx] = cpsdPeriodogram(a, b, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
+    case 'tfe'
+      [Sxx, Sxy, Syy] = tfePeriodogram(a, b, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
+    case {'mscohere','cohere'}
+      [Sxx, Sxy, Syy] = tfePeriodogram(a, b, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
+    otherwise
+      error('Unknown estimation type %s', esttype);
+  end
+  
+  % Scale to PSD
+  switch lower(esttype)
+    case {'psd','cpsd'}
+      [P, Pvxx] = scaleToPSD(Sxx, Svxx, nfft, fs);
+      % the 1/nSegments factor should come after welchscale if we don't
+      % want to apply sqrt() to it.
+      % We correct for that here. It is only needed for 'asd','as' in
+      % psd/cpsd, the other cases go always through 'PSD'.
+      if (strcmpi(scale,'PSD') || strcmpi(scale,'PS'))
+        dP = Pvxx;
+      elseif (strcmpi(scale,'ASD') || strcmpi(scale,'AS'))
+        dP = Pvxx/nSegments;
+      else
+        error('### Unknown scale')
+      end
+    case 'tfe'
+      % Compute the 1-sided or 2-sided PSD [Power/freq] or mean-square [Power].
+      % Also, corresponding freq vector and freq units.
+      % In the Cross PSD, the frequency vector and xunits are not used.
+      Pxx = scaleToPSD(Sxx, [], nfft, fs);
+      Pxy = scaleToPSD(Sxy, [], nfft, fs);
+      Pyy = scaleToPSD(Syy, [], nfft, fs);
+      % mean and std
+      P = Pxy ./ Pxx; % Txy
+      if nSegments == 1
+        dP =[];
+      else
+        dP = (nSegments/(nSegments-1)^2)*(Pyy./Pxx).*(1 - (abs(Pxy).^2)./(Pxx.*Pyy));
+      end
+    case 'mscohere'
+      % Magnitude Square Coherence estimate.
+      % Auto PSD for 2nd input vector. The freq vector & xunits are not
+      % used.
+      Pxx = scaleToPSD(Sxx, [], nfft, fs);
+      Pxy = scaleToPSD(Sxy, [], nfft, fs);
+      Pyy = scaleToPSD(Syy, [], nfft, fs);
+      % mean and std
+      P = (abs(Pxy).^2)./(Pxx.*Pyy); % Magnitude-squared coherence
+      dP = (2*P/nSegments).*(1-P).^2;
+    case 'cohere'
+      % Complex Coherence estimate.
+      % Auto PSD for 2nd input vector. The freq vector & xunits are not
+      % used.
+      Pxx = scaleToPSD(Sxx, [], nfft, fs);
+      Pxy = scaleToPSD(Sxy, [], nfft, fs);
+      Pyy = scaleToPSD(Syy, [], nfft, fs);
+      P = Pxy./sqrt(Pxx.*Pyy); % Complex coherence
+      dP = (2*abs(P)/nSegments).*(1-abs(P)).^2;
+  
+  end
+  
+  % Compute frequencies
+  freqs = psdfreqvec('npts', nfft,'Fs', fs, 'Range', 'half').';
+  
+  % Scale to required units
+  [Pxx, dP, info] = utils.math.welchscale(P, dP, winVals, fs, scale, inunits);
+  info.navs = nSegments;
+  
+  if nSegments ==1
+    dev = [];
+  else
+    dev = sqrt(dP);
+  end
+  
+  % Set outputs
+  varargout = {Pxx, freqs, info, dev};
+    
+end
+
+% scale averaged periodogram to PSD
+function [Pxx, Pvxx] = scaleToPSD(Sxx, Svxx, nfft, fs)
+  
+  % Take 1-sided spectrum which means we double the power in the
+  % appropriate bins
+  if rem(nfft,2),
+    indices = 1:(nfft+1)/2;  % ODD
+    Sxx1sided = Sxx(indices,:);
+    % double the power except for the DC bin
+    Sxx = [Sxx1sided(1,:); 2*Sxx1sided(2:end,:)];  
+    if ~isempty(Svxx)
+      Svxx1sided = Svxx(indices,:);
+      Svxx = [Svxx1sided(1,:); 4*Svxx1sided(2:end,:)];
+    end
+  else
+    indices = 1:nfft/2+1;    % EVEN
+    Sxx1sided = Sxx(indices,:);
+    % Double power except the DC bin and the Nyquist bin
+    Sxx = [Sxx1sided(1,:); 2*Sxx1sided(2:end-1,:); Sxx1sided(end,:)];
+    if ~isempty(Svxx)
+      Svxx1sided = Svxx(indices,:); % Take only [0,pi] or [0,pi)
+      Svxx = [Svxx1sided(1,:); 4*Svxx1sided(2:end-1,:); Svxx1sided(end,:)];
+    end
+  end
+
+  % Now scale to PSD
+  Pxx   = Sxx./fs;
+  Pvxx  = Svxx./fs^2;
+  
+end
+
+% compute tfe
+function [Sxx, Sxy, Syy] = tfePeriodogram(x, y, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder)
+  
+  nfft = segmentEnds(1);
+  Sxx = zeros(nfft,1); % Initialize Sxx
+  Sxy = zeros(nfft,1); % Initialize Sxy
+  Syy = zeros(nfft,1); % Initialize Syy
+  % loop over segments
+  for ii = 1:nSegments
+    if detrendOrder < 0
+      xseg = x.y(segmentStarts(ii):segmentEnds(ii));
+      yseg = y.y(segmentStarts(ii):segmentEnds(ii));
+    else
+      [xseg,coeffs] = ltpda_polyreg(x.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
+      [yseg,coeffs] = ltpda_polyreg(y.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
+    end
+
+    % Compute periodograms
+    Sxxk = wosa_periodogram(xseg, [], winVals, nfft);
+    Sxyk = wosa_periodogram(yseg, xseg, winVals, nfft);
+    Syyk = wosa_periodogram(yseg, [], winVals, nfft);
+      
+    Sxx = Sxx + Sxxk;
+    Sxy = Sxy + Sxyk;
+    Syy = Syy + Syyk;
+    % don't need to be divided by k because only rations are used here
+  end
+  
+end
+
+% compute cpsd
+function [Sxx, Svxx] = cpsdPeriodogram(x, y, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder)
+  
+  Mnxx = 0; 
+  Mn2xx = 0;
+  nfft = segmentEnds(1);
+  for ii = 1:nSegments
+    if detrendOrder < 0
+      xseg = x.y(segmentStarts(ii):segmentEnds(ii));
+      yseg = y.y(segmentStarts(ii):segmentEnds(ii));
+    else
+      [xseg,coeffs] = ltpda_polyreg(x.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
+      [yseg,coeffs] = ltpda_polyreg(y.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
+    end
+    
+    % Compute periodogram
+    Sxxk = wosa_periodogram(xseg, yseg, winVals, nfft);
+    
+    % Welford's algorithm to update mean and variance
+    Qxx = Sxxk - Mnxx;
+    Mnxx = Mnxx +Qxx/ii;
+    Mn2xx = Mn2xx + abs(Qxx.*conj(Sxxk - Mnxx));
+  end
+  Sxx = Mnxx;
+  if nSegments ==1
+    Svxx = [];
+  else
+    Svxx = Mn2xx/(nSegments-1)/nSegments;
+  end
+  
+  
+end
+
+% compute psd
+function [Sxx, Svxx] = psdPeriodogram(x, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder)
+  
+  Mnxx = 0; 
+  Mn2xx = 0;
+  nfft = segmentEnds(1);
+  % Loop over the segments
+  for ii = 1:nSegments
+    % Detrend if desired
+    if detrendOrder < 0
+      seg = x.y(segmentStarts(ii):segmentEnds(ii));
+    else
+      [seg,coeffs] = ltpda_polyreg(x.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
+    end
+    % Compute periodogram
+    Sxxk = wosa_periodogram(seg, [], winVals,nfft);
+    % Welford's algorithm for updating mean and variance
+    if ii == 1
+      Mnxx = Sxxk;
+    else
+      Qxx = Sxxk - Mnxx;
+      Mnxx = Mnxx + Qxx/ii;
+      Mn2xx = Mn2xx + Qxx.*(Sxxk - Mnxx);
+    end
+  end
+  Sxx = Mnxx;
+  if nSegments == 1
+    Svxx = [];
+  else
+    Svxx = Mn2xx/(nSegments-1)/nSegments;
+  end
+  
+end
+
+% Scaled periodogram of one or two input signals
+function Sxx = wosa_periodogram(x, y, win, nfft)
+  
+  % window data
+  xwin = x.*win;
+  isCross = false;
+  if ~isempty(y)
+    ywin = y.*win;
+    isCross = true;
+  end
+  
+  % take fft
+  X = fft(xwin, nfft);
+  if isCross
+    Y = fft(ywin, nfft);
+  end
+  
+  % Compute scale factor to compensate for the window power
+  K = win'*win;
+  
+  % Compute scaled power
+  Sxx = X.*conj(X)/K;
+  if isCross,
+    Sxx = X.*conj(Y)/K;
+  end  
+  
+end
--- a/m-toolbox/classes/@ao/xspec.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ao/xspec.m	Tue Dec 06 11:09:25 2011 +0100
@@ -21,7 +21,7 @@
 %
 % OUTPUTS:     b  - output AO
 %
-% VERSION:     $Id: xspec.m,v 1.61 2011/05/22 22:09:44 mauro Exp $
+% VERSION:     $Id: xspec.m,v 1.62 2011/12/01 08:51:11 hewitson Exp $
 %
 % HISTORY:     30-05-2007 M Hewitson
 %                 Creation
@@ -39,7 +39,7 @@
   mi     = varargin{4};
   invars = varargin{5};
 
-  VERSION  = '$Id: xspec.m,v 1.61 2011/05/22 22:09:44 mauro Exp $';
+  VERSION  = '$Id: xspec.m,v 1.62 2011/12/01 08:51:11 hewitson Exp $';
 
   % Set the method version string in the minfo object
   mi.setMversion([VERSION '-->' mi.mversion]);
@@ -133,8 +133,12 @@
   % -------- Make Xspec estimate
 
   % Compute xspec using welch and always scale to PSD
-  [txy, f, info, dev] = welch(tsao(1), tsao(2), method, usepl.pset('Scale', 'PSD'));
-
+  if usepl.find('newmethod')
+    [txy, f, info, dev] = ao.wosa(tsao(1), tsao(2), method, usepl.pset('Scale', 'PSD'));
+  else
+    [txy, f, info, dev] = welch(tsao(1), tsao(2), method, usepl.pset('Scale', 'PSD'));
+  end
+  
   % Keep the data shape of the first input AO
   if size(tsao(1).data.y,1) == 1
     txy = txy.';
--- a/m-toolbox/classes/@ltpda_uo/fromModel.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ltpda_uo/fromModel.m	Tue Dec 06 11:09:25 2011 +0100
@@ -12,7 +12,7 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-function [obj, ii, fcnname] = fromModel(obj, pl)
+function [obj, ii, fcnname, pl] = fromModel(obj, pl)
   
   % Get the name of the model the user wants
   model = find(pl, 'built-in');
--- a/m-toolbox/classes/@ltpda_uo/retrieve.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ltpda_uo/retrieve.m	Tue Dec 06 11:09:25 2011 +0100
@@ -28,7 +28,7 @@
 % If only a single object is requested, it is returned as an object,
 % not packed in a cell array.
 %
-% VERSION:     $Id: retrieve.m,v 1.28 2011/07/01 14:38:57 ingo Exp $
+% VERSION:     $Id: retrieve.m,v 1.29 2011/11/30 06:49:01 hewitson Exp $
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -141,6 +141,7 @@
         % Retrieve the bytes
         q =  sprintf('select mat from bobjs where obj_id="%d"', ids(j));
         results = conn.query(q);
+%         dd = '';
         while results.next
           dd = results.getObject(1);
         end
@@ -205,7 +206,7 @@
     
   catch ME
     fprintf(2, [ME.message, '\n\n']);
-    utils.helper.msg(msg.PROC1, '### Retrieve error.');
+    conn.setLocked(false);
     rethrow(ME)
   end
   
@@ -246,7 +247,7 @@
     pl   = getDefaultPlist;
   end
   % Build info object
-  ii = minfo(mfilename, 'ltpda_uo', 'ltpda', utils.const.categories.internal, '$Id: retrieve.m,v 1.28 2011/07/01 14:38:57 ingo Exp $', sets, pl);
+  ii = minfo(mfilename, 'ltpda_uo', 'ltpda', utils.const.categories.internal, '$Id: retrieve.m,v 1.29 2011/11/30 06:49:01 hewitson Exp $', sets, pl);
   ii.setModifier(false);
 end
 
--- a/m-toolbox/classes/@ltpda_uoh/fromModel.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ltpda_uoh/fromModel.m	Tue Dec 06 11:09:25 2011 +0100
@@ -10,13 +10,13 @@
 %
 % PARAMETER:   pl: Parameter list object
 %
-% VERSION:     $Id: fromModel.m,v 1.13 2011/08/15 12:59:04 hewitson Exp $
+% VERSION:     $Id: fromModel.m,v 1.14 2011/12/02 09:05:15 hewitson Exp $
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 function obj = fromModel(obj, pl)
 
-  [obj, ii, fcnname] = fromModel@ltpda_uo(obj, pl);
+  [obj, ii, fcnname, pl] = fromModel@ltpda_uo(obj, pl);
   % Add history
   try 
     feval(fcnname, 'info');
@@ -26,6 +26,20 @@
     obj.addHistory(ii, pl, [], []);
   end
   
+  % Here we make a check on the default plist keys because these properties
+  % may have been set in the model. We make the assumption that if they are
+  % empty, the user had no intention of setting them so we can remove them
+  % and keep the value set in the model.
+  if isempty(pl.find('name'))
+    pl.remove('name');
+  end
+  if isempty(pl.find('description'))
+    pl.remove('description');
+  end
+  if isempty(pl.find('plotinfo'))
+    pl.remove('plotinfo');
+  end
+  
   % Set object properties
   obj.setObjectProperties(pl);
   
--- a/m-toolbox/classes/@ltpda_uoh/ltpda_uoh.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ltpda_uoh/ltpda_uoh.m	Tue Dec 06 11:09:25 2011 +0100
@@ -20,7 +20,7 @@
 %       plotinfo      - plist with the plot information
 %       procinfo      - plist with additional information for an object.
 %
-% VERSION:     $Id: ltpda_uoh.m,v 1.64 2011/09/16 04:59:23 hewitson Exp $
+% VERSION:     $Id: ltpda_uoh.m,v 1.65 2011/12/02 09:03:42 hewitson Exp $
 %
 % SEE ALSO:    ltpda_obj, ltpda_uo, ltpda_tf, ltpda_filter, 
 %              ao, miir, mfir, filterbank, timespan, pzmodel, history, ssm, 
@@ -157,7 +157,7 @@
     end
     
     function out = VEROUT()
-      out = '$Id: ltpda_uoh.m,v 1.64 2011/09/16 04:59:23 hewitson Exp $';
+      out = '$Id: ltpda_uoh.m,v 1.65 2011/12/02 09:03:42 hewitson Exp $';
     end
     
     function out = SETS()
@@ -202,6 +202,7 @@
     function pl = removeGlobalKeys(pl)
       pl.remove('name');
       pl.remove('description');
+      pl.remove('plotinfo');
     end
   end
   
--- a/m-toolbox/classes/@matrix/fromInput.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@matrix/fromInput.m	Tue Dec 06 11:09:25 2011 +0100
@@ -7,20 +7,34 @@
 %
 % CALL:        matrix = matrix.fromInput(inobjs)
 %
-% VERSION:     $Id: fromInput.m,v 1.9 2011/02/14 13:20:24 ingo Exp $
+% VERSION:     $Id: fromInput.m,v 1.10 2011/12/05 07:40:23 hewitson Exp $
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function obj = fromInput(obj, pli, callerIsMethod)
   
+  VERSION = '$Id: fromInput.m,v 1.10 2011/12/05 07:40:23 hewitson Exp $';
+  
   import utils.const.*
   
   if callerIsMethod
     % do nothing
-    pl = pli;
   else
-    utils.helper.msg(msg.OPROC1, 'Construct from inputs');
-    pl = combine(pli, matrix.getDefaultPlist('From Input'));
+    % get AO info
+    ii = matrix.getInfo('matrix', 'From Input');
+    
+    % Set the method version string in the minfo object
+    ii.setMversion([VERSION '-->' ii.mversion]);
   end
+  
+  % Combine input plist with default values
+  if callerIsMethod
+    pl        = pli;
+  else
+    % Combine input plist with default values
+    % TODO: the parse step should be removed and included somehow into plist/applyDefaults
+    pl = applyDefaults(ii.plists, pli);
+  end  
+  
   shape = pl.find('shape');
   objs  = pl.find('objs');
   
@@ -53,17 +67,9 @@
       inhists = [inobjs(:).hist];
     end
     obj.addHistory(matrix.getInfo('matrix', 'None'), pl, [], inhists);
-    warning('off', utils.const.warnings.METHOD_NOT_FOUND);
-    % remove parameters we already used
-    pl_set = copy(pl,1);
-    if pl_set.isparam('objs')
-      pl_set.remove('objs');
-    end
-    if pl_set.isparam('shape')
-      pl_set.remove('shape');
-    end
-    obj.setProperties(pl_set);
-    warning('on', utils.const.warnings.METHOD_NOT_FOUND);
+
+    % Set any remaining object properties which exist in the default plist
+    obj.setObjectProperties(pl);
   end
   
 end
--- a/m-toolbox/classes/@ssm/ssm.m	Wed Nov 23 19:22:13 2011 +0100
+++ b/m-toolbox/classes/@ssm/ssm.m	Tue Dec 06 11:09:25 2011 +0100
@@ -107,7 +107,7 @@
 %
 % <a href="matlab:utils.helper.displayMethodInfo('ssm', 'ssm')">Parameters Description</a>
 %
-% VERSION: $Id: ssm.m,v 1.199 2011/05/13 15:14:39 ingo Exp $
+% VERSION: $Id: ssm.m,v 1.200 2011/12/02 09:02:35 hewitson Exp $
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -716,7 +716,7 @@
     end
     
     function out = VEROUT()
-      out = '$Id: ssm.m,v 1.199 2011/05/13 15:14:39 ingo Exp $';
+      out = '$Id: ssm.m,v 1.200 2011/12/02 09:02:35 hewitson Exp $';
     end
     
     function out = SETS()
@@ -753,7 +753,7 @@
         case 'from built-in model'
           % Built-in
           % This is inherited
-          pl = plist.FROM_BUILT_IN;
+          pl = combine(pl, plist.FROM_BUILT_IN);
           % withParams --> withparams changed to 'symbolic params'
           p = param({'symbolic params',['Give a cell-array of parameter names to keep in the expression.<br>',...
             'By default this is empty and the model will be returned fully numeric.',...