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
0054
0055
0056
0057
0058
0059
0060 ALGONAME = mfilename;
0061 VERSION = '$Id: ltpda_lpsd.html,v 1.14 2008/03/31 10:27:39 hewitson Exp $';
0062 CATEGORY = 'Signal Processing';
0063
0064
0065
0066
0067 if nargin == 1 && ischar(varargin{1})
0068 in = char(varargin{1});
0069 if strcmp(in, 'Params')
0070 varargout{1} = getDefaultPL();
0071 return
0072 elseif strcmp(in, 'Version')
0073 varargout{1} = VERSION;
0074 return
0075 elseif strcmp(in, 'Category')
0076 varargout{1} = CATEGORY;
0077 return
0078 end
0079 end
0080
0081
0082 invars = {};
0083 as = [];
0084 ps = [];
0085 for j=1:nargin
0086 invars = [invars cellstr(inputname(j))];
0087 if isa(varargin{j}, 'ao')
0088 as = [as varargin{j}];
0089 end
0090 if isa(varargin{j}, 'plist')
0091 ps = [ps varargin{j}];
0092 end
0093 end
0094
0095
0096
0097
0098 if isempty(ps)
0099 pl = getDefaultPL();
0100 else
0101 pl = combine(ps, getDefaultPL);
0102 end
0103
0104
0105 bs = [];
0106 bxx = [];
0107
0108
0109 na = length(as);
0110 for i=1:na
0111 a = as(i);
0112 d = a.data;
0113
0114
0115 if ~isa(d, 'tsdata')
0116 error('### lpsd requires tsdata (time-series) inputs.');
0117 end
0118
0119
0120 Kdes = find(pl, 'Kdes');
0121 if isempty(Kdes)
0122 Kdes = 100;
0123 disp(sprintf('! Using default Kdes of %d', Kdes))
0124 pl = append(pl, param('Kdes', Kdes));
0125 end
0126
0127
0128 Kmin = find(pl, 'Kmin');
0129 if isempty(Kmin)
0130 Kmin = 1;
0131 disp(sprintf('! Using default Kmin of %d', Kmin))
0132 pl = append(pl, param('Kmin', Kmin));
0133 end
0134
0135
0136 Jdes = find(pl, 'Jdes');
0137 if isempty(Jdes)
0138 Jdes = round(d.fs/4);
0139 disp(sprintf('! Using default Jdes of %d', Jdes))
0140 pl = append(pl, param('Jdes', Jdes));
0141 end
0142 Win = find(pl, 'Win');
0143 if isempty(Win)
0144 Win = specwin('Kaiser', 100, 250);
0145 pl = append(pl, param('Win', Win));
0146 end
0147 Nolap = find(pl, 'Nolap');
0148 if isempty(Nolap) || Nolap == -1
0149 pl = remove(pl, 'Nolap');
0150
0151 Nolap = Win.rov/100;
0152 disp(sprintf('! Using recommended Nolap of %d', Nolap))
0153 pl = append(pl, param('Nolap', Nolap));
0154 end
0155 Order = find(pl, 'Order');
0156 if isempty(Order)
0157 Order = 0;
0158 disp(sprintf('! Using default detrending order of %d', Order))
0159 pl = append(pl, param('Order', Order));
0160 end
0161 Lmin = find(pl, 'Lmin');
0162 if isempty(Lmin)
0163 Lmin = 0;
0164 end
0165 Debug = find(pl, 'DEBUG');
0166
0167
0168 if find(pl, 'OLD SCHEDULER')
0169 [f,r,m,L,rr,rrr] = ltpda_compute_f(d.fs, length(d.y), Kdes, Kmin, Jdes, Nolap, Debug);
0170 else
0171 [f, r, m, L, K] = ltpda_ltf_plan(length(d.y), d.fs, Nolap, 1, Lmin, Jdes, Kdes);
0172 end
0173
0174
0175 try
0176 if find(pl, 'M-FILE ONLY')
0177
0178 [P, Pxx, ENBW] = mlpsd(d.y, f, r, m, L, d.fs, Win, Order, Nolap);
0179 else
0180 [P, Pxx, ENBW] = mlpsd2(d.y, f, r, m, L, d.fs, Win, Order, Win.rov, Lmin);
0181 end
0182 catch
0183
0184 [P, Pxx, ENBW] = mlpsd(d.y, f, r, m, L, d.fs, Win, Order, Nolap);
0185 end
0186
0187
0188 scale = find(pl, 'Scale');
0189 switch scale
0190 case 'AS'
0191 fs = fsdata(f, sqrt(P), d.fs);
0192 fs = set(fs, 'yunits', [d.yunits]);
0193 case 'ASD'
0194 fs = fsdata(f, sqrt(Pxx), d.fs);
0195 fs = set(fs, 'yunits', [d.yunits '/\surdHz']);
0196 case 'PS'
0197 fs = fsdata(f, P, d.fs);
0198 fs = set(fs, 'yunits', [d.yunits '^2']);
0199 case 'PSD'
0200 fs = fsdata(f, Pxx, d.fs);
0201 fs = set(fs, 'yunits', [d.yunits '^2/Hz']);
0202 otherwise
0203 error(['### Unknown scaling:' scale]);
0204 end
0205 fs = set(fs, 'name', sprintf('lpsd(%s)', d.name));
0206 fs = set(fs, 'enbw', ENBW);
0207
0208
0209 h = history(ALGONAME, VERSION, pl, a.hist);
0210 h = set(h, 'invars', invars);
0211
0212
0213 b = ao(fs, h);
0214
0215
0216 if isempty(invars{i})
0217 n = a.name;
0218 else
0219 n = invars{i};
0220 end
0221 b = setnh(b, 'name', sprintf('lpsd(%s)', n), ...
0222 'mfilename', ALGONAME);
0223 bs = [bs b];
0224
0225 end
0226
0227 varargout{1} = bs;
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238 function varargout = mlpsd2(varargin)
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249 disp('%%%%%%%% Using ltpda_dft.mex to compute core DFT ');
0250
0251
0252 x = varargin{1};
0253 f = varargin{2};
0254 r = varargin{3};
0255 m = varargin{4};
0256 L = varargin{5};
0257 fs = varargin{6};
0258 win = varargin{7};
0259 order = varargin{8};
0260 olap = varargin{9};
0261 Lmin = varargin{10};
0262
0263 twopi = 2.0*pi;
0264
0265 nx = length(x);
0266 nf = length(f);
0267 Sxx = zeros(nf,1);
0268 S = zeros(nf,1);
0269 ENBW = zeros(nf,1);
0270
0271 disp_each = round(nf/100)*10;
0272
0273 minReached = 0;
0274
0275 for j=1:nf
0276
0277 if mod(j,disp_each) == 0 || j == 1 || j == nf
0278 disp(sprintf('++ computing frequency %04d of %d: %f Hz', j, nf, f(j)));
0279 end
0280
0281
0282 l = L(j);
0283
0284
0285
0286 if ~minReached
0287 switch win.name
0288 case 'Kaiser'
0289 win = specwin(win.name, l, win.psll);
0290 otherwise
0291 win = specwin(win.name, l);
0292 end
0293 if l==Lmin
0294 minReached = 1;
0295 end
0296 end
0297
0298 p = 1i * twopi * m(j)/l.*(0:l-1);
0299 C = win.win .* exp(p);
0300
0301
0302 [A, nsegs] = ltpda_dft(x, l, C, olap, order);
0303
0304 if mod(j,disp_each) == 0 || j == 1 || j == nf
0305 disp(sprintf(' averaged %d segments', nsegs));
0306 end
0307
0308 A2ns = 2.0*A/nsegs;
0309 S1 = win.ws;
0310 S12 = S1*S1;
0311 S2 = win.ws2;
0312 ENBW(j) = fs*S2/S12;
0313 Sxx(j) = A2ns/fs/S2;
0314 S(j) = A2ns/S12;
0315
0316 end
0317
0318 varargout{1} = S;
0319 varargout{2} = Sxx;
0320 varargout{3} = ENBW;
0321
0322
0323
0324
0325 function varargout = mlpsd(varargin)
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336 x = varargin{1};
0337 f = varargin{2};
0338 r = varargin{3};
0339 m = varargin{4};
0340 L = varargin{5};
0341 fs = varargin{6};
0342 win = varargin{7};
0343 order = varargin{8};
0344 olap = varargin{9};
0345
0346 twopi = 2.0*pi;
0347
0348 nx = length(x);
0349 nf = length(f);
0350 Sxx = zeros(nf,1);
0351 S = zeros(nf,1);
0352 ENBW = zeros(nf,1);
0353
0354 disp_each = round(nf/100)*10;
0355
0356 for j=1:nf
0357
0358 if mod(j,disp_each) == 0 || j == 1 || j == nf
0359 disp(sprintf('++ computing frequency %d of %d: %f Hz', j, nf, f(j)));
0360 end
0361
0362
0363 l = L(j);
0364 switch win.name
0365 case 'Kaiser'
0366 win = specwin(win.name, l, win.psll);
0367 otherwise
0368 win = specwin(win.name, l);
0369 end
0370
0371 p = 1i * twopi * m(j)/l.*[0:l-1];
0372 C = win.win .* exp(p);
0373 Xolap = (1-olap);
0374
0375 A = 0.0;
0376 start = 1;
0377 nsegs = 0;
0378 while start + l - 1 <= nx
0379
0380
0381 xs = x(start:start+l-1);
0382
0383
0384 switch order
0385 case -1
0386
0387 case 0
0388 xs = xs - mean(xs);
0389 case 1
0390 xs = detrend(xs);
0391 otherwise
0392 xs = polydetrend(xs, order);
0393 end
0394
0395
0396 a = C*xs;
0397 A = A + a*conj(a);
0398
0399
0400 start = floor(start + l*Xolap);
0401 nsegs = nsegs + 1;
0402 end
0403
0404 if mod(j,disp_each) == 0 || j == 1 || j == nf
0405 disp(sprintf(' averaged %d segments', nsegs));
0406 end
0407
0408 A2ns = 2.0*A/nsegs;
0409 S1 = win.ws;
0410 S12 = S1*S1;
0411 S2 = win.ws2;
0412 ENBW(j) = fs*S2/S12;
0413 Sxx(j) = A2ns/fs/S2;
0414 S(j) = A2ns/S12;
0415
0416 end
0417
0418 varargout{1} = S;
0419 varargout{2} = Sxx;
0420 varargout{3} = ENBW;
0421
0422
0423
0424
0425
0426 function plo = getDefaultPL()
0427
0428 disp('* creating default plist...');
0429 plo = plist('Kdes', 100, ...
0430 'Kmin', 1, ...
0431 'Lmin', 0, ...
0432 'Jdes', 100, ...
0433 'Win', specwin('Kaiser', 10, 200), ...
0434 'Olap', -1, ...
0435 'Order', -1, ...
0436 'Scale', 'PSD');
0437
0438 disp('* done.');
0439
0440