LTPDA_PWELCH makes power spectral density estimates of the time-series objects %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION: LTPDA_PWELCH makes power spectral density estimates of the time-series objects in the input analysis objects. CALL: b = ltpda_pwelch(a1,a2,a3,...,pl) INPUTS: b - output analysis objects aN - input analysis objects pl - input parameter list Makes PSD estimates using ltpda_psd() of each input AO. 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: fs of input data] '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 VERSION: $Id: ltpda_pwelch.m,v 1.18 2008/02/20 07:39:01 hewitson Exp $ HISTORY: 07-02-2007 M Hewitson Creation The following call returns a parameter list object that contains the default parameter values: >> pl = ltpda_pwelch('Params') The following call returns a string that contains the routine CVS version: >> version = ltpda_pwelch(ao,'Version') The following call returns a string that contains the routine category: >> category = ltpda_pwelch(ao,'Category') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 function varargout = ltpda_pwelch(varargin) 0002 % LTPDA_PWELCH makes power spectral density estimates of the time-series objects 0003 % 0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0005 % 0006 % DESCRIPTION: LTPDA_PWELCH makes power spectral density estimates of the 0007 % time-series objects in the input analysis objects. 0008 % 0009 % CALL: b = ltpda_pwelch(a1,a2,a3,...,pl) 0010 % 0011 % INPUTS: b - output analysis objects 0012 % aN - input analysis objects 0013 % pl - input parameter list 0014 % 0015 % Makes PSD estimates using ltpda_psd() of each input AO. 0016 % 0017 % If the last input argument is a parameter list (plist) it is used. 0018 % The following parameters are recognised. 0019 % 0020 % PARAMETERS: 'Win' - a specwin window object [default: Kaiser -200dB psll] 0021 % 'olap' - segment percent overlap [default: taken from window function] 0022 % 'Nfft' - number of samples in each fft [default: fs of input data] 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 % VERSION: $Id: ltpda_pwelch.m,v 1.18 2008/02/20 07:39:01 hewitson Exp $ 0035 % 0036 % HISTORY: 07-02-2007 M Hewitson 0037 % Creation 0038 % 0039 % The following call returns a parameter list object that contains the 0040 % default parameter values: 0041 % 0042 % >> pl = ltpda_pwelch('Params') 0043 % 0044 % The following call returns a string that contains the routine CVS version: 0045 % 0046 % >> version = ltpda_pwelch(ao,'Version') 0047 % 0048 % The following call returns a string that contains the routine category: 0049 % 0050 % >> category = ltpda_pwelch(ao,'Category') 0051 % 0052 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0053 0054 ALGONAME = mfilename; 0055 VERSION = '$Id: ltpda_pwelch.m,v 1.18 2008/02/20 07:39:01 hewitson Exp $'; 0056 CATEGORY = 'Signal Processing'; 0057 0058 % Initialise output 0059 bs = []; 0060 0061 % Check if this is a call for parameters 0062 if nargin == 1 0063 in = char(varargin{1}); 0064 if strcmpi(in, 'Params') 0065 varargout{1} = getDefaultPL(); 0066 return 0067 elseif strcmpi(in, 'Version') 0068 varargout{1} = VERSION; 0069 return 0070 elseif strcmpi(in, 'Category') 0071 varargout{1} = CATEGORY; 0072 return 0073 end 0074 end 0075 0076 % Collect input ao's, plist's and ao variable names 0077 in_names = {}; 0078 for ii = 1:nargin 0079 in_names{end+1} = inputname(ii); 0080 end 0081 [as, pls, invars] = collect_inputs(varargin, in_names); 0082 na = length(as); 0083 0084 % Loop over input AOs 0085 for j=1:na 0086 0087 % Unpack AO 0088 a = as(j); 0089 d = a.data; 0090 0091 if isa(a, 'ao') 0092 % check input data 0093 if isa(d, 'tsdata') 0094 0095 % combine plists 0096 plo = combine(pls, getDefaultPL()); 0097 0098 % check all defaults are sensible 0099 Nfft = find(plo, 'Nfft'); % Number of points in each fft 0100 setWindow = 0; 0101 if Nfft <= 0 0102 Nfft = length(d.y); 0103 disp(sprintf('! Using default Nfft of %g', Nfft)) 0104 setWindow = 1; 0105 end 0106 plo = pset(plo, 'Nfft', Nfft); 0107 0108 Win = find(plo, 'Win'); % Window object 0109 if setWindow || length(Win.win)~=Nfft 0110 switch Win.name 0111 case {'Kaiser', 'Flattop'} 0112 Win = specwin(Win.name, Nfft, Win.psll); 0113 otherwise 0114 Win = specwin(Win.name, Nfft); 0115 end 0116 disp(sprintf('! Reset window to %s(%d)', strrep(Win.name, '_', '\_'), length(Win.win))) 0117 end 0118 plo = pset(plo, 'Win', Win); 0119 0120 olap = find(plo, 'olap'); % Amount to overlap each fft 0121 if olap <= 0 0122 olap = Win.rov; 0123 disp(sprintf('! Using default overlap of %d /%', olap)) 0124 end 0125 plo = pset(plo, 'olap', olap); 0126 0127 % Get detrending order 0128 order = find(plo, 'Order'); 0129 0130 % scaling of spectrum 0131 scale = find(plo, 'Scale'); 0132 0133 % check parameters 0134 if length(Win.win) ~= Nfft 0135 error('### Window length should be same as Nfft.'); 0136 end 0137 0138 % Compute PSD using pwelch 0139 [f, pxx, info] = ltpda_psd(d.y, Win.win, round(olap/100*Nfft), Nfft, d.fs, scale, d.yunits, order); 0140 0141 % create new output fsdata 0142 fs = fsdata(f, pxx, d.fs); 0143 fs = set(fs, 'name', sprintf('%s(%s)', scale, d.name)); 0144 fs = set(fs, 'yunits', info.units); 0145 fs = set(fs, 'enbw', Win.nenbw*d.fs/Nfft); 0146 fs = set(fs, 'navs', info.navs); 0147 0148 % create new output history 0149 h = history(ALGONAME, VERSION, plo, a.hist); 0150 h = set(h, 'invars', invars); 0151 0152 % make output analysis object 0153 b = ao(fs, h); 0154 0155 % set name 0156 if na~=length(invars) || isempty(invars{j}) 0157 nj = a.name; 0158 else 0159 nj = invars{j}; 0160 end 0161 b = setnh(b, 'name', sprintf('%s(%s)', scale, nj)); 0162 0163 % add to outputs 0164 bs = [bs b]; 0165 0166 else 0167 warning('### Ignoring input AO number %d (%s); it is not a time-series.', j, a.name) 0168 end 0169 end 0170 end % loop over analysis objects 0171 0172 % set output 0173 varargout{1} = bs; 0174 0175 0176 %-------------------------------------------------------------------------- 0177 % Get default params 0178 function plo = getDefaultPL() 0179 0180 plo = plist(); 0181 plo = append(plo, param('Nfft', -1)); 0182 plo = append(plo, param('Win', specwin('Kaiser', 1, 200))); 0183 plo = append(plo, param('Olap', -1)); 0184 plo = append(plo, param('Scale', 'PSD')); 0185 plo = append(plo, param('Order', 0)); 0186 0187 0188