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) 
 
   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)
     Lmin  - minimum segment length   (default 0)
     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.
     Olap - desired overlap percentage (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  
     Scale - Scaling of output. Choose from:
               AS  - Amplitude (linear) Spectrum
               ASD - Amplitude (linear) Spectral Density
               PS  - Power Spectrum
               PSD - Power Spectral Density [default]
 
 The following call returns a parameter list object that contains the
 default parameter values:
 
 >> pl = ltpda_lpsd('Params')
 
 The following call returns a string that contains the routine CVS version:

 >> version = ltpda_lpsd('Version')

 The following call returns a string that contains the routine category:

 >> category = ltpda_lpsd('Category')

 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.html,v 1.14 2008/03/31 10:27:39 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)
0008 %
0009 %   Inputs:
0010 %     as  - array of analysis objects
0011 %     pl  - parameter list (see below)
0012 %
0013 %   Outputs:
0014 %     bs  - array of analysis objects, one for each input
0015 %
0016 %   Parameter list:
0017 %     Kdes  - desired number of averages   (default 100)
0018 %     Lmin  - minimum segment length   (default 0)
0019 %     Jdes  - number of spectral frequencies to compute (default fs/2)
0020 %     win   - a specwin window object
0021 %             Only the design parameters of the window object are used; the
0022 %             window is recomputed for each DFT length inside the lpsd_core
0023 %             algorithm.
0024 %     Olap - desired overlap percentage (default: taken from window)
0025 %     Order - order of detrending
0026 %             -1 - no detrending
0027 %              0 - subtract mean
0028 %              1 - subtract linear fit
0029 %              N - subtract fit of polynomial, order N
0030 %     Scale - Scaling of output. Choose from:
0031 %               AS  - Amplitude (linear) Spectrum
0032 %               ASD - Amplitude (linear) Spectral Density
0033 %               PS  - Power Spectrum
0034 %               PSD - Power Spectral Density [default]
0035 %
0036 % The following call returns a parameter list object that contains the
0037 % default parameter values:
0038 %
0039 % >> pl = ltpda_lpsd('Params')
0040 %
0041 % The following call returns a string that contains the routine CVS version:
0042 %
0043 % >> version = ltpda_lpsd('Version')
0044 %
0045 % The following call returns a string that contains the routine category:
0046 %
0047 % >> category = ltpda_lpsd('Category')
0048 %
0049 % References:
0050 %     "Improved spectrum estimation from digitized time series
0051 %      on a logarithmic frequency axis", Michael Troebs, Gerhard Heinzel,
0052 %      Measurement 39 (2006) 120-129.
0053 %
0054 % M Hewitson 02-02-07
0055 %
0056 % $Id: ltpda_lpsd.html,v 1.14 2008/03/31 10:27:39 hewitson Exp $
0057 %
0058 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0059 
0060 ALGONAME = mfilename;
0061 VERSION  = '$Id: ltpda_lpsd.html,v 1.14 2008/03/31 10:27:39 hewitson Exp $';
0062 CATEGORY = 'Signal Processing';
0063 
0064 
0065 %% Check if this is a call for parameters, the CVS version string
0066 % or the function category
0067 if nargin == 1 && ischar(varargin{1})
0068   in = char(varargin{1});
0069   if strcmp(in, 'Params')
0070     varargout{1} = getDefaultPL();
0071     return
0072   elseif strcmp(in, 'Version')
0073     varargout{1} = VERSION;
0074     return
0075   elseif strcmp(in, 'Category')
0076     varargout{1} = CATEGORY;
0077     return
0078   end
0079 end
0080 
0081 % capture input variable names
0082 invars = {};
0083 as     = [];
0084 ps     = [];
0085 for j=1:nargin
0086   invars = [invars cellstr(inputname(j))];
0087   if isa(varargin{j}, 'ao')
0088     as = [as varargin{j}];
0089   end
0090   if isa(varargin{j}, 'plist')
0091     ps = [ps varargin{j}];
0092   end
0093 end
0094 
0095 
0096 
0097 % check plist
0098 if isempty(ps)
0099   pl = getDefaultPL();
0100 else
0101   pl = combine(ps, getDefaultPL);
0102 end
0103 
0104 % initialise output array
0105 bs  = []; 
0106 bxx = [];
0107 
0108 % Go through each input AO
0109 na = length(as);
0110 for i=1:na
0111   a = as(i);
0112   d = a.data;
0113   
0114   % check this is a time-series object
0115   if ~isa(d, 'tsdata')
0116     error('### lpsd requires tsdata (time-series) inputs.');
0117   end
0118   
0119   %--- check input parameters
0120   Kdes = find(pl, 'Kdes');   % desired averages
0121   if isempty(Kdes)
0122     Kdes = 100;
0123     disp(sprintf('! Using default Kdes of %d', Kdes))
0124     pl = append(pl, param('Kdes', Kdes));
0125   end
0126   
0127   % *** Only use for the old scheduler ***
0128   Kmin = find(pl, 'Kmin');    % minimum averages
0129   if isempty(Kmin)
0130     Kmin = 1;
0131     disp(sprintf('! Using default Kmin of %d', Kmin))
0132     pl = append(pl, param('Kmin', Kmin));
0133   end
0134   % *** Only use for the old scheduler ***
0135 
0136   Jdes = find(pl, 'Jdes');    % num desired spectral frequencies
0137   if isempty(Jdes)
0138     Jdes = round(d.fs/4);
0139     disp(sprintf('! Using default Jdes of %d', Jdes))
0140     pl = append(pl, param('Jdes', Jdes));
0141   end
0142   Win = find(pl, 'Win');
0143   if isempty(Win)
0144     Win = specwin('Kaiser', 100, 250);
0145     pl = append(pl, param('Win', Win));
0146   end  
0147   Nolap = find(pl, 'Nolap');    % desired overlap
0148   if isempty(Nolap) || Nolap == -1
0149     pl = remove(pl, 'Nolap'); % remove old value
0150     % use window function rov
0151     Nolap = Win.rov/100;
0152     disp(sprintf('! Using recommended Nolap of %d', Nolap))
0153     pl = append(pl, param('Nolap', Nolap));
0154   end
0155   Order = find(pl, 'Order');    % desired overlap
0156   if isempty(Order)
0157     Order = 0;
0158     disp(sprintf('! Using default detrending order of %d', Order))
0159     pl = append(pl, param('Order', Order));
0160   end
0161   Lmin = find(pl, 'Lmin');
0162   if isempty(Lmin)
0163     Lmin = 0;
0164   end
0165   Debug = find(pl, 'DEBUG');
0166   
0167   % Get frequency vector
0168   if find(pl, 'OLD SCHEDULER')
0169     [f,r,m,L,rr,rrr] = ltpda_compute_f(d.fs, length(d.y), Kdes, Kmin, Jdes, Nolap, Debug);
0170   else
0171     [f, r, m, L, K] = ltpda_ltf_plan(length(d.y), d.fs, Nolap, 1, Lmin, Jdes, Kdes);
0172   end
0173   
0174   % compute LPSD
0175   try
0176     if find(pl, 'M-FILE ONLY')
0177       % Using pure m-file version
0178       [P, Pxx, ENBW] = mlpsd(d.y, f, r, m, L, d.fs, Win, Order, Nolap);
0179     else
0180       [P, Pxx, ENBW] = mlpsd2(d.y, f, r, m, L, d.fs, Win, Order, Win.rov, Lmin);
0181     end
0182   catch
0183     % Using pure m-file version
0184     [P, Pxx, ENBW] = mlpsd(d.y, f, r, m, L, d.fs, Win, Order, Nolap);
0185   end
0186   
0187   % create new output fsdata
0188   scale = find(pl, 'Scale');
0189   switch scale
0190     case 'AS'
0191       fs = fsdata(f, sqrt(P), d.fs);
0192       fs = set(fs, 'yunits', [d.yunits]);      
0193     case 'ASD'
0194       fs = fsdata(f, sqrt(Pxx), d.fs);
0195       fs = set(fs, 'yunits', [d.yunits '/\surdHz']);
0196     case 'PS'
0197       fs = fsdata(f, P, d.fs);
0198       fs = set(fs, 'yunits', [d.yunits '^2']);      
0199     case 'PSD'
0200       fs = fsdata(f, Pxx, d.fs);
0201       fs = set(fs, 'yunits', [d.yunits '^2/Hz']);
0202     otherwise
0203       error(['### Unknown scaling:' scale]);
0204   end
0205   fs = set(fs, 'name', sprintf('lpsd(%s)', d.name));
0206   fs = set(fs, 'enbw', ENBW);
0207   
0208   % create new output history
0209   h = history(ALGONAME, VERSION, pl, a.hist);
0210   h = set(h, 'invars', invars);
0211 
0212   % make output analysis object
0213   b = ao(fs, h);
0214   
0215   % set name
0216   if isempty(invars{i})
0217     n = a.name;
0218   else
0219     n = invars{i};
0220   end
0221   b = setnh(b, 'name', sprintf('lpsd(%s)', n), ...
0222         'mfilename', ALGONAME);  
0223   bs = [bs b];
0224     
0225 end
0226 
0227 varargout{1} = bs;
0228 
0229 % include AS
0230 % if nargout > 1
0231 %   varargout{2} = bxx;
0232 % end
0233 
0234 
0235 
0236 %--------------------------------------------------------------------------
0237 %
0238 function varargout = mlpsd2(varargin)
0239 
0240 % MLPSD2
0241 %
0242 % [S,Sxx,ENBW] = mlpsd2(x,f,r,m,L,fs,win,order,olap)
0243 %
0244 % M Hewitson 19-05-07
0245 %
0246 %
0247 
0248 
0249 disp('%%%%%%%%   Using ltpda_dft.mex to compute core DFT ');
0250 
0251 % Get inputs
0252 x     = varargin{1};
0253 f     = varargin{2};
0254 r     = varargin{3};
0255 m     = varargin{4};
0256 L     = varargin{5};
0257 fs    = varargin{6};
0258 win   = varargin{7};
0259 order = varargin{8};
0260 olap  = varargin{9};
0261 Lmin  = varargin{10};
0262 
0263 twopi    = 2.0*pi;
0264 
0265 nx   = length(x);
0266 nf   = length(f);
0267 Sxx  = zeros(nf,1);
0268 S    = zeros(nf,1);
0269 ENBW = zeros(nf,1);
0270 
0271 disp_each = round(nf/100)*10;
0272 
0273 minReached = 0;
0274 
0275 for j=1:nf
0276   
0277   if mod(j,disp_each) == 0 || j == 1 || j == nf
0278     disp(sprintf('++ computing frequency %04d of %d: %f Hz', j, nf, f(j)));
0279   end
0280   
0281   % compute DFT exponent and window
0282   l = L(j);
0283   
0284   % Check if we need to update the window values
0285   % - once we reach Lmin, the window values don't change.
0286   if ~minReached
0287     switch win.name
0288       case 'Kaiser'
0289         win = specwin(win.name, l, win.psll);
0290       otherwise
0291         win = specwin(win.name, l);
0292     end
0293     if l==Lmin
0294       minReached = 1;
0295     end
0296   end
0297 
0298   p     = 1i * twopi * m(j)/l.*(0:l-1);
0299   C     = win.win .* exp(p);
0300   
0301   % Core DFT part in C-mex file
0302   [A, nsegs] = ltpda_dft(x, l, C, olap, order);
0303   
0304   if mod(j,disp_each) == 0 || j == 1 || j == nf
0305     disp(sprintf('   averaged %d segments', nsegs));
0306   end
0307 
0308   A2ns    = 2.0*A/nsegs;
0309   S1      = win.ws;
0310   S12     = S1*S1;
0311   S2      = win.ws2;
0312   ENBW(j) = fs*S2/S12;
0313   Sxx(j)  = A2ns/fs/S2;
0314   S(j)    = A2ns/S12;
0315   
0316 end
0317 
0318 varargout{1} = S;
0319 varargout{2} = Sxx;
0320 varargout{3} = ENBW;
0321 
0322 
0323 %--------------------------------------------------------------------------
0324 %
0325 function varargout = mlpsd(varargin)
0326 
0327 % MLPSD
0328 %
0329 % [S,Sxx,ENBW] = mlpsd(x,f,r,m,L,fs,win,order,olap)
0330 %
0331 % M Hewitson 19-05-07
0332 %
0333 %
0334 
0335 % Get inputs
0336 x     = varargin{1};
0337 f     = varargin{2};
0338 r     = varargin{3};
0339 m     = varargin{4};
0340 L     = varargin{5};
0341 fs    = varargin{6};
0342 win   = varargin{7};
0343 order = varargin{8};
0344 olap  = varargin{9};
0345 
0346 twopi    = 2.0*pi;
0347 
0348 nx   = length(x);
0349 nf   = length(f);
0350 Sxx  = zeros(nf,1);
0351 S    = zeros(nf,1);
0352 ENBW = zeros(nf,1);
0353 
0354 disp_each = round(nf/100)*10;
0355 
0356 for j=1:nf
0357   
0358   if mod(j,disp_each) == 0 || j == 1 || j == nf
0359     disp(sprintf('++ computing frequency %d of %d: %f Hz', j, nf, f(j)));
0360   end
0361 
0362   % compute DFT exponent and window
0363   l = L(j);
0364   switch win.name
0365     case 'Kaiser'
0366       win = specwin(win.name, l, win.psll);
0367     otherwise
0368       win = specwin(win.name, l);
0369   end
0370 
0371   p     = 1i * twopi * m(j)/l.*[0:l-1];
0372   C     = win.win .* exp(p);
0373   Xolap = (1-olap);
0374   % do segments
0375   A  = 0.0;
0376   start = 1;
0377   nsegs = 0;
0378   while start + l - 1 <= nx
0379 
0380     % get segment
0381     xs  = x(start:start+l-1);
0382     
0383     % detrend segment
0384     switch order
0385       case -1
0386         % do nothing
0387       case 0
0388         xs = xs - mean(xs);
0389       case 1
0390         xs = detrend(xs);
0391       otherwise
0392         xs = polydetrend(xs, order);
0393     end
0394 
0395     % make DFT
0396     a   = C*xs;
0397     A   = A + a*conj(a);
0398     
0399     % make next step
0400     start = floor(start + l*Xolap);
0401     nsegs = nsegs + 1;
0402   end
0403 
0404   if mod(j,disp_each) == 0 || j == 1 || j == nf
0405     disp(sprintf('   averaged %d segments', nsegs));
0406   end
0407 
0408   A2ns    = 2.0*A/nsegs;
0409   S1      = win.ws;
0410   S12     = S1*S1;
0411   S2      = win.ws2;
0412   ENBW(j) = fs*S2/S12;
0413   Sxx(j)  = A2ns/fs/S2;
0414   S(j)    = A2ns/S12;
0415   
0416 end
0417 
0418 varargout{1} = S;
0419 varargout{2} = Sxx;
0420 varargout{3} = ENBW;
0421 
0422 
0423 
0424 %--------------------------------------------------------------------------
0425 % Get default params
0426 function plo = getDefaultPL()
0427 
0428 disp('* creating default plist...');
0429 plo = plist('Kdes', 100, ...
0430   'Kmin', 1, ...
0431   'Lmin', 0, ...
0432   'Jdes', 100, ...
0433   'Win', specwin('Kaiser', 10, 200), ...
0434   'Olap', -1, ...
0435   'Order', -1, ...
0436   'Scale', 'PSD');
0437 
0438 disp('* done.');
0439 
0440 % END

Generated on Mon 31-Mar-2008 12:20:24 by m2html © 2003