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.8 2007/07/10 05:44: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
0128 Nolap = Win.rov/100;
0129 disp(sprintf('! Using recommended Nolap of %d', Nolap))
0130 pl = append(pl, param('Nolap', Nolap));
0131 end
0132 Order = find(pl, 'Order');
0133 if isempty(Order)
0134 Order = 0;
0135 disp(sprintf('! Using default detrending order of %d', Order))
0136 pl = append(pl, param('Order', Order));
0137 end
0138
0139
0140 [f,r,m,L,rr,rrr] = ltpda_compute_f(d.fs, length(d.x), Kdes, Kmin, Jdes, Nolap);
0141
0142
0143 [Q, Qxx, ENBW] = mlpsd(d.x, f, r, m, L, d.fs, Win, Order);
0144
0145
0146 fs = fsdata(f, sqrt(Qxx), d.fs);
0147 fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0148 fs = set(fs, 'yunits', [d.yunits '/\surdHz']);
0149 fs = set(fs, 'enbw', ENBW);
0150
0151
0152 h = history(ALGONAME, VERSION, pl, a.hist);
0153 h = set(h, 'invars', invars);
0154
0155
0156 b = ao(fs, h);
0157
0158
0159 if isempty(invars{i})
0160 n = a.name;
0161 else
0162 n = invars{i};
0163 end
0164 b = set(b, 'name', sprintf('lpsd(%s)', n));
0165 bs = [bs b];
0166
0167
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 end
0193
0194 varargout{1} = bs;
0195
0196
0197
0198
0199
0200
0201
0202
0203 function varargout = mlpsd(varargin)
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214 x = varargin{1};
0215 f = varargin{2};
0216 r = varargin{3};
0217 m = varargin{4};
0218 L = varargin{5};
0219 fs = varargin{6};
0220 win = varargin{7};
0221 order = varargin{8};
0222
0223 twopi = 2.0*pi;
0224
0225 nx = length(x);
0226 nf = length(f);
0227 Sxx = zeros(nf,1);
0228 S = zeros(nf,1);
0229 ENBW = zeros(nf,1);
0230
0231 disp_each = round(nf/100)*10;
0232
0233 for j=1:nf
0234
0235 if mod(j,disp_each) == 0 || j == 1 || j == nf
0236 disp(sprintf('++ computing frequency %d of %d: %f Hz', j, nf, f(j)));
0237 end
0238
0239
0240 l = L(j);
0241 switch win.name
0242 case 'Kaiser'
0243 win = specwin(win.name, l, win.psll);
0244 otherwise
0245 win = specwin(win.name, l);
0246 end
0247
0248 p = 1i * twopi * m(j)/l.*[0:l-1];
0249 C = win.win .* exp(p);
0250 Xolap = (1-win.rov/100);
0251
0252 A = 0.0;
0253 start = 1;
0254 nsegs = 0;
0255 while start + l - 1 <= nx
0256
0257
0258 xs = x(start:start+l-1);
0259
0260
0261 switch order
0262 case -1
0263
0264 case 0
0265 xs = xs - mean(xs);
0266 case 1
0267 xs = detrend(xs);
0268 otherwise
0269 xs = polydetrend(xs, order);
0270 end
0271
0272
0273 a = C*xs;
0274 A = A + a*conj(a);
0275
0276
0277 start = floor(start + l*Xolap);
0278 nsegs = nsegs + 1;
0279 end
0280
0281 if mod(j,disp_each) == 0 || j == 1 || j == nf
0282 disp(sprintf(' averaged %d segments', nsegs));
0283 end
0284
0285 A2ns = 2.0*A/nsegs;
0286 S1 = win.ws;
0287 S12 = S1*S1;
0288 S2 = win.ws2;
0289 ENBW(j) = fs*S2/S12;
0290 Sxx(j) = A2ns/fs/S2;
0291 S(j) = A2ns/S12;
0292
0293 end
0294
0295 varargout{1} = S;
0296 varargout{2} = Sxx;
0297 varargout{3} = ENBW;
0298
0299
0300
0301
0302
0303 function plo = getDefaultPL()
0304
0305 disp('* creating default plist...');
0306 plo = plist();
0307 plo = append(plo, param('Nfft', -1));
0308 plo = append(plo, param('Win', specwin('Kaiser', 10, 100)));
0309 plo = append(plo, param('Nolap', -1));
0310 plo = append(plo, param('Order', 0));
0311 disp('* done.');
0312
0313