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: length 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.html,v 1.14 2008/03/31 10:27:39 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('Version') The following call returns a string that contains the routine category: >> category = ltpda_pwelch('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: length 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.html,v 1.14 2008/03/31 10:27:39 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('Version') 0047 % 0048 % The following call returns a string that contains the routine category: 0049 % 0050 % >> category = ltpda_pwelch('Category') 0051 % 0052 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0053 0054 ALGONAME = mfilename; 0055 VERSION = '$Id: ltpda_pwelch.html,v 1.14 2008/03/31 10:27:39 hewitson Exp $'; 0056 CATEGORY = 'Signal Processing'; 0057 0058 % Initialise output 0059 bs = []; 0060 0061 %% Check if this is a call for parameters, the CVS version string 0062 % or the function category 0063 if nargin == 1 && ischar(varargin{1}) 0064 in = char(varargin{1}); 0065 if strcmp(in, 'Params') 0066 varargout{1} = getDefaultPL(); 0067 return 0068 elseif strcmp(in, 'Version') 0069 varargout{1} = VERSION; 0070 return 0071 elseif strcmp(in, 'Category') 0072 varargout{1} = CATEGORY; 0073 return 0074 end 0075 end 0076 0077 % Collect input ao's, plist's and ao variable names 0078 in_names = {}; 0079 for ii = 1:nargin 0080 in_names{end+1} = inputname(ii); 0081 end 0082 [as, pls, invars] = collect_inputs(varargin, in_names); 0083 na = length(as); 0084 0085 % Loop over input AOs 0086 for j=1:na 0087 0088 % Unpack AO 0089 a = as(j); 0090 d = a.data; 0091 0092 if isa(a, 'ao') 0093 % check input data 0094 if isa(d, 'tsdata') 0095 0096 % combine plists 0097 plo = combine(pls, getDefaultPL()); 0098 0099 % check all defaults are sensible 0100 Nfft = find(plo, 'Nfft'); % Number of points in each fft 0101 setWindow = 0; 0102 if Nfft <= 0 0103 Nfft = length(d.y); 0104 disp(sprintf('! Using default Nfft of %g', Nfft)) 0105 setWindow = 1; 0106 end 0107 plo = pset(plo, 'Nfft', Nfft); 0108 0109 Win = find(plo, 'Win'); % Window object 0110 if setWindow || length(Win.win)~=Nfft 0111 switch Win.name 0112 case {'Kaiser', 'Flattop'} 0113 Win = specwin(Win.name, Nfft, Win.psll); 0114 otherwise 0115 Win = specwin(Win.name, Nfft); 0116 end 0117 disp(sprintf('! Reset window to %s(%d)', strrep(Win.name, '_', '\_'), length(Win.win))) 0118 end 0119 plo = pset(plo, 'Win', Win); 0120 0121 Olap = find(plo, 'Olap'); % Amount to overlap each fft 0122 if Olap <= 0 0123 Olap = Win.rov; 0124 disp(sprintf('! Using default overlap of %2.1f%%', Olap)) 0125 end 0126 plo = pset(plo, 'Olap', Olap); 0127 0128 % Get detrending order 0129 order = find(plo, 'Order'); 0130 0131 % scaling of spectrum 0132 scale = find(plo, 'Scale'); 0133 0134 % check parameters 0135 if length(Win.win) ~= Nfft 0136 error('### Window length should be same as Nfft.'); 0137 end 0138 0139 % Compute PSD using pwelch 0140 [f, pxx, info] = ltpda_psd(d.y, Win.win, round(Olap/100*Nfft), Nfft, d.fs, scale, d.yunits, order); 0141 0142 % create new output fsdata 0143 fs = fsdata(f, pxx, d.fs); 0144 fs = set(fs, 'name', sprintf('%s(%s)', scale, d.name)); 0145 fs = set(fs, 'yunits', info.units); 0146 fs = set(fs, 'enbw', Win.nenbw*d.fs/Nfft); 0147 fs = set(fs, 'navs', info.navs); 0148 0149 % create new output history 0150 h = history(ALGONAME, VERSION, plo, a.hist); 0151 h = set(h, 'invars', invars); 0152 0153 % make output analysis object 0154 b = ao(fs, h); 0155 0156 % set name 0157 if na~=length(invars) || isempty(invars{j}) 0158 nj = a.name; 0159 else 0160 nj = invars{j}; 0161 end 0162 b = setnh(b, 'name', sprintf('%s(%s)', scale, nj)); 0163 0164 % add to outputs 0165 bs = [bs b]; 0166 0167 else 0168 warning('### Ignoring input AO number %d (%s); it is not a time-series.', j, a.name) 0169 end 0170 end 0171 end % loop over analysis objects 0172 0173 % set output 0174 varargout{1} = bs; 0175 0176 0177 %-------------------------------------------------------------------------- 0178 % Get default params 0179 function plo = getDefaultPL() 0180 0181 plo = plist('Nfft', -1, ... 0182 'Win', specwin('Kaiser', 1, 200), ... 0183 'Olap', -1, ... 0184 'Scale', 'PSD',... 0185 'Order', 0); 0186 0187 0188