LPSD implement the LPSD algorithm for analysis objects. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION: LPSD implement the LPSD algorithm for analysis objects. CALL: bs = lpsd(as) INPUTS: as - array of analysis objects pl - parameter list (see below) OUTPUTS: bs - array of analysis objects, one for each input PARAMETER LIST: Kdes - desired number of averages (default 100) Lmin - minimum segment length (default 0) Jdes - number of spectral frequencies to compute (default fs/2) win - a specwin window object Only the design parameters of the window object are used; the window is recomputed for each DFT length inside the lpsd_core algorithm. Olap - desired overlap percentage (default: taken from window) Order - order of detrending -1 - no detrending 0 - subtract mean 1 - subtract linear fit N - subtract fit of polynomial, order N Scale - Scaling of output. Choose from: AS - Amplitude (linear) Spectrum ASD - Amplitude (linear) Spectral Density PS - Power Spectrum PSD - Power Spectral Density [default] M-FILE INFO: Get information about this methods by calling >> ao.getInfo('lpsd') Get information about a specified set-plist by calling: >> ao.getInfo('lpsd', 'None') VERSION: $Id: lpsd.m,v 1.13 2008/09/05 14:14:58 hewitson Exp $ HISTORY: 02-02-07 M Hewitson Created References: "Improved spectrum estimation from digitized time series on a logarithmic frequency axis", Michael Troebs, Gerhard Heinzel, Measurement 39 (2006) 120-129. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 % LPSD implement the LPSD algorithm for analysis objects. 0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0003 % 0004 % DESCRIPTION: LPSD implement the LPSD algorithm for analysis objects. 0005 % 0006 % CALL: bs = lpsd(as) 0007 % 0008 % INPUTS: as - array of analysis objects 0009 % pl - parameter list (see below) 0010 % 0011 % OUTPUTS: bs - array of analysis objects, one for each input 0012 % 0013 % PARAMETER LIST: 0014 % 0015 % Kdes - desired number of averages (default 100) 0016 % Lmin - minimum segment length (default 0) 0017 % Jdes - number of spectral frequencies to compute (default fs/2) 0018 % win - a specwin window object 0019 % Only the design parameters of the window object are used; the 0020 % window is recomputed for each DFT length inside the lpsd_core 0021 % algorithm. 0022 % Olap - desired overlap percentage (default: taken from window) 0023 % Order - order of detrending 0024 % -1 - no detrending 0025 % 0 - subtract mean 0026 % 1 - subtract linear fit 0027 % N - subtract fit of polynomial, order N 0028 % Scale - Scaling of output. Choose from: 0029 % AS - Amplitude (linear) Spectrum 0030 % ASD - Amplitude (linear) Spectral Density 0031 % PS - Power Spectrum 0032 % PSD - Power Spectral Density [default] 0033 % 0034 % M-FILE INFO: Get information about this methods by calling 0035 % >> ao.getInfo('lpsd') 0036 % 0037 % Get information about a specified set-plist by calling: 0038 % >> ao.getInfo('lpsd', 'None') 0039 % 0040 % VERSION: $Id: lpsd.m,v 1.13 2008/09/05 14:14:58 hewitson Exp $ 0041 % 0042 % HISTORY: 02-02-07 M Hewitson 0043 % Created 0044 % 0045 % References: "Improved spectrum estimation from digitized time series 0046 % on a logarithmic frequency axis", Michael Troebs, Gerhard Heinzel, 0047 % Measurement 39 (2006) 120-129. 0048 % 0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0050 0051 function varargout = lpsd(varargin) 0052 0053 import utils.const.* 0054 0055 % Check if this is a call for parameters 0056 if utils.helper.isinfocall(varargin{:}) 0057 varargout{1} = getInfo(varargin{3}); 0058 return 0059 end 0060 0061 utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename); 0062 0063 % Collect input variable names 0064 in_names = cell(size(varargin)); 0065 for ii = 1:nargin,in_names{ii} = inputname(ii);end 0066 0067 % Collect all AOs and plists 0068 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); 0069 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names); 0070 0071 % Decide on a deep copy or a modify 0072 bs = copy(as, nargout); 0073 0074 % Combine plists 0075 pl = combine(pl, getDefaultPlist); 0076 0077 % Go through each input AO 0078 for j=1:numel(bs) 0079 % check this is a time-series object 0080 if ~isa(bs(j).data, 'tsdata') 0081 warning('!!! lpsd requires tsdata (time-series) inputs. Skipping AO %s', ao_invars{j}); 0082 else 0083 % Desired number of averages 0084 Kdes = find(pl, 'Kdes'); 0085 % num desired spectral frequencies 0086 Jdes = find(pl, 'Jdes'); 0087 % Window function 0088 Win = find(pl, 'Win'); 0089 % Overlap 0090 Nolap = find(pl, 'Olap'); 0091 if Nolap < 0 0092 % use window function rov 0093 Nolap = Win.rov/100; 0094 utils.helper.msg(msg.PROC1, 'using recommended overlap of %d', Nolap); 0095 pl.pset('Olap', Nolap); 0096 end 0097 % Order of detrending 0098 Order = find(pl, 'Order'); 0099 % Minimum segment length 0100 Lmin = find(pl, 'Lmin'); 0101 % Get frequency vector 0102 [f, r, m, L, K] = ao.ltf_plan(length(bs(j).data.y), bs(j).data.fs, Nolap, 1, Lmin, Jdes, Kdes); 0103 % compute LPSD 0104 try 0105 if find(pl, 'M-FILE ONLY') 0106 % Using pure m-file version 0107 [P, Pxx, ENBW] = ao.mlpsd_m(bs(j).data.y, f, r, m, L, bs(j).data.fs, Win, Order, Nolap); 0108 else 0109 [P, Pxx, ENBW] = ao.mlpsd_mex(bs(j).data.y, f, r, m, L, bs(j).data.fs, Win, Order, Win.rov, Lmin); 0110 end 0111 catch 0112 warning('!!! mex file dft failed. Using m-file version of lpsd.'); 0113 % Using pure m-file version 0114 [P, Pxx, ENBW] = ao.mlpsd_m(bs(j).data.y, f, r, m, L, bs(j).data.fs, Win, Order, Nolap); 0115 end 0116 0117 % Keep the data shape of the input AO 0118 if size(bs(j).data.y,1) == 1 0119 P = P.'; 0120 Pxx = Pxx.'; 0121 f = f.'; 0122 end 0123 0124 % create new output fsdata 0125 scale = find(pl, 'Scale'); 0126 switch scale 0127 case 'AS' 0128 fsd = fsdata(f, sqrt(P), bs(j).data.fs); 0129 fsd.setYunits(bs(j).data.yunits); 0130 case 'ASD' 0131 fsd = fsdata(f, sqrt(Pxx), bs(j).data.fs); 0132 fsd.setYunits(bs(j).data.yunits / unit('Hz^0.5')); 0133 case 'PS' 0134 fsd = fsdata(f, P, bs(j).data.fs); 0135 fsd.setYunits(bs(j).data.yunits.^2); 0136 case 'PSD' 0137 fsd = fsdata(f, Pxx, bs(j).data.fs); 0138 fsd.setYunits(bs(j).data.yunits.^2/unit('Hz')); 0139 otherwise 0140 error(['### Unknown scaling:' scale]); 0141 end 0142 fsd.setXunits('Hz'); 0143 fsd.setEnbw(ENBW); 0144 % make output analysis object 0145 bs(j).data = fsd; 0146 % set name 0147 bs(j).setName(sprintf('l%s(%s)', scale, ao_invars{j}), 'internal'); 0148 % Add history 0149 bs(j).addHistory(getInfo, pl, ao_invars(j), bs(j).hist); 0150 % Add processing info 0151 bs(j).procinfo = plist('r', r, 'm', m, 'l', L, 'k', K); 0152 0153 end % End tsdata if/else 0154 end % End AO loop 0155 0156 % Set output 0157 varargout{1} = bs; 0158 0159 end 0160 0161 %-------------------------------------------------------------------------- 0162 % Get Info Object 0163 %-------------------------------------------------------------------------- 0164 function ii = getInfo(varargin) 0165 if nargin == 1 && strcmpi(varargin{1}, 'None') 0166 sets = {}; 0167 pl = []; 0168 else 0169 sets = {'Default'}; 0170 pl = getDefaultPlist; 0171 end 0172 % Build info object 0173 ii = minfo(mfilename, 'ao', '', 'Signal Processing', '$Id: lpsd.m,v 1.13 2008/09/05 14:14:58 hewitson Exp $', sets, pl); 0174 end 0175 0176 %-------------------------------------------------------------------------- 0177 % Get Default Plist 0178 %-------------------------------------------------------------------------- 0179 function pl = getDefaultPlist() 0180 pl = plist('Kdes', 100, 'Jdes', 1000, 'Lmin', 0,... 0181 'Win', getappdata(0, 'ltpda_default_spectral_window'), ... 0182 'Olap', -1, 'Order', 0, 'Scale', 'PSD'); 0183 end 0184