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. CALL: b = psd(a1,a2,a3,...,pl) INPUTS: b - output analysis objects aN - input analysis objects pl - input parameter list Makes PSD estimates of each input AO using the Welch Overlap method. If the last input argument is a parameter list (plist) it is used. The following parameters are recognised. 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] A string value containing the variable 'fs' can also be used, e.g., plist('Nfft', '2*fs') 'Scale' - one of 'ASD' - amplitude spectral density 'PSD' - power spectral density [default] 'AS' - amplitude spectrum 'PS' - power spectrum 'Order' - order of segment detrending -1 - no detrending 0 - subtract mean [default] 1 - subtract linear fit N - subtract fit of polynomial, order N M-FILE INFO: Get information about this methods by calling >> ao.getInfo('psd') Get information about a specified set-plist by calling: >> ao.getInfo('psd', 'None') VERSION: $Id: psd.m,v 1.6 2008/09/05 11:05:29 ingo Exp $ HISTORY: 07-02-2007 M Hewitson Creation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 % PSD makes power spectral density estimates of the time-series objects 0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0003 % 0004 % DESCRIPTION: PSD makes power spectral density estimates of the 0005 % time-series objects in the input analysis objects. 0006 % 0007 % CALL: b = psd(a1,a2,a3,...,pl) 0008 % 0009 % INPUTS: b - output analysis objects 0010 % aN - input analysis objects 0011 % pl - input parameter list 0012 % 0013 % Makes PSD estimates of each input AO using the Welch Overlap method. 0014 % 0015 % If the last input argument is a parameter list (plist) it is used. 0016 % The following parameters are recognised. 0017 % 0018 % PARAMETERS: 'Win' - a specwin window object [default: Kaiser -200dB psll] 0019 % 'Olap' - segment percent overlap [default: taken from window function] 0020 % 'Nfft' - number of samples in each fft [default: length of input data] 0021 % A string value containing the variable 'fs' can 0022 % also be used, e.g., plist('Nfft', '2*fs') 0023 % 'Scale' - one of 0024 % 'ASD' - amplitude spectral density 0025 % 'PSD' - power spectral density [default] 0026 % 'AS' - amplitude spectrum 0027 % 'PS' - power spectrum 0028 % 'Order' - order of segment detrending 0029 % -1 - no detrending 0030 % 0 - subtract mean [default] 0031 % 1 - subtract linear fit 0032 % N - subtract fit of polynomial, order N 0033 % 0034 % M-FILE INFO: Get information about this methods by calling 0035 % >> ao.getInfo('psd') 0036 % 0037 % Get information about a specified set-plist by calling: 0038 % >> ao.getInfo('psd', 'None') 0039 % 0040 % VERSION: $Id: psd.m,v 1.6 2008/09/05 11:05:29 ingo Exp $ 0041 % 0042 % HISTORY: 07-02-2007 M Hewitson 0043 % Creation 0044 % 0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0046 0047 function varargout = psd(varargin) 0048 0049 % Check if this is a call for parameters 0050 if utils.helper.isinfocall(varargin{:}) 0051 varargout{1} = getInfo(varargin{3}); 0052 return 0053 end 0054 0055 import utils.const.* 0056 utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename); 0057 0058 % Collect input variable names 0059 in_names = cell(size(varargin)); 0060 for ii = 1:nargin,in_names{ii} = inputname(ii);end 0061 0062 % Collect all AOs and plists 0063 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); 0064 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names); 0065 0066 % Decide on a deep copy or a modify 0067 bs = copy(as, nargout); 0068 0069 % combine plists 0070 pl = combine(pl, getDefaultPlist()); 0071 0072 % Loop over input AOs 0073 for j=1:numel(bs) 0074 % check input data 0075 if isa(bs(j).data, 'tsdata') 0076 fs = bs(j).data.fs; 0077 %--------- Check all defaults are sensible 0078 % Check the number of points in FFT. If this is not set (<0) we set it 0079 % to be the length of the input data. 0080 Nfft = find(pl, 'Nfft'); 0081 setWindow = 0; 0082 if ischar(Nfft) 0083 nNfft = floor(eval(Nfft)); 0084 utils.helper.msg(msg.PROC1, 'setting Nfft to %s = %d', Nfft, nNfft); 0085 Nfft = nNfft; 0086 else 0087 if Nfft <= 0 0088 Nfft = length(bs(j).data.y); 0089 utils.helper.msg(msg.PROC1, 'using default Nfft of %g', Nfft); 0090 setWindow = 1; 0091 end 0092 end 0093 pl.pset('Nfft', Nfft); 0094 0095 % Check the window function. If the length of the window doesn't match 0096 % NFFT then we resize it. 0097 Win = find(pl, 'Win'); % Window object 0098 if setWindow || length(Win.win)~=Nfft 0099 switch Win.type 0100 case {'Kaiser', 'Flattop'} 0101 Win = specwin(Win.type, Nfft, Win.psll); 0102 otherwise 0103 Win = specwin(Win.type, Nfft); 0104 end 0105 utils.helper.msg(msg.PROC1, 'reset window to %s(%d)', strrep(Win.type, '_', '\_'), length(Win.win)); 0106 end 0107 pl.pset('Win', Win); 0108 0109 % Check the overlap. If this is not set, we take the overlap from that 0110 % recommended by the window function. 0111 Olap = find(pl, 'Olap'); 0112 if Olap <= 0 0113 Olap = Win.rov; 0114 utils.helper.msg(msg.PROC1, 'using default overlap of %2.1f%%', Olap); 0115 end 0116 pl.pset('Olap', Olap); 0117 0118 % scaling of spectrum 0119 scale = find(pl, 'Scale'); 0120 0121 % check parameters 0122 if length(Win.win) ~= Nfft 0123 error('### Window length should be same as Nfft.'); 0124 end 0125 % Compute PSD using pwelch 0126 [pxx, f, info] = ao.welch(bs(j), 'psd', pl); 0127 % Keep the data shape of the input AO 0128 if size(bs(j).data.y,1) == 1 0129 pxx = pxx.'; 0130 f = f.'; 0131 end 0132 % create new output fsdata 0133 fsd = fsdata(f, pxx, fs); 0134 fsd.setXunits('Hz'); 0135 fsd.setYunits(info.units); 0136 fsd.setEnbw(Win.nenbw*fs/Nfft); 0137 fsd.setNavs(info.navs); 0138 fsd.setT0(bs(j).data.t0); 0139 % update AO 0140 bs(j).data = fsd; 0141 % set name 0142 bs(j).setName(sprintf('%s(%s)', scale, ao_invars{j}), 'internal'); 0143 % Add history 0144 bs(j).addHistory(getInfo, pl, ao_invars(j), bs(j).hist); 0145 else 0146 warning('### Ignoring input AO number %d (%s); it is not a time-series.', j, bs(j).name) 0147 end 0148 end % loop over analysis objects 0149 0150 % set output 0151 varargout{1} = bs; 0152 end 0153 0154 %-------------------------------------------------------------------------- 0155 % Get Info Object 0156 %-------------------------------------------------------------------------- 0157 function ii = getInfo(varargin) 0158 if nargin == 1 && strcmpi(varargin{1}, 'None') 0159 sets = {}; 0160 pl = []; 0161 else 0162 sets = {'Default'}; 0163 pl = getDefaultPlist; 0164 end 0165 % Build info object 0166 ii = minfo(mfilename, 'ao', '', utils.const.categories.sigproc, '$Id: psd.m,v 1.6 2008/09/05 11:05:29 ingo Exp $', sets, pl); 0167 end 0168 0169 %-------------------------------------------------------------------------- 0170 % Get Default Plist 0171 %-------------------------------------------------------------------------- 0172 function pl = getDefaultPlist() 0173 pl = plist(... 0174 'Nfft', -1, ... 0175 'Win', getappdata(0, 'ltpda_default_spectral_window'), ... 0176 'Olap', -1, ... 0177 'Scale', 'PSD', ... 0178 'Order', 0); 0179 end 0180 % END 0181