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
 
 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.7 2007/07/02 08:28:02 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 %
0031 % The following call returns a parameter list object that contains the
0032 % default parameter values:
0033 %
0034 % >> pl = ltpda_lpsd('Params')
0035 %
0036 %
0037 % References:
0038 %     "Improved spectrum estimation from digitized time series
0039 %      on a logarithmic frequency axis", Michael Troebs, Gerhard Heinzel,
0040 %      Measurement 39 (2006) 120-129.
0041 %
0042 % M Hewitson 02-02-07
0043 %
0044 % $Id: ltpda_lpsd.m,v 1.7 2007/07/02 08:28:02 hewitson Exp $
0045 %
0046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0047 
0048 % Check if this is a call for parameters
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 % capture input variable names
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 % check plist
0076 if isempty(ps)
0077   pl = getDefaultPL();
0078 else
0079   pl = combine(ps, getDefaultPL);
0080 end
0081 
0082 % initialise output array
0083 bs  = []; 
0084 bxx = [];
0085 
0086 % Go through each input AO
0087 na = length(as);
0088 for i=1:na
0089   a = as(i);
0090   d = a.data;
0091   
0092   % check this is a time-series object
0093   if ~isa(d, 'tsdata')
0094     error('### lpsd requires tsdata (time-series) inputs.');
0095   end
0096   
0097   %--- check input parameters
0098   Kdes = find(pl, 'Kdes');   % desired averages
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');    % minimum averages
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');    % num desired spectral frequencies
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');    % desired overlap
0122   if isempty(Nolap) || Nolap == -1
0123     % use window function rov
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');    % desired overlap
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   % Get frequency vector
0136   [f,r,m,L,rr,rrr] = ltpda_compute_f(d.fs, length(d.x), Kdes, Kmin, Jdes, Nolap);
0137   
0138   % compute LPSD
0139   [Q, Qxx, ENBW] = mlpsd(d.x, f, r, m, L, d.fs, Win, Order);
0140 
0141   % create new output fsdata
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   % create new output history
0148   h = history(ALGONAME, VERSION, pl, a.hist);
0149   h = set(h, 'invars', invars);
0150 
0151   % make output analysis object
0152   b = ao(fs, h);
0153   
0154   % set name
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 %   % add AS as well
0164 %   if nargout == 2
0165 %     % create new output fsdata
0166 %     fs = fsdata(f, sqrt(Q), d.fs);
0167 %     fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0168 %     fs = set(fs, 'yunits', [d.yunits 'rms']);
0169 %     fs = set(fs, 'enbw', ENBW);
0170 %
0171 %     % create new output history
0172 %     h = history(ALGONAME, VERSION, pl, a.hist);
0173 %     h = set(h, 'invars', invars);
0174 %
0175 %     % make output analysis object
0176 %     b = ao(fs, h);
0177 %
0178 %     % set name
0179 %     if isempty(invars{i})
0180 %       n = a.name;
0181 %     else
0182 %       n = invars{i};
0183 %     end
0184 %     b = set(b, 'name', sprintf('lpsd(%s)', n));
0185 %     bxx = [bxx b];
0186 %   end
0187   
0188 end
0189 
0190 varargout{1} = bs;
0191 
0192 % include AS
0193 % if nargout > 1
0194 %   varargout{2} = bxx;
0195 % end
0196 
0197 %--------------------------------------------------------------------------
0198 %
0199 function varargout = mlpsd(varargin)
0200 
0201 % MLPSD
0202 %
0203 % [S,Sxx,ENBW] = mlpsd(x,f,r,m,L,fs,win,order)
0204 %
0205 % M Hewitson 19-05-07
0206 %
0207 %
0208 
0209 % Get inputs
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   % compute DFT exponent and window
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   % do segments
0248   A  = 0.0;
0249   start = 1;
0250   nsegs = 0;
0251   while start + l - 1 <= nx
0252 
0253     % get segment
0254     xs  = x(start:start+l-1);
0255     
0256     % detrend segment
0257     switch order
0258       case -1
0259         % do nothing
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     % make DFT
0269     a   = C*xs;
0270     A   = A + a*conj(a);
0271     
0272     % make next step
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 % Get default params
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 % END

Generated on Mon 02-Jul-2007 12:19:41 by m2html © 2003