Home > classes > @ao > lpsd.m

lpsd

PURPOSE ^

LPSD implement the LPSD algorithm for analysis objects.

SYNOPSIS ^

function varargout = lpsd(varargin)

DESCRIPTION ^

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

 DESCRIPTION: LPSD implement the LPSD algorithm for analysis objects.

 CALL:        bs = 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]

 M-FILE INFO: Get information about this methods by calling
              >> ao.getInfo('lpsd')

              Get information about a specified set-plist by calling:
              >> ao.getInfo('lpsd', 'None')

 VERSION:     $Id: lpsd.m,v 1.13 2008/09/05 14:14:58 hewitson Exp $

 HISTORY:     02-02-07 M Hewitson
                 Created

 References:  "Improved spectrum estimation from digitized time series
               on a logarithmic frequency axis", Michael Troebs, Gerhard Heinzel,
               Measurement 39 (2006) 120-129.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % LPSD implement the LPSD algorithm for analysis objects.
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 % DESCRIPTION: LPSD implement the LPSD algorithm for analysis objects.
0005 %
0006 % CALL:        bs = lpsd(as)
0007 %
0008 % INPUTS:      as  - array of analysis objects
0009 %              pl  - parameter list (see below)
0010 %
0011 % OUTPUTS:     bs  - array of analysis objects, one for each input
0012 %
0013 % PARAMETER LIST:
0014 %
0015 %     Kdes  - desired number of averages   (default 100)
0016 %     Lmin  - minimum segment length   (default 0)
0017 %     Jdes  - number of spectral frequencies to compute (default fs/2)
0018 %     win   - a specwin window object
0019 %             Only the design parameters of the window object are used; the
0020 %             window is recomputed for each DFT length inside the lpsd_core
0021 %             algorithm.
0022 %     Olap - desired overlap percentage (default: taken from window)
0023 %     Order - order of detrending
0024 %             -1 - no detrending
0025 %              0 - subtract mean
0026 %              1 - subtract linear fit
0027 %              N - subtract fit of polynomial, order N
0028 %     Scale - Scaling of output. Choose from:
0029 %               AS  - Amplitude (linear) Spectrum
0030 %               ASD - Amplitude (linear) Spectral Density
0031 %               PS  - Power Spectrum
0032 %               PSD - Power Spectral Density [default]
0033 %
0034 % M-FILE INFO: Get information about this methods by calling
0035 %              >> ao.getInfo('lpsd')
0036 %
0037 %              Get information about a specified set-plist by calling:
0038 %              >> ao.getInfo('lpsd', 'None')
0039 %
0040 % VERSION:     $Id: lpsd.m,v 1.13 2008/09/05 14:14:58 hewitson Exp $
0041 %
0042 % HISTORY:     02-02-07 M Hewitson
0043 %                 Created
0044 %
0045 % References:  "Improved spectrum estimation from digitized time series
0046 %               on a logarithmic frequency axis", Michael Troebs, Gerhard Heinzel,
0047 %               Measurement 39 (2006) 120-129.
0048 %
0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0050 
0051 function varargout = lpsd(varargin)
0052 
0053   import utils.const.*
0054 
0055   % Check if this is a call for parameters
0056   if utils.helper.isinfocall(varargin{:})
0057     varargout{1} = getInfo(varargin{3});
0058     return
0059   end
0060 
0061   utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
0062 
0063   % Collect input variable names
0064   in_names = cell(size(varargin));
0065   for ii = 1:nargin,in_names{ii} = inputname(ii);end
0066 
0067   % Collect all AOs and plists
0068   [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
0069   pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
0070 
0071   % Decide on a deep copy or a modify
0072   bs = copy(as, nargout);
0073 
0074   % Combine plists
0075   pl = combine(pl, getDefaultPlist);
0076 
0077   % Go through each input AO
0078   for j=1:numel(bs)
0079     % check this is a time-series object
0080     if ~isa(bs(j).data, 'tsdata')
0081       warning('!!! lpsd requires tsdata (time-series) inputs. Skipping AO %s', ao_invars{j});
0082     else
0083       % Desired number of averages
0084       Kdes = find(pl, 'Kdes');
0085       % num desired spectral frequencies
0086       Jdes = find(pl, 'Jdes');
0087       % Window function
0088       Win = find(pl, 'Win');
0089       % Overlap
0090       Nolap = find(pl, 'Olap');
0091       if Nolap < 0
0092         % use window function rov
0093         Nolap = Win.rov/100;
0094         utils.helper.msg(msg.PROC1, 'using recommended overlap of %d', Nolap);
0095         pl.pset('Olap', Nolap);
0096       end
0097       % Order of detrending
0098       Order = find(pl, 'Order');
0099       % Minimum segment length
0100       Lmin = find(pl, 'Lmin');
0101       % Get frequency vector
0102       [f, r, m, L, K] = ao.ltf_plan(length(bs(j).data.y), bs(j).data.fs, Nolap, 1, Lmin, Jdes, Kdes);
0103       % compute LPSD
0104       try
0105         if find(pl, 'M-FILE ONLY')
0106           % Using pure m-file version
0107           [P, Pxx, ENBW] = ao.mlpsd_m(bs(j).data.y, f, r, m, L, bs(j).data.fs, Win, Order, Nolap);
0108         else
0109           [P, Pxx, ENBW] = ao.mlpsd_mex(bs(j).data.y, f, r, m, L, bs(j).data.fs, Win, Order, Win.rov, Lmin);
0110         end
0111       catch
0112         warning('!!! mex file dft failed. Using m-file version of lpsd.');
0113         % Using pure m-file version
0114         [P, Pxx, ENBW] = ao.mlpsd_m(bs(j).data.y, f, r, m, L, bs(j).data.fs, Win, Order, Nolap);
0115       end
0116 
0117       % Keep the data shape of the input AO
0118       if size(bs(j).data.y,1) == 1
0119         P   = P.';
0120         Pxx = Pxx.';
0121         f   = f.';
0122       end
0123 
0124       % create new output fsdata
0125       scale = find(pl, 'Scale');
0126       switch scale
0127         case 'AS'
0128           fsd = fsdata(f, sqrt(P), bs(j).data.fs);
0129           fsd.setYunits(bs(j).data.yunits);
0130         case 'ASD'
0131           fsd = fsdata(f, sqrt(Pxx), bs(j).data.fs);
0132           fsd.setYunits(bs(j).data.yunits / unit('Hz^0.5'));
0133         case 'PS'
0134           fsd = fsdata(f, P, bs(j).data.fs);
0135           fsd.setYunits(bs(j).data.yunits.^2);
0136         case 'PSD'
0137           fsd = fsdata(f, Pxx, bs(j).data.fs);
0138           fsd.setYunits(bs(j).data.yunits.^2/unit('Hz'));
0139         otherwise
0140           error(['### Unknown scaling:' scale]);
0141       end
0142       fsd.setXunits('Hz');
0143       fsd.setEnbw(ENBW);
0144       % make output analysis object
0145       bs(j).data = fsd;
0146       % set name
0147       bs(j).setName(sprintf('l%s(%s)', scale, ao_invars{j}), 'internal');
0148       % Add history
0149       bs(j).addHistory(getInfo, pl, ao_invars(j), bs(j).hist);
0150       % Add processing info
0151       bs(j).procinfo = plist('r', r, 'm', m, 'l', L, 'k', K);
0152 
0153     end % End tsdata if/else
0154   end % End AO loop
0155 
0156   % Set output
0157   varargout{1} = bs;
0158 
0159 end
0160 
0161 %--------------------------------------------------------------------------
0162 % Get Info Object
0163 %--------------------------------------------------------------------------
0164 function ii = getInfo(varargin)
0165   if nargin == 1 && strcmpi(varargin{1}, 'None')
0166     sets = {};
0167     pl   = [];
0168   else
0169     sets = {'Default'};
0170     pl   = getDefaultPlist;
0171   end
0172   % Build info object
0173   ii = minfo(mfilename, 'ao', '', 'Signal Processing', '$Id: lpsd.m,v 1.13 2008/09/05 14:14:58 hewitson Exp $', sets, pl);
0174 end
0175 
0176 %--------------------------------------------------------------------------
0177 % Get Default Plist
0178 %--------------------------------------------------------------------------
0179 function pl = getDefaultPlist()
0180   pl = plist('Kdes', 100, 'Jdes', 1000, 'Lmin', 0,...
0181     'Win', getappdata(0, 'ltpda_default_spectral_window'), ...
0182     'Olap', -1, 'Order', 0, 'Scale', 'PSD');
0183 end
0184

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003