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