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
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037 ALGONAME = mfilename;
0038 VERSION = '$Id: ltpda_pwelch.m,v 1.12 2007/07/16 12:52:21 ingo Exp $';
0039
0040
0041 bs = [];
0042
0043
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
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
0071 for j=1:na
0072
0073
0074 a = as(j);
0075 d = a.data;
0076
0077 if isa(a, 'ao')
0078
0079 if isa(d, 'tsdata')
0080
0081
0082 plo = combine(pls, getDefaultPL());
0083
0084
0085 Nfft = find(plo, 'Nfft');
0086 setWindow = 0;
0087 if Nfft <= 0
0088 Nfft = length(d.x);
0089 disp(sprintf('! Using default Nfft of %g', Nfft))
0090 setWindow = 1;
0091 end
0092 plo = set(plo, 'Nfft', Nfft);
0093
0094 Win = find(plo, 'Win');
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 = set(plo, 'Win', Win);
0105
0106 Nolap = find(plo, 'Nolap');
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 = set(plo, 'Nolap', Nolap);
0112
0113
0114 if length(Win.win) ~= Nfft
0115 error('### Window length should be same as Nfft.');
0116 end
0117
0118
0119 [f, pxx, info] = ltpda_psd(d.x, Win.win, Nolap, Nfft, d.fs, 'ASD', d.yunits);
0120
0121
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
0129 h = history(ALGONAME, VERSION, plo, a.hist);
0130 h = set(h, 'invars', invars);
0131
0132
0133 b = ao(fs, h);
0134
0135
0136 if na~=length(invars) || isempty(invars{j})
0137 nj = a.name;
0138 else
0139 nj = invars{j};
0140 end
0141 b = set(b, 'name', sprintf('PSD(%s)', nj));
0142
0143
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
0151
0152
0153 varargout{1} = bs;
0154
0155
0156
0157
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