line source
+ − % LPSD implements the LPSD algorithm for analysis objects.
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION: LPSD implements the LPSD algorithm for analysis objects.
+ − %
+ − % CALL: bs = lpsd(a1,a2,a3,...,pl)
+ − % bs = lpsd(as,pl)
+ − % bs = as.lpsd(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', 'lpsd')">Parameters Description</a>
+ − %
+ − % VERSION: $Id: lpsd.m,v 1.55 2011/05/22 21:22:09 mauro Exp $
+ − %
+ − % References: "Improved spectrum estimation from digitized time series
+ − % on a logarithmic frequency axis", Michael Troebs, Gerhard Heinzel,
+ − % Measurement 39 (2006) 120-129.
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function varargout = lpsd(varargin)
+ −
+ − % Check if this is a call for parameters
+ − if utils.helper.isinfocall(varargin{:})
+ − varargout{1} = getInfo(varargin{3});
+ − return
+ − end
+ −
+ − import utils.const.*
+ − 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
+ − [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+ −
+ − % Decide on a deep copy or a modify
+ − bs = copy(as, nargout);
+ −
+ − % Apply defaults to plist
+ − pl = applyDefaults(getDefaultPlist, varargin{:});
+ −
+ − inhists = [];
+ −
+ − % Loop over input AOs
+ − for jj = 1 : numel(bs)
+ − % gather the input history objects
+ − inhists = [inhists bs(jj).hist];
+ −
+ − % check this is a time-series object
+ − if ~isa(bs(jj).data, 'tsdata')
+ − warning('!!! lpsd requires tsdata (time-series) inputs. Skipping AO %s', ao_invars{jj});
+ − else
+ −
+ − % Check the time range.
+ − time_range = mfind(pl, 'split', 'times');
+ − if ~isempty(time_range)
+ − switch class(time_range)
+ − case 'double'
+ − bs(jj) = split(bs(jj), plist(...
+ − 'times', time_range));
+ − case 'timespan'
+ − bs(jj) = split(bs(jj), plist(...
+ − 'timespan', time_range));
+ − case 'time'
+ − bs(jj) = split(bs(jj), plist(...
+ − 'start_time', time_range(1), ...
+ − 'end_time', time_range(2)));
+ − case 'cell'
+ − bs(jj) = split(bs(jj), plist(...
+ − 'start_time', time_range{1}, ...
+ − 'end_time', time_range{2}));
+ − otherwise
+ − end
+ − 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');
+ −
+ − % Desired number of averages
+ − Kdes = find(pl, 'Kdes');
+ − % num desired spectral frequencies
+ − Jdes = find(pl, 'Jdes');
+ − % Minimum segment length
+ − Lmin = find(pl, 'Lmin');
+ − % Window function
+ − Win = find(pl, 'Win');
+ − % Overlap
+ − Nolap = find(pl, 'Olap')/100;
+ − % Order of detrending
+ − Order = find(pl, 'Order');
+ −
+ − % Get frequency vector
+ − [f, r, m, L, K] = ao.ltf_plan(length(bs(jj).data.y), bs(jj).data.fs, Nolap, 1, Lmin, Jdes, Kdes);
+ −
+ − % compute LPSD
+ − try
+ − if find(pl, 'M-FILE ONLY')
+ − % Using pure m-file version
+ − [P, Pxx, ENBW] = ao.mlpsd_m(bs(jj).data.y, f, r, m, L, bs(jj).data.fs, Win, Order, Nolap);
+ − else
+ − [P, Pxx, dev, devxx, ENBW] = ao.mlpsd_mex(bs(jj).data.y, f, r, m, L, bs(jj).data.fs, Win, Order, Nolap*100, Lmin);
+ − end
+ − catch ME
+ − warning('!!! mex file dft failed. Using m-file version of lpsd.');
+ − % Using pure m-file version
+ − [P, Pxx, ENBW] = ao.mlpsd_m(bs(jj).data.y, f, r, m, L, bs(jj).data.fs, Win, Order, Nolap);
+ − end
+ −
+ − % Keep the data shape of the input AO
+ − if size(bs(jj).data.y,1) == 1
+ − P = P.';
+ − Pxx = Pxx.';
+ − dev = dev.';
+ − devxx = devxx.';
+ − f = f.';
+ − end
+ −
+ − % create new output fsdata
+ − scale = find(pl, 'Scale');
+ − switch lower(scale)
+ − case 'as'
+ − fsd = fsdata(f, sqrt(P), bs(jj).data.fs);
+ − fsd.setYunits(bs(jj).data.yunits);
+ − std = sqrt(dev);
+ − case 'asd'
+ − fsd = fsdata(f, sqrt(Pxx), bs(jj).data.fs);
+ − fsd.setYunits(bs(jj).data.yunits / unit('Hz^0.5'));
+ − std = sqrt(devxx);
+ − case 'ps'
+ − fsd = fsdata(f, P, bs(jj).data.fs);
+ − fsd.setYunits(bs(jj).data.yunits.^2);
+ − std = dev;
+ − case 'psd'
+ − fsd = fsdata(f, Pxx, bs(jj).data.fs);
+ − fsd.setYunits(bs(jj).data.yunits.^2/unit('Hz'));
+ − std = devxx;
+ − otherwise
+ − error(['### Unknown scaling:' scale]);
+ − end
+ − fsd.setXunits('Hz');
+ − fsd.setEnbw(ENBW);
+ − fsd.setT0(bs(jj).data.t0);
+ − % make output analysis object
+ − bs(jj).data = fsd;
+ − % set name
+ − bs(jj).name = sprintf('L%s(%s)', upper(scale), ao_invars{jj});
+ − % Add processing info
+ − bs(jj).procinfo = plist('r', r, 'm', m, 'l', L, 'k', K);
+ − % Add standard deviation
+ − bs(jj).data.dy = std;
+ − % Add history
+ − bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), inhists(jj));
+ −
+ − end % End tsdata if/else
+ − end % loop over analysis objects
+ −
+ − % Set output
+ − varargout = utils.helper.setoutputs(nargout, bs);
+ − 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: lpsd.m,v 1.55 2011/05/22 21:22:09 mauro Exp $', sets, pl);
+ − end
+ −
+ − %--------------------------------------------------------------------------
+ − % Get Default Plist
+ − %--------------------------------------------------------------------------
+ − function plout = getDefaultPlist()
+ − persistent pl;
+ − if ~exist('pl', 'var') || isempty(pl)
+ − pl = buildplist();
+ − end
+ − plout = pl;
+ − end
+ −
+ − function pl = buildplist()
+ −
+ − % General plist for Welch-based, log-scale spaced spectral estimators
+ − pl = plist.LPSD_PLIST;
+ −
+ − % Scale
+ − p = param({'Scale',['The scaling of output. Choose from:<ul>', ...
+ − '<li>PSD - Power Spectral Density</li>', ...
+ − '<li>ASD - Amplitude (linear) Spectral Density</li>', ...
+ − '<li>PS - Power Spectrum</li>', ...
+ − '<li>AS - Amplitude (linear) Spectrum</li></ul>']}, {1, {'PSD', 'ASD', 'PS', 'AS'}, paramValue.SINGLE});
+ − pl.append(p);
+ −
+ − end
+ −
+ − % PARAMETERS:
+ − %
+ − % 'Kdes' - desired number of averages to perform [default: 100]
+ − % 'Jdes' - number of spectral frequencies to compute [default: 1000]
+ − % 'Lmin' - minimum segment length [default: 0]
+ − % 'Win' - the window to be applied to the data to remove the
+ − % discontinuities at edges of segments. [default: taken from
+ − % user prefs]
+ − % Only the design parameters of the window object are
+ − % used. Enter either:
+ − % - a specwin window object OR
+ − % - a string value containing the window name
+ − % e.g., plist('Win', 'Kaiser', 'psll', 200)
+ − % 'Olap' - segment percent overlap [default: -1, (taken from window function)]
+ − % 'Scale' - scaling of output. Choose from:
+ − % PSD - Power Spectral Density [default]
+ − % ASD - Amplitude (linear) Spectral Density
+ − % PS - Power Spectrum
+ − % AS - Amplitude (linear) Spectrum
+ − % 'Order' - order of segment detrending
+ − % -1 - no detrending
+ − % 0 - subtract mean [default]
+ − % 1 - subtract linear fit
+ − % N - subtract fit of polynomial, order N