0001 function varargout = ltpda_pwelch(varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 ALGONAME = mfilename;
0027 VERSION = '$Id: ltpda_pwelch.html,v 1.1 2007/06/08 14:15:11 hewitson Exp $';
0028
0029
0030 bs = [];
0031
0032
0033 if nargin == 1
0034 in = char(varargin{1});
0035 if strcmp(in, 'Params')
0036 varargout{1} = getDefaultPL();
0037 return
0038 end
0039 end
0040
0041
0042 invars = {};
0043 as = [];
0044 pls = [];
0045 for j=1:nargin
0046 if isa(varargin{j}, 'ao')
0047 as = [as varargin{j}];
0048 invars = [invars cellstr(inputname(j))];
0049 end
0050 if isa(varargin{j}, 'plist')
0051 pls = [pls varargin{j}];
0052 end
0053 end
0054 if isempty(pls)
0055 pls = plist();
0056 end
0057 na = length(as);
0058
0059
0060 for j=1:na
0061
0062
0063 a = as(j);
0064 d = a.data;
0065
0066 if isa(a, 'ao')
0067
0068 if isa(d, 'tsdata')
0069
0070
0071 plo = combine(pls, getDefaultPL());
0072
0073
0074 Nfft = find(plo, 'Nfft');
0075 setWindow = 0;
0076 if Nfft <= 0
0077 Nfft = length(d.x);
0078 disp(sprintf('! Using default Nfft of %g', Nfft))
0079 setWindow = 1;
0080 end
0081 plo = set(plo, 'Nfft', Nfft);
0082
0083 Win = find(plo, 'Win');
0084 if setWindow || length(Win.win)~=Nfft
0085 switch Win.name
0086 case {'Kaiser', 'Flattop'}
0087 Win = specwin(Win.name, Nfft, Win.psll);
0088 otherwise
0089 Win = specwin(Win.name, Nfft);
0090 end
0091 disp(sprintf('! Reset window to %s(%d)', strrep(Win.name, '_', '\_'), length(Win.win)))
0092 end
0093 plo = set(plo, 'Win', Win);
0094
0095 Nolap = find(plo, 'Nolap');
0096 if Nolap <= 0
0097 Nolap = round(Win.rov*Nfft/100);
0098 disp(sprintf('! Using default overlap of %d samples', Nolap))
0099 end
0100 plo = set(plo, 'Nolap', Nolap);
0101
0102
0103 if length(Win.win) ~= Nfft
0104 error('### Window length should be same as Nfft.');
0105 end
0106
0107
0108 [f, pxx, info] = ltpda_psd(d.x, Win.win, Nolap, Nfft, d.fs, 'ASD', d.yunits);
0109
0110
0111 fs = fsdata(f, pxx, d.fs);
0112 fs = set(fs, 'name', sprintf('PSD(%s)', d.name));
0113 fs = set(fs, 'yunits', info.units);
0114 fs = set(fs, 'enbw', Win.nenbw);
0115 fs = set(fs, 'navs', info.navs);
0116
0117
0118 h = history(ALGONAME, VERSION, plo, a.hist);
0119 h = set(h, 'invars', invars);
0120
0121
0122 b = ao(fs, h);
0123
0124
0125 if na~=length(invars) || isempty(invars{j})
0126 nj = a.name;
0127 else
0128 nj = invars{j};
0129 end
0130 b = set(b, 'name', sprintf('PSD(%s)', nj));
0131
0132
0133 bs = [bs b];
0134
0135 else
0136 warning('### Ignoring input AO number %d (%s); it is not a time-series.', j, a.name)
0137 end
0138 end
0139 end
0140
0141
0142 varargout{1} = bs;
0143
0144
0145
0146
0147 function plo = getDefaultPL()
0148
0149 disp('* creating default plist...');
0150 plo = plist();
0151 plo = append(plo, param('Nfft', -1));
0152 plo = append(plo, param('Win', specwin('Kaiser', 1, 200)));
0153 plo = append(plo, param('Nolap', -1));
0154 disp('* done.');
0155
0156
0157