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.15 2008/01/21 10:07:37 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.y), Kdes, Kmin, Jdes, Nolap, Debug);
0143
0144
0145 try
0146 if find(pl, 'M-FILE ONLY')
0147
0148 [Q, Qxx, ENBW] = mlpsd(d.y, f, r, m, L, d.fs, Win, Order, Nolap);
0149 else
0150 [Q, Qxx, ENBW] = mlpsd2(d.y, f, r, m, L, d.fs, Win, Order, Win.rov);
0151 end
0152 catch
0153
0154 [Q, Qxx, ENBW] = mlpsd(d.y, f, r, m, L, d.fs, Win, Order, Nolap);
0155 end
0156
0157
0158 fs = fsdata(f, sqrt(Qxx), d.fs);
0159 fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0160 fs = set(fs, 'yunits', [d.yunits '/\surdHz']);
0161 fs = set(fs, 'enbw', ENBW);
0162
0163
0164 h = history(ALGONAME, VERSION, pl, a.hist);
0165 h = set(h, 'invars', invars);
0166
0167
0168 b = ao(fs, h);
0169
0170
0171 if isempty(invars{i})
0172 n = a.name;
0173 else
0174 n = invars{i};
0175 end
0176 b = setnh(b, 'name', sprintf('lpsd(%s)', n));
0177 bs = [bs b];
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204 end
0205
0206 varargout{1} = bs;
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217 function varargout = mlpsd2(varargin)
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228 disp('%%%%%%%% Using ltpda_dft.mex to compute core DFT ');
0229
0230
0231 x = varargin{1};
0232 f = varargin{2};
0233 r = varargin{3};
0234 m = varargin{4};
0235 L = varargin{5};
0236 fs = varargin{6};
0237 win = varargin{7};
0238 order = varargin{8};
0239 olap = varargin{9};
0240
0241 twopi = 2.0*pi;
0242
0243 nx = length(x);
0244 nf = length(f);
0245 Sxx = zeros(nf,1);
0246 S = zeros(nf,1);
0247 ENBW = zeros(nf,1);
0248
0249 disp_each = round(nf/100)*10;
0250
0251 for j=1:nf
0252
0253 if mod(j,disp_each) == 0 || j == 1 || j == nf
0254 disp(sprintf('++ computing frequency %d of %d: %f Hz', j, nf, f(j)));
0255 end
0256
0257
0258 l = L(j);
0259 switch win.name
0260 case 'Kaiser'
0261 win = specwin(win.name, l, win.psll);
0262 otherwise
0263 win = specwin(win.name, l);
0264 end
0265
0266 p = 1i * twopi * m(j)/l.*(0:l-1);
0267 C = win.win .* exp(p);
0268
0269
0270 [A, nsegs] = ltpda_dft(x, l, C, olap, order);
0271
0272 if mod(j,disp_each) == 0 || j == 1 || j == nf
0273 disp(sprintf(' averaged %d segments', nsegs));
0274 end
0275
0276 A2ns = 2.0*A/nsegs;
0277 S1 = win.ws;
0278 S12 = S1*S1;
0279 S2 = win.ws2;
0280 ENBW(j) = fs*S2/S12;
0281 Sxx(j) = A2ns/fs/S2;
0282 S(j) = A2ns/S12;
0283
0284 end
0285
0286 varargout{1} = S;
0287 varargout{2} = Sxx;
0288 varargout{3} = ENBW;
0289
0290
0291
0292
0293 function varargout = mlpsd(varargin)
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304 x = varargin{1};
0305 f = varargin{2};
0306 r = varargin{3};
0307 m = varargin{4};
0308 L = varargin{5};
0309 fs = varargin{6};
0310 win = varargin{7};
0311 order = varargin{8};
0312 olap = varargin{9};
0313
0314 twopi = 2.0*pi;
0315
0316 nx = length(x);
0317 nf = length(f);
0318 Sxx = zeros(nf,1);
0319 S = zeros(nf,1);
0320 ENBW = zeros(nf,1);
0321
0322 disp_each = round(nf/100)*10;
0323
0324 for j=1:nf
0325
0326 if mod(j,disp_each) == 0 || j == 1 || j == nf
0327 disp(sprintf('++ computing frequency %d of %d: %f Hz', j, nf, f(j)));
0328 end
0329
0330
0331 l = L(j);
0332 switch win.name
0333 case 'Kaiser'
0334 win = specwin(win.name, l, win.psll);
0335 otherwise
0336 win = specwin(win.name, l);
0337 end
0338
0339 p = 1i * twopi * m(j)/l.*[0:l-1];
0340 C = win.win .* exp(p);
0341 Xolap = (1-olap);
0342
0343 A = 0.0;
0344 start = 1;
0345 nsegs = 0;
0346 while start + l - 1 <= nx
0347
0348
0349 xs = x(start:start+l-1);
0350
0351
0352 switch order
0353 case -1
0354
0355 case 0
0356 xs = xs - mean(xs);
0357 case 1
0358 xs = detrend(xs);
0359 otherwise
0360 xs = polydetrend(xs, order);
0361 end
0362
0363
0364 a = C*xs;
0365 A = A + a*conj(a);
0366
0367
0368 start = floor(start + l*Xolap);
0369 nsegs = nsegs + 1;
0370 end
0371
0372 if mod(j,disp_each) == 0 || j == 1 || j == nf
0373 disp(sprintf(' averaged %d segments', nsegs));
0374 end
0375
0376 A2ns = 2.0*A/nsegs;
0377 S1 = win.ws;
0378 S12 = S1*S1;
0379 S2 = win.ws2;
0380 ENBW(j) = fs*S2/S12;
0381 Sxx(j) = A2ns/fs/S2;
0382 S(j) = A2ns/S12;
0383
0384 end
0385
0386 varargout{1} = S;
0387 varargout{2} = Sxx;
0388 varargout{3} = ENBW;
0389
0390
0391
0392
0393
0394 function plo = getDefaultPL()
0395
0396 disp('* creating default plist...');
0397 plo = plist();
0398 plo = append(plo, param('Nfft', -1));
0399 plo = append(plo, param('Win', specwin('Kaiser', 10, 100)));
0400 plo = append(plo, param('Nolap', -1));
0401 plo = append(plo, param('Order', -1));
0402 disp('* done.');
0403
0404