


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]
'Nolap' - segment overlap (number of samples) [default: taken from window function]
'Nfft' - number of samples in each fft [default: fs of input data]
'Debug' - debug level for terminal output (0-1)
VERSION: $Id: ltpda_pwelch.m,v 1.15 2007/11/08 10:35:04 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')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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 % 'Nolap' - segment overlap (number of samples) [default: taken from window function] 0022 % 'Nfft' - number of samples in each fft [default: fs of input data] 0023 % 'Debug' - debug level for terminal output (0-1) 0024 % 0025 % VERSION: $Id: ltpda_pwelch.m,v 1.15 2007/11/08 10:35:04 hewitson Exp $ 0026 % 0027 % HISTORY: 07-02-2007 M Hewitson 0028 % Creation 0029 % 0030 % The following call returns a parameter list object that contains the 0031 % default parameter values: 0032 % 0033 % >> pl = ltpda_pwelch('Params') 0034 % 0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0036 0037 ALGONAME = mfilename; 0038 VERSION = '$Id: ltpda_pwelch.m,v 1.15 2007/11/08 10:35:04 hewitson Exp $'; 0039 0040 % Initialise output 0041 bs = []; 0042 0043 % Check if this is a call for parameters 0044 if nargin == 1 0045 in = char(varargin{1}); 0046 if strcmp(in, 'Params') 0047 varargout{1} = getDefaultPL(); 0048 return 0049 end 0050 end 0051 0052 % Read inputs 0053 invars = {}; 0054 as = []; 0055 pls = []; 0056 for j=1:nargin 0057 if isa(varargin{j}, 'ao') 0058 as = [as varargin{j}]; 0059 invars = [invars cellstr(inputname(j))]; 0060 end 0061 if isa(varargin{j}, 'plist') 0062 pls = [pls varargin{j}]; 0063 end 0064 end 0065 if isempty(pls) 0066 pls = plist(); 0067 end 0068 na = length(as); 0069 0070 % Loop over input AOs 0071 for j=1:na 0072 0073 % Unpack AO 0074 a = as(j); 0075 d = a.data; 0076 0077 if isa(a, 'ao') 0078 % check input data 0079 if isa(d, 'tsdata') 0080 0081 % combine plists 0082 plo = combine(pls, getDefaultPL()); 0083 0084 % check all defaults are sensible 0085 Nfft = find(plo, 'Nfft'); % Number of points in each fft 0086 setWindow = 0; 0087 if Nfft <= 0 0088 Nfft = length(d.y); 0089 disp(sprintf('! Using default Nfft of %g', Nfft)) 0090 setWindow = 1; 0091 end 0092 plo = pset(plo, 'Nfft', Nfft); 0093 0094 Win = find(plo, 'Win'); % Window object 0095 if setWindow || length(Win.win)~=Nfft 0096 switch Win.name 0097 case {'Kaiser', 'Flattop'} 0098 Win = specwin(Win.name, Nfft, Win.psll); 0099 otherwise 0100 Win = specwin(Win.name, Nfft); 0101 end 0102 disp(sprintf('! Reset window to %s(%d)', strrep(Win.name, '_', '\_'), length(Win.win))) 0103 end 0104 plo = pset(plo, 'Win', Win); 0105 0106 Nolap = find(plo, 'Nolap'); % Amount to overlap each fft 0107 if Nolap <= 0 0108 Nolap = round(Win.rov*Nfft/100); 0109 disp(sprintf('! Using default overlap of %d samples', Nolap)) 0110 end 0111 plo = pset(plo, 'Nolap', Nolap); 0112 0113 % check parameters 0114 if length(Win.win) ~= Nfft 0115 error('### Window length should be same as Nfft.'); 0116 end 0117 0118 % Compute PSD using pwelch 0119 [f, pxx, info] = ltpda_psd(d.y, Win.win, Nolap, Nfft, d.fs, 'ASD', d.yunits); 0120 0121 % create new output fsdata 0122 fs = fsdata(f, pxx, d.fs); 0123 fs = set(fs, 'name', sprintf('PSD(%s)', d.name)); 0124 fs = set(fs, 'yunits', info.units); 0125 fs = set(fs, 'enbw', Win.nenbw*d.fs/Nfft); 0126 fs = set(fs, 'navs', info.navs); 0127 0128 % create new output history 0129 h = history(ALGONAME, VERSION, plo, a.hist); 0130 h = set(h, 'invars', invars); 0131 0132 % make output analysis object 0133 b = ao(fs, h); 0134 0135 % set name 0136 if na~=length(invars) || isempty(invars{j}) 0137 nj = a.name; 0138 else 0139 nj = invars{j}; 0140 end 0141 b = setnh(b, 'name', sprintf('PSD(%s)', nj)); 0142 0143 % add to outputs 0144 bs = [bs b]; 0145 0146 else 0147 warning('### Ignoring input AO number %d (%s); it is not a time-series.', j, a.name) 0148 end 0149 end 0150 end % loop over analysis objects 0151 0152 % set output 0153 varargout{1} = bs; 0154 0155 0156 %-------------------------------------------------------------------------- 0157 % Get default params 0158 function plo = getDefaultPL() 0159 0160 disp('* creating default plist...'); 0161 plo = plist(); 0162 plo = append(plo, param('Nfft', -1)); 0163 plo = append(plo, param('Win', specwin('Kaiser', 1, 200))); 0164 plo = append(plo, param('Nolap', -1)); 0165 disp('* done.'); 0166 0167 0168