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.10 2007/10/28 10:23: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.10 2007/10/28 10:23: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.10 2007/10/28 10:23: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     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.x), Kdes, Kmin, Jdes, Nolap, Debug);
0143   
0144   % compute LPSD
0145   [Q, Qxx, ENBW] = mlpsd(d.x, f, r, m, L, d.fs, Win, Order);
0146 
0147   % create new output fsdata
0148   fs = fsdata(f, sqrt(Qxx), d.fs);
0149   fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0150   fs = set(fs, 'yunits', [d.yunits '/\surdHz']);
0151   fs = set(fs, 'enbw', ENBW);
0152   
0153   % create new output history
0154   h = history(ALGONAME, VERSION, pl, a.hist);
0155   h = set(h, 'invars', invars);
0156 
0157   % make output analysis object
0158   b = ao(fs, h);
0159   
0160   % set name
0161   if isempty(invars{i})
0162     n = a.name;
0163   else
0164     n = invars{i};
0165   end
0166   b = set(b, 'name', sprintf('lpsd(%s)', n));  
0167   bs = [bs b];
0168   
0169 %   % add AS as well
0170 %   if nargout == 2
0171 %     % create new output fsdata
0172 %     fs = fsdata(f, sqrt(Q), d.fs);
0173 %     fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0174 %     fs = set(fs, 'yunits', [d.yunits 'rms']);
0175 %     fs = set(fs, 'enbw', ENBW);
0176 %
0177 %     % create new output history
0178 %     h = history(ALGONAME, VERSION, pl, a.hist);
0179 %     h = set(h, 'invars', invars);
0180 %
0181 %     % make output analysis object
0182 %     b = ao(fs, h);
0183 %
0184 %     % set name
0185 %     if isempty(invars{i})
0186 %       n = a.name;
0187 %     else
0188 %       n = invars{i};
0189 %     end
0190 %     b = set(b, 'name', sprintf('lpsd(%s)', n));
0191 %     bxx = [bxx b];
0192 %   end
0193   
0194 end
0195 
0196 varargout{1} = bs;
0197 
0198 % include AS
0199 % if nargout > 1
0200 %   varargout{2} = bxx;
0201 % end
0202 
0203 %--------------------------------------------------------------------------
0204 %
0205 function varargout = mlpsd(varargin)
0206 
0207 % MLPSD
0208 %
0209 % [S,Sxx,ENBW] = mlpsd(x,f,r,m,L,fs,win,order)
0210 %
0211 % M Hewitson 19-05-07
0212 %
0213 %
0214 
0215 % Get inputs
0216 x     = varargin{1};
0217 f     = varargin{2};
0218 r     = varargin{3};
0219 m     = varargin{4};
0220 L     = varargin{5};
0221 fs    = varargin{6};
0222 win   = varargin{7};
0223 order = varargin{8};
0224 
0225 twopi    = 2.0*pi;
0226 
0227 nx   = length(x);
0228 nf   = length(f);
0229 Sxx  = zeros(nf,1);
0230 S    = zeros(nf,1);
0231 ENBW = zeros(nf,1);
0232 
0233 disp_each = round(nf/100)*10;
0234 
0235 for j=1:nf
0236   
0237   if mod(j,disp_each) == 0 || j == 1 || j == nf
0238     disp(sprintf('++ computing frequency %d of %d: %f Hz', j, nf, f(j)));
0239   end
0240 
0241   % compute DFT exponent and window
0242   l = L(j);
0243   switch win.name
0244     case 'Kaiser'
0245       win = specwin(win.name, l, win.psll);
0246     otherwise
0247       win = specwin(win.name, l);
0248   end
0249 
0250   p     = 1i * twopi * m(j)/l.*[0:l-1];
0251   C     = win.win .* exp(p);
0252   Xolap = (1-win.rov/100);
0253   % do segments
0254   A  = 0.0;
0255   start = 1;
0256   nsegs = 0;
0257   while start + l - 1 <= nx
0258 
0259     % get segment
0260     xs  = x(start:start+l-1);
0261     
0262     % detrend segment
0263     switch order
0264       case -1
0265         % do nothing
0266       case 0
0267         xs = xs - mean(xs);
0268       case 1
0269         xs = detrend(xs);
0270       otherwise
0271         xs = polydetrend(xs, order);
0272     end
0273 
0274     % make DFT
0275     a   = C*xs;
0276     A   = A + a*conj(a);
0277     
0278     % make next step
0279     start = floor(start + l*Xolap);
0280     nsegs = nsegs + 1;
0281   end
0282 
0283   if mod(j,disp_each) == 0 || j == 1 || j == nf
0284     disp(sprintf('   averaged %d segments', nsegs));
0285   end
0286 
0287   A2ns    = 2.0*A/nsegs;
0288   S1      = win.ws;
0289   S12     = S1*S1;
0290   S2      = win.ws2;
0291   ENBW(j) = fs*S2/S12;
0292   Sxx(j)  = A2ns/fs/S2;
0293   S(j)    = A2ns/S12;
0294   
0295 end
0296 
0297 varargout{1} = S;
0298 varargout{2} = Sxx;
0299 varargout{3} = ENBW;
0300 
0301 
0302 
0303 %--------------------------------------------------------------------------
0304 % Get default params
0305 function plo = getDefaultPL()
0306 
0307 disp('* creating default plist...');
0308 plo = plist();
0309 plo = append(plo, param('Nfft', -1));
0310 plo = append(plo, param('Win', specwin('Kaiser', 10, 100)));
0311 plo = append(plo, param('Nolap', -1));
0312 plo = append(plo, param('Order', 0));
0313 disp('* done.');
0314 
0315 % END

Generated on Fri 02-Nov-2007 19:39:27 by m2html © 2003