line source
+ − % PSD makes power spectral density estimates of the time-series objects
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION: PSD makes power spectral density estimates of the
+ − % time-series objects in the input analysis objects
+ − % using the Welch Overlap method. PSD is computed
+ − % using a modified version of MATLAB's welch (>> help welch).
+ − %
+ − % CALL: bs = psd(a1,a2,a3,...,pl)
+ − % bs = psd(as,pl)
+ − % bs = as.psd(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', 'psd')">Parameters Description</a>
+ − %
+ − % VERSION: $Id: psd.m,v 1.69 2011/08/24 07:29:02 hewitson Exp $
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ −
+ − function varargout = psd(varargin)
+ −
+ − callerIsMethod = utils.helper.callerIsMethod;
+ −
+ − % 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
+ − usepl = applyDefaults(getDefaultPlist, varargin{:});
+ −
+ − inhists = [];
+ −
+ − % Loop over input AOs
+ − for jj = 1 : numel(bs)
+ − % gather the input history objects
+ − inhists = [inhists bs(jj).hist];
+ −
+ − % check input data
+ − if isa(bs(jj).data, 'tsdata')
+ −
+ − % Check the time range.
+ − time_range = mfind(usepl, '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
+ −
+ − % Utility to deal with nfft, win, olap etc
+ − use_pl = utils.helper.process_spectral_options(usepl, 'lin', len(bs(jj)), bs(jj).fs);
+ −
+ − % Compute PSD using pwelch
+ − [pxx, f, info, dev] = welch(bs(jj), 'psd', use_pl);
+ −
+ − % Keep the data shape of the input AO
+ − if size(bs(jj).data.y,1) == 1
+ − pxx = pxx.';
+ − f = f.';
+ − end
+ − % create new output fsdata
+ − fs = bs(jj).fs;
+ − fsd = fsdata(f, pxx, fs);
+ − fsd.setXunits('Hz');
+ − fsd.setYunits(info.units);
+ − fsd.setEnbw(info.enbw);
+ − fsd.setNavs(info.navs);
+ − fsd.setT0(bs(jj).data.t0);
+ − % update AO
+ − bs(jj).data = fsd;
+ − % add std deviation
+ − bs(jj).data.dy = dev;
+ − % set name: scaling of spectrum
+ − scale = upper(find(use_pl, 'Scale'));
+ − bs(jj).name = sprintf('%s(%s)', lower(scale), ao_invars{jj});
+ − % Add history
+ − if ~callerIsMethod
+ − bs(jj).addHistory(getInfo('None'), use_pl, ao_invars(jj), inhists(jj));
+ − end
+ − else
+ − warning('### Ignoring input AO number %d (%s); it is not a time-series.', jj, bs(jj).name)
+ − end
+ − 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: psd.m,v 1.69 2011/08/24 07:29:02 hewitson 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, linearly spaced spectral estimators
+ − pl = plist.WELCH_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
+ −