Home > classes > @ao > psd.m

psd

PURPOSE ^

PSD makes power spectral density estimates of the time-series objects

SYNOPSIS ^

function varargout = psd(varargin)

DESCRIPTION ^

 PSD makes power spectral density estimates of the time-series objects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: PSD makes power spectral density estimates of the
              time-series objects in the input analysis objects.

 CALL:        b = psd(a1,a2,a3,...,pl)

 INPUTS:      b    - output analysis objects
              aN   - input analysis objects
              pl   - input parameter list

              Makes PSD estimates of each input AO using the Welch Overlap method.

              If the last input argument is a parameter list (plist) it is used.
              The following parameters are recognised.

 PARAMETERS: 'Win'   - a specwin window object [default: Kaiser -200dB psll]
             'Olap'  - segment percent overlap [default: taken from window function]
             'Nfft'  - number of samples in each fft [default: length of input data]
                       A string value containing the variable 'fs' can
                       also be used, e.g., plist('Nfft', '2*fs')
             'Scale' - one of
                                'ASD' - amplitude spectral density
                                'PSD' - power spectral density [default]
                                'AS'  - amplitude spectrum
                                'PS'  - power spectrum
             'Order' - order of segment detrending
                        -1 - no detrending
                         0 - subtract mean [default]
                         1 - subtract linear fit
                         N - subtract fit of polynomial, order N

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

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

 VERSION:    $Id: psd.m,v 1.6 2008/09/05 11:05:29 ingo Exp $

 HISTORY:    07-02-2007 M Hewitson
                Creation

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % PSD makes power spectral density estimates of the time-series objects
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 % DESCRIPTION: PSD makes power spectral density estimates of the
0005 %              time-series objects in the input analysis objects.
0006 %
0007 % CALL:        b = psd(a1,a2,a3,...,pl)
0008 %
0009 % INPUTS:      b    - output analysis objects
0010 %              aN   - input analysis objects
0011 %              pl   - input parameter list
0012 %
0013 %              Makes PSD estimates of each input AO using the Welch Overlap method.
0014 %
0015 %              If the last input argument is a parameter list (plist) it is used.
0016 %              The following parameters are recognised.
0017 %
0018 % PARAMETERS: 'Win'   - a specwin window object [default: Kaiser -200dB psll]
0019 %             'Olap'  - segment percent overlap [default: taken from window function]
0020 %             'Nfft'  - number of samples in each fft [default: length of input data]
0021 %                       A string value containing the variable 'fs' can
0022 %                       also be used, e.g., plist('Nfft', '2*fs')
0023 %             'Scale' - one of
0024 %                                'ASD' - amplitude spectral density
0025 %                                'PSD' - power spectral density [default]
0026 %                                'AS'  - amplitude spectrum
0027 %                                'PS'  - power spectrum
0028 %             'Order' - order of segment detrending
0029 %                        -1 - no detrending
0030 %                         0 - subtract mean [default]
0031 %                         1 - subtract linear fit
0032 %                         N - subtract fit of polynomial, order N
0033 %
0034 % M-FILE INFO: Get information about this methods by calling
0035 %              >> ao.getInfo('psd')
0036 %
0037 %              Get information about a specified set-plist by calling:
0038 %              >> ao.getInfo('psd', 'None')
0039 %
0040 % VERSION:    $Id: psd.m,v 1.6 2008/09/05 11:05:29 ingo Exp $
0041 %
0042 % HISTORY:    07-02-2007 M Hewitson
0043 %                Creation
0044 %
0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0046 
0047 function varargout = psd(varargin)
0048 
0049   % Check if this is a call for parameters
0050   if utils.helper.isinfocall(varargin{:})
0051     varargout{1} = getInfo(varargin{3});
0052     return
0053   end
0054 
0055   import utils.const.*
0056   utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
0057   
0058   % Collect input variable names
0059   in_names = cell(size(varargin));
0060   for ii = 1:nargin,in_names{ii} = inputname(ii);end
0061 
0062   % Collect all AOs and plists
0063   [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
0064   pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
0065 
0066   % Decide on a deep copy or a modify
0067   bs = copy(as, nargout);
0068 
0069   % combine plists
0070   pl = combine(pl, getDefaultPlist());
0071 
0072   % Loop over input AOs
0073   for j=1:numel(bs)
0074     % check input data
0075     if isa(bs(j).data, 'tsdata')
0076       fs = bs(j).data.fs;
0077       %--------- Check all defaults are sensible
0078       % Check the number of points in FFT. If this is not set (<0) we set it
0079       % to be the length of the input data.
0080       Nfft = find(pl, 'Nfft');
0081       setWindow = 0;
0082       if ischar(Nfft)
0083         nNfft = floor(eval(Nfft));
0084         utils.helper.msg(msg.PROC1, 'setting Nfft to %s = %d', Nfft, nNfft);
0085         Nfft = nNfft;
0086       else
0087         if Nfft <= 0
0088           Nfft = length(bs(j).data.y);
0089           utils.helper.msg(msg.PROC1, 'using default Nfft of %g', Nfft);
0090           setWindow = 1;
0091         end
0092       end
0093       pl.pset('Nfft', Nfft);
0094 
0095       % Check the window function. If the length of the window doesn't match
0096       % NFFT then we resize it.
0097       Win = find(pl, 'Win');   % Window object
0098       if setWindow || length(Win.win)~=Nfft
0099         switch Win.type
0100           case {'Kaiser', 'Flattop'}
0101             Win = specwin(Win.type, Nfft, Win.psll);
0102           otherwise
0103             Win = specwin(Win.type, Nfft);
0104         end
0105         utils.helper.msg(msg.PROC1, 'reset window to %s(%d)', strrep(Win.type, '_', '\_'), length(Win.win));
0106       end
0107       pl.pset('Win', Win);
0108 
0109       % Check the overlap. If this is not set, we take the overlap from that
0110       % recommended by the window function.
0111       Olap = find(pl, 'Olap');
0112       if Olap <= 0
0113         Olap = Win.rov;
0114         utils.helper.msg(msg.PROC1, 'using default overlap of %2.1f%%', Olap);
0115       end
0116       pl.pset('Olap', Olap);
0117 
0118       % scaling of spectrum
0119       scale = find(pl, 'Scale');
0120 
0121       % check parameters
0122       if length(Win.win) ~= Nfft
0123         error('### Window length should be same as Nfft.');
0124       end
0125       % Compute PSD using pwelch
0126       [pxx, f, info] = ao.welch(bs(j), 'psd', pl);
0127       % Keep the data shape of the input AO
0128       if size(bs(j).data.y,1) == 1
0129         pxx = pxx.';
0130         f   = f.';
0131       end
0132       % create new output fsdata
0133       fsd = fsdata(f, pxx, fs);
0134       fsd.setXunits('Hz');
0135       fsd.setYunits(info.units);
0136       fsd.setEnbw(Win.nenbw*fs/Nfft);
0137       fsd.setNavs(info.navs);
0138       fsd.setT0(bs(j).data.t0);
0139       % update AO
0140       bs(j).data = fsd;
0141       % set name
0142       bs(j).setName(sprintf('%s(%s)', scale, ao_invars{j}), 'internal');
0143       % Add history
0144       bs(j).addHistory(getInfo, pl, ao_invars(j), bs(j).hist);
0145     else
0146       warning('### Ignoring input AO number %d (%s); it is not a time-series.', j, bs(j).name)
0147     end
0148   end % loop over analysis objects
0149 
0150   % set output
0151   varargout{1} = bs;
0152 end
0153 
0154 %--------------------------------------------------------------------------
0155 % Get Info Object
0156 %--------------------------------------------------------------------------
0157 function ii = getInfo(varargin)
0158   if nargin == 1 && strcmpi(varargin{1}, 'None')
0159     sets = {};
0160     pl   = [];
0161   else
0162     sets = {'Default'};
0163     pl   = getDefaultPlist;
0164   end
0165   % Build info object
0166   ii = minfo(mfilename, 'ao', '', utils.const.categories.sigproc, '$Id: psd.m,v 1.6 2008/09/05 11:05:29 ingo Exp $', sets, pl);
0167 end
0168 
0169 %--------------------------------------------------------------------------
0170 % Get Default Plist
0171 %--------------------------------------------------------------------------
0172 function pl = getDefaultPlist()
0173   pl = plist(...
0174     'Nfft', -1, ...
0175     'Win', getappdata(0, 'ltpda_default_spectral_window'), ...
0176     'Olap', -1, ...
0177     'Scale', 'PSD', ...
0178     'Order', 0);
0179 end
0180 % END
0181

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