0001 function bs = 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 ALGONAME = mfilename;
0036 VERSION = '$Id: ltpda_lpsd.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $';
0037
0038
0039 invars = {};
0040 as = [];
0041 pl = [];
0042 for j=1:nargin
0043 invars = [invars cellstr(inputname(j))];
0044 if isa(varargin{j}, 'ao')
0045 as = [as varargin{j}];
0046 end
0047 if isa(varargin{j}, 'plist')
0048 pl = varargin{j};
0049 end
0050 end
0051
0052
0053 bs = [];
0054
0055
0056 if isempty(pl)
0057 pl = plist();
0058 end
0059
0060
0061 na = length(as);
0062 for i=1:na
0063 a = as(i);
0064 d = a.data;
0065
0066
0067 if ~isa(d, 'tsdata')
0068 error('### lpsd requires tsdata (time-series) inputs.');
0069 end
0070
0071
0072
0073
0074 Kdes = find(pl, 'Kdes');
0075 if isempty(Kdes)
0076 Kdes = 100;
0077 disp(sprintf('! Using default Kdes of %d', Kdes))
0078 pl = append(pl, param('Kdes', Kdes));
0079 end
0080 Kmin = find(pl, 'Kmin');
0081 if isempty(Kmin)
0082 Kmin = 1;
0083 disp(sprintf('! Using default Kmin of %d', Kmin))
0084 pl = append(pl, param('Kmin', Kmin));
0085 end
0086
0087 Jdes = find(pl, 'Jdes');
0088 if isempty(Jdes)
0089 Jdes = round(d.fs/4);
0090 disp(sprintf('! Using default Jdes of %d', Jdes))
0091 pl = append(pl, param('Jdes', Jdes));
0092 end
0093
0094 Win = find(pl, 'Win');
0095 if isempty(Win)
0096 Win = specwin('Kaiser', 100, 250);
0097 pl = append(pl, param('Win', Win));
0098 end
0099
0100 Nolap = find(pl, 'Nolap');
0101 if isempty(Nolap)
0102
0103 Nolap = Win.rov/100;
0104 disp(sprintf('! Using recommended Nolap of %d', Nolap))
0105 pl = append(pl, param('Nolap', Nolap));
0106 end
0107
0108 Order = find(pl, 'Order');
0109 if isempty(Order)
0110 Order = 0;
0111 disp(sprintf('! Using default detrending order of %d', Order))
0112 pl = append(pl, param('Order', Order));
0113 end
0114
0115
0116 [f,r,m,L,rr,rrr] = ltpda_compute_f(d.fs, length(d.x), Kdes, Kmin, Jdes, Nolap);
0117
0118
0119 [Q, Qxx, ENBW] = mlpsd(d.x, f, r, m, L, d.fs, Win, Order);
0120
0121
0122 fs = fsdata(f, sqrt(Qxx), d.fs);
0123 fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0124 fs = set(fs, 'yunits', [d.yunits '/\surdHz']);
0125 fs = set(fs, 'enbw', ENBW);
0126
0127
0128 h = history(ALGONAME, VERSION, pl, a.hist);
0129 h = set(h, 'invars', invars);
0130
0131
0132 b = ao(fs, h);
0133
0134
0135 if isempty(invars{1})
0136 n = a.name;
0137 else
0138 n = invars{1};
0139 end
0140 b = set(b, 'name', sprintf('lpsd(%s)', n));
0141
0142 bs = [bs b];
0143 end
0144
0145
0146
0147 function varargout = mlpsd(varargin)
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158 x = varargin{1};
0159 f = varargin{2};
0160 r = varargin{3};
0161 m = varargin{4};
0162 L = varargin{5};
0163 fs = varargin{6};
0164 win = varargin{7};
0165 order = varargin{8};
0166
0167 twopi = 2.0*pi;
0168
0169 nx = length(x);
0170 nf = length(f);
0171 Sxx = zeros(nf,1);
0172 S = zeros(nf,1);
0173 ENBW = zeros(nf,1);
0174
0175 disp_each = round(nf/100)*10;
0176
0177 for j=1:nf
0178
0179 if mod(j,disp_each) == 0 || j == 1 || j == nf
0180 disp(sprintf('++ computing frequency %d of %d: %f Hz', j, nf, f(j)));
0181 end
0182
0183
0184 l = L(j);
0185 switch win.name
0186 case 'Kaiser'
0187 win = specwin(win.name, l, win.psll);
0188 otherwise
0189 win = specwin(win.name, l);
0190 end
0191
0192 p = 1i * twopi * m(j)/l.*[0:l-1];
0193 C = win.win .* exp(p);
0194 Xolap = (1-win.rov/100);
0195
0196 A = 0.0;
0197 start = 1;
0198 nsegs = 0;
0199 while start + l - 1 <= nx
0200
0201
0202 xs = x(start:start+l-1);
0203
0204
0205 switch order
0206 case -1
0207
0208 case 0
0209 xs = xs - mean(xs);
0210 case 1
0211 xs = detrend(xs);
0212 otherwise
0213 xs = polydetrend(xs, order);
0214 end
0215
0216
0217 a = C*xs;
0218 A = A + a*conj(a);
0219
0220
0221 start = floor(start + l*Xolap);
0222 nsegs = nsegs + 1;
0223 end
0224
0225 if mod(j,disp_each) == 0 || j == 1 || j == nf
0226 disp(sprintf(' averaged %d segments', nsegs));
0227 end
0228
0229 A2ns = 2.0*A/nsegs;
0230 S1 = win.ws;
0231 S12 = S1*S1;
0232 S2 = win.ws2;
0233 ENBW(j) = fs*S2/S12;
0234 Sxx(j) = A2ns/fs/S2;
0235 S(j) = A2ns/S12;
0236
0237 end
0238
0239 varargout{1} = S;
0240 varargout{2} = Sxx;
0241 varargout{3} = ENBW;
0242
0243
0244
0245
0246