Home > m > sigproc > frequency_domain > ltpda_lpsd.m

ltpda_lpsd

PURPOSE ^

LTPDA_LPSD implement LPSD algorithm for analysis objects.

SYNOPSIS ^

function varargout = ltpda_lpsd(varargin)

DESCRIPTION ^

 LTPDA_LPSD implement LPSD algorithm for analysis objects.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 LTPDA_LPSD implement LPSD algorithm for analysis objects.
 
   >> bs = ltpda_lpsd(as)   % returns ASD
 or
   >> bs = ltpda_lpsd(as, pl)   % returns ASD
 or
   >> [a, axx] = ltpda_lpsd(as, pl)  % returns AS and ASD
 
   Inputs:
     as  - array of analysis objects
     pl  - parameter list (see below)
 
   Outputs:
     bs  - array of analysis objects, one for each input
 
   Parameter list:
     Kdes  - desired number of averages   (default 100)
     Kmin  - minimum number of averages   (default 1)
     Jdes  - number of spectral frequencies to compute (default fs/2)
     win   - a specwin window object
             Only the design parameters of the window object are used; the
             window is recomputed for each DFT length inside the lpsd_core
             algorithm.
     Nolap - desired overlap (default: taken from window)
     Order - order of detrending
             -1 - no detrending
              0 - subtract mean
              1 - subtract linear fit
              N - subtract fit of polynomial, order N  
 
 The following call returns a parameter list object that contains the
 default parameter values:
 
 >> pl = ltpda_lpsd('Params')
 
 
 References: 
     "Improved spectrum estimation from digitized time series 
      on a logarithmic frequency axis", Michael Troebs, Gerhard Heinzel,
      Measurement 39 (2006) 120-129.
 
 M Hewitson 02-02-07
 
 $Id: ltpda_lpsd.m,v 1.15 2008/01/21 10:07:37 hewitson Exp $
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_lpsd(varargin)
0002 % LTPDA_LPSD implement LPSD algorithm for analysis objects.
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % LTPDA_LPSD implement LPSD algorithm for analysis objects.
0006 %
0007 %   >> bs = ltpda_lpsd(as)   % returns ASD
0008 % or
0009 %   >> bs = ltpda_lpsd(as, pl)   % returns ASD
0010 % or
0011 %   >> [a, axx] = ltpda_lpsd(as, pl)  % returns AS and ASD
0012 %
0013 %   Inputs:
0014 %     as  - array of analysis objects
0015 %     pl  - parameter list (see below)
0016 %
0017 %   Outputs:
0018 %     bs  - array of analysis objects, one for each input
0019 %
0020 %   Parameter list:
0021 %     Kdes  - desired number of averages   (default 100)
0022 %     Kmin  - minimum number of averages   (default 1)
0023 %     Jdes  - number of spectral frequencies to compute (default fs/2)
0024 %     win   - a specwin window object
0025 %             Only the design parameters of the window object are used; the
0026 %             window is recomputed for each DFT length inside the lpsd_core
0027 %             algorithm.
0028 %     Nolap - desired overlap (default: taken from window)
0029 %     Order - order of detrending
0030 %             -1 - no detrending
0031 %              0 - subtract mean
0032 %              1 - subtract linear fit
0033 %              N - subtract fit of polynomial, order N
0034 %
0035 % The following call returns a parameter list object that contains the
0036 % default parameter values:
0037 %
0038 % >> pl = ltpda_lpsd('Params')
0039 %
0040 %
0041 % References:
0042 %     "Improved spectrum estimation from digitized time series
0043 %      on a logarithmic frequency axis", Michael Troebs, Gerhard Heinzel,
0044 %      Measurement 39 (2006) 120-129.
0045 %
0046 % M Hewitson 02-02-07
0047 %
0048 % $Id: ltpda_lpsd.m,v 1.15 2008/01/21 10:07:37 hewitson Exp $
0049 %
0050 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0051 
0052 % Check if this is a call for parameters
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 % capture input variable names
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 % check plist
0080 if isempty(ps)
0081   pl = getDefaultPL();
0082 else
0083   pl = combine(ps, getDefaultPL);
0084 end
0085 
0086 % initialise output array
0087 bs  = []; 
0088 bxx = [];
0089 
0090 % Go through each input AO
0091 na = length(as);
0092 for i=1:na
0093   a = as(i);
0094   d = a.data;
0095   
0096   % check this is a time-series object
0097   if ~isa(d, 'tsdata')
0098     error('### lpsd requires tsdata (time-series) inputs.');
0099   end
0100   
0101   %--- check input parameters
0102   Kdes = find(pl, 'Kdes');   % desired averages
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');    % minimum averages
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');    % num desired spectral frequencies
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');    % desired overlap
0126   if isempty(Nolap) || Nolap == -1
0127     pl = remove(pl, 'Nolap'); % remove old value
0128     % use window function rov
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');    % desired overlap
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   % Get frequency vector
0142   [f,r,m,L,rr,rrr] = ltpda_compute_f(d.fs, length(d.y), Kdes, Kmin, Jdes, Nolap, Debug);
0143   
0144   % compute LPSD
0145   try
0146     if find(pl, 'M-FILE ONLY')
0147       % Using pure m-file version
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     % Using pure m-file version
0154     [Q, Qxx, ENBW] = mlpsd(d.y, f, r, m, L, d.fs, Win, Order, Nolap);
0155   end
0156   
0157   % create new output fsdata
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   % create new output history
0164   h = history(ALGONAME, VERSION, pl, a.hist);
0165   h = set(h, 'invars', invars);
0166 
0167   % make output analysis object
0168   b = ao(fs, h);
0169   
0170   % set name
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 %   % add AS as well
0180 %   if nargout == 2
0181 %     % create new output fsdata
0182 %     fs = fsdata(f, sqrt(Q), d.fs);
0183 %     fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0184 %     fs = set(fs, 'yunits', [d.yunits 'rms']);
0185 %     fs = set(fs, 'enbw', ENBW);
0186 %
0187 %     % create new output history
0188 %     h = history(ALGONAME, VERSION, pl, a.hist);
0189 %     h = set(h, 'invars', invars);
0190 %
0191 %     % make output analysis object
0192 %     b = ao(fs, h);
0193 %
0194 %     % set name
0195 %     if isempty(invars{i})
0196 %       n = a.name;
0197 %     else
0198 %       n = invars{i};
0199 %     end
0200 %     b = set(b, 'name', sprintf('lpsd(%s)', n));
0201 %     bxx = [bxx b];
0202 %   end
0203   
0204 end
0205 
0206 varargout{1} = bs;
0207 
0208 % include AS
0209 % if nargout > 1
0210 %   varargout{2} = bxx;
0211 % end
0212 
0213 
0214 
0215 %--------------------------------------------------------------------------
0216 %
0217 function varargout = mlpsd2(varargin)
0218 
0219 % MLPSD2
0220 %
0221 % [S,Sxx,ENBW] = mlpsd2(x,f,r,m,L,fs,win,order,olap)
0222 %
0223 % M Hewitson 19-05-07
0224 %
0225 %
0226 
0227 
0228 disp('%%%%%%%%   Using ltpda_dft.mex to compute core DFT ');
0229 
0230 % Get inputs
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   % compute DFT exponent and window
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   % Core DFT part in C-mex file
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 % MLPSD
0296 %
0297 % [S,Sxx,ENBW] = mlpsd(x,f,r,m,L,fs,win,order,olap)
0298 %
0299 % M Hewitson 19-05-07
0300 %
0301 %
0302 
0303 % Get inputs
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   % compute DFT exponent and window
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   % do segments
0343   A  = 0.0;
0344   start = 1;
0345   nsegs = 0;
0346   while start + l - 1 <= nx
0347 
0348     % get segment
0349     xs  = x(start:start+l-1);
0350     
0351     % detrend segment
0352     switch order
0353       case -1
0354         % do nothing
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     % make DFT
0364     a   = C*xs;
0365     A   = A + a*conj(a);
0366     
0367     % make next step
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 % Get default params
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 % END

Generated on Tue 22-Jan-2008 10:39:13 by m2html © 2003