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