# HG changeset patch # User Daniele Nicolodi # Date 1323166165 -3600 # Node ID bc767aaa99a84f5de0b1cf8fc2a143d97489e862 # Parent f0afece42f48a49931a4cd3d4c3a15b498020579 CVS Update diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/+utils/@math/math.m --- 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) diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/+utils/@math/rjsample.m --- 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) diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/+utils/@math/welchscale.m --- /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 + diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/+utils/@models/built_in_model_template.m --- 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 diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/+utils/@models/mainFnc.m --- 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 diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ao/ao.m --- 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 @@ % % Parameters Description % -% 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) diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ao/cohere.m --- 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 @@ % % Parameters Description % -% 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 @@ '
  • C - Complex Coherence:
    Sxy / sqrt(Sxx * Syy)
    • ']}, {1, {'C', 'MS'}, paramValue.SINGLE}); pl.append(p); + p = param('newMethod', paramValue.TRUE_FALSE); + pl.append(p); + end diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ao/cpsd.m --- 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 @@ % % Parameters Description % -% 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 diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ao/tfe.m --- 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 @@ % % Parameters Description % -% 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 diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ao/wosa.m --- /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 diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ao/xspec.m --- 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.'; diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ltpda_uo/fromModel.m --- 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'); diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ltpda_uo/retrieve.m --- 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 diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ltpda_uoh/fromModel.m --- 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); diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ltpda_uoh/ltpda_uoh.m --- 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 diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@matrix/fromInput.m --- 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 diff -r f0afece42f48 -r bc767aaa99a8 m-toolbox/classes/@ssm/ssm.m --- 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 @@ % % Parameters Description % -% 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.
    ',... 'By default this is empty and the model will be returned fully numeric.',...