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.8 2007/07/10 05:44:40 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.8 2007/07/10 05:44:40 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.8 2007/07/10 05:44:40 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     % use window function rov
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');    % desired overlap
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   % Get frequency vector
0140   [f,r,m,L,rr,rrr] = ltpda_compute_f(d.fs, length(d.x), Kdes, Kmin, Jdes, Nolap);
0141   
0142   % compute LPSD
0143   [Q, Qxx, ENBW] = mlpsd(d.x, f, r, m, L, d.fs, Win, Order);
0144 
0145   % create new output fsdata
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   % create new output history
0152   h = history(ALGONAME, VERSION, pl, a.hist);
0153   h = set(h, 'invars', invars);
0154 
0155   % make output analysis object
0156   b = ao(fs, h);
0157   
0158   % set name
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 %   % add AS as well
0168 %   if nargout == 2
0169 %     % create new output fsdata
0170 %     fs = fsdata(f, sqrt(Q), d.fs);
0171 %     fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0172 %     fs = set(fs, 'yunits', [d.yunits 'rms']);
0173 %     fs = set(fs, 'enbw', ENBW);
0174 %
0175 %     % create new output history
0176 %     h = history(ALGONAME, VERSION, pl, a.hist);
0177 %     h = set(h, 'invars', invars);
0178 %
0179 %     % make output analysis object
0180 %     b = ao(fs, h);
0181 %
0182 %     % set name
0183 %     if isempty(invars{i})
0184 %       n = a.name;
0185 %     else
0186 %       n = invars{i};
0187 %     end
0188 %     b = set(b, 'name', sprintf('lpsd(%s)', n));
0189 %     bxx = [bxx b];
0190 %   end
0191   
0192 end
0193 
0194 varargout{1} = bs;
0195 
0196 % include AS
0197 % if nargout > 1
0198 %   varargout{2} = bxx;
0199 % end
0200 
0201 %--------------------------------------------------------------------------
0202 %
0203 function varargout = mlpsd(varargin)
0204 
0205 % MLPSD
0206 %
0207 % [S,Sxx,ENBW] = mlpsd(x,f,r,m,L,fs,win,order)
0208 %
0209 % M Hewitson 19-05-07
0210 %
0211 %
0212 
0213 % Get inputs
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   % compute DFT exponent and window
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   % do segments
0252   A  = 0.0;
0253   start = 1;
0254   nsegs = 0;
0255   while start + l - 1 <= nx
0256 
0257     % get segment
0258     xs  = x(start:start+l-1);
0259     
0260     % detrend segment
0261     switch order
0262       case -1
0263         % do nothing
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     % make DFT
0273     a   = C*xs;
0274     A   = A + a*conj(a);
0275     
0276     % make next step
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 % Get default params
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 % END

Generated on Mon 03-Sep-2007 12:12:34 by m2html © 2003