0001 function varargout = ltpda_lpsd(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
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053 if nargin == 1
0054 in = char(varargin{1});
0055 if strcmp(in, 'Params')
0056 varargout{1} = getDefaultPL();
0057 return
0058 end
0059 end
0060
0061
0062 invars = {};
0063 as = [];
0064 ps = [];
0065 for j=1:nargin
0066 invars = [invars cellstr(inputname(j))];
0067 if isa(varargin{j}, 'ao')
0068 as = [as varargin{j}];
0069 end
0070 if isa(varargin{j}, 'plist')
0071 ps = [ps varargin{j}];
0072 end
0073 end
0074
0075
0076 ALGONAME = mfilename;
0077 VERSION = '$Id: ltpda_lpsd.m,v 1.10 2007/10/28 10:23:40 hewitson Exp $';
0078
0079
0080 if isempty(ps)
0081 pl = getDefaultPL();
0082 else
0083 pl = combine(ps, getDefaultPL);
0084 end
0085
0086
0087 bs = [];
0088 bxx = [];
0089
0090
0091 na = length(as);
0092 for i=1:na
0093 a = as(i);
0094 d = a.data;
0095
0096
0097 if ~isa(d, 'tsdata')
0098 error('### lpsd requires tsdata (time-series) inputs.');
0099 end
0100
0101
0102 Kdes = find(pl, 'Kdes');
0103 if isempty(Kdes)
0104 Kdes = 100;
0105 disp(sprintf('! Using default Kdes of %d', Kdes))
0106 pl = append(pl, param('Kdes', Kdes));
0107 end
0108 Kmin = find(pl, 'Kmin');
0109 if isempty(Kmin)
0110 Kmin = 1;
0111 disp(sprintf('! Using default Kmin of %d', Kmin))
0112 pl = append(pl, param('Kmin', Kmin));
0113 end
0114 Jdes = find(pl, 'Jdes');
0115 if isempty(Jdes)
0116 Jdes = round(d.fs/4);
0117 disp(sprintf('! Using default Jdes of %d', Jdes))
0118 pl = append(pl, param('Jdes', Jdes));
0119 end
0120 Win = find(pl, 'Win');
0121 if isempty(Win)
0122 Win = specwin('Kaiser', 100, 250);
0123 pl = append(pl, param('Win', Win));
0124 end
0125 Nolap = find(pl, 'Nolap');
0126 if isempty(Nolap) || Nolap == -1
0127 pl = remove(pl, 'Nolap');
0128
0129 Nolap = Win.rov/100;
0130 disp(sprintf('! Using recommended Nolap of %d', Nolap))
0131 pl = append(pl, param('Nolap', Nolap));
0132 end
0133 Order = find(pl, 'Order');
0134 if isempty(Order)
0135 Order = 0;
0136 disp(sprintf('! Using default detrending order of %d', Order))
0137 pl = append(pl, param('Order', Order));
0138 end
0139 Debug = find(pl, 'DEBUG');
0140
0141
0142 [f,r,m,L,rr,rrr] = ltpda_compute_f(d.fs, length(d.x), Kdes, Kmin, Jdes, Nolap, Debug);
0143
0144
0145 [Q, Qxx, ENBW] = mlpsd(d.x, f, r, m, L, d.fs, Win, Order);
0146
0147
0148 fs = fsdata(f, sqrt(Qxx), d.fs);
0149 fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0150 fs = set(fs, 'yunits', [d.yunits '/\surdHz']);
0151 fs = set(fs, 'enbw', ENBW);
0152
0153
0154 h = history(ALGONAME, VERSION, pl, a.hist);
0155 h = set(h, 'invars', invars);
0156
0157
0158 b = ao(fs, h);
0159
0160
0161 if isempty(invars{i})
0162 n = a.name;
0163 else
0164 n = invars{i};
0165 end
0166 b = set(b, 'name', sprintf('lpsd(%s)', n));
0167 bs = [bs b];
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194 end
0195
0196 varargout{1} = bs;
0197
0198
0199
0200
0201
0202
0203
0204
0205 function varargout = mlpsd(varargin)
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216 x = varargin{1};
0217 f = varargin{2};
0218 r = varargin{3};
0219 m = varargin{4};
0220 L = varargin{5};
0221 fs = varargin{6};
0222 win = varargin{7};
0223 order = varargin{8};
0224
0225 twopi = 2.0*pi;
0226
0227 nx = length(x);
0228 nf = length(f);
0229 Sxx = zeros(nf,1);
0230 S = zeros(nf,1);
0231 ENBW = zeros(nf,1);
0232
0233 disp_each = round(nf/100)*10;
0234
0235 for j=1:nf
0236
0237 if mod(j,disp_each) == 0 || j == 1 || j == nf
0238 disp(sprintf('++ computing frequency %d of %d: %f Hz', j, nf, f(j)));
0239 end
0240
0241
0242 l = L(j);
0243 switch win.name
0244 case 'Kaiser'
0245 win = specwin(win.name, l, win.psll);
0246 otherwise
0247 win = specwin(win.name, l);
0248 end
0249
0250 p = 1i * twopi * m(j)/l.*[0:l-1];
0251 C = win.win .* exp(p);
0252 Xolap = (1-win.rov/100);
0253
0254 A = 0.0;
0255 start = 1;
0256 nsegs = 0;
0257 while start + l - 1 <= nx
0258
0259
0260 xs = x(start:start+l-1);
0261
0262
0263 switch order
0264 case -1
0265
0266 case 0
0267 xs = xs - mean(xs);
0268 case 1
0269 xs = detrend(xs);
0270 otherwise
0271 xs = polydetrend(xs, order);
0272 end
0273
0274
0275 a = C*xs;
0276 A = A + a*conj(a);
0277
0278
0279 start = floor(start + l*Xolap);
0280 nsegs = nsegs + 1;
0281 end
0282
0283 if mod(j,disp_each) == 0 || j == 1 || j == nf
0284 disp(sprintf(' averaged %d segments', nsegs));
0285 end
0286
0287 A2ns = 2.0*A/nsegs;
0288 S1 = win.ws;
0289 S12 = S1*S1;
0290 S2 = win.ws2;
0291 ENBW(j) = fs*S2/S12;
0292 Sxx(j) = A2ns/fs/S2;
0293 S(j) = A2ns/S12;
0294
0295 end
0296
0297 varargout{1} = S;
0298 varargout{2} = Sxx;
0299 varargout{3} = ENBW;
0300
0301
0302
0303
0304
0305 function plo = getDefaultPL()
0306
0307 disp('* creating default plist...');
0308 plo = plist();
0309 plo = append(plo, param('Nfft', -1));
0310 plo = append(plo, param('Win', specwin('Kaiser', 10, 100)));
0311 plo = append(plo, param('Nolap', -1));
0312 plo = append(plo, param('Order', 0));
0313 disp('* done.');
0314
0315