Home > m > sigproc > frequency_domain > ltpda_pwelch.m

ltpda_pwelch

PURPOSE ^

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

SYNOPSIS ^

function varargout = ltpda_pwelch(varargin)

DESCRIPTION ^

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

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

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

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

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

             Makes PSD estimates using ltpda_psd() of each input AO.

             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: fs of input data]
             '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  

 VERSION:    $Id: ltpda_pwelch.m,v 1.18 2008/02/20 07:39:01 hewitson Exp $

 HISTORY:    07-02-2007 M Hewitson
                Creation

 The following call returns a parameter list object that contains the
 default parameter values:

 >> pl = ltpda_pwelch('Params')
 
 The following call returns a string that contains the routine CVS version:

 >> version = ltpda_pwelch(ao,'Version')

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

 >> category = ltpda_pwelch(ao,'Category')

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_pwelch(varargin)
0002 % LTPDA_PWELCH makes power spectral density estimates of the time-series objects
0003 %
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % DESCRIPTION: LTPDA_PWELCH makes power spectral density estimates of the
0007 %              time-series objects in the input analysis objects.
0008 %
0009 % CALL:       b = ltpda_pwelch(a1,a2,a3,...,pl)
0010 %
0011 % INPUTS:     b    - output analysis objects
0012 %             aN   - input analysis objects
0013 %             pl   - input parameter list
0014 %
0015 %             Makes PSD estimates using ltpda_psd() of each input AO.
0016 %
0017 %             If the last input argument is a parameter list (plist) it is used.
0018 %             The following parameters are recognised.
0019 %
0020 % PARAMETERS: 'Win'   - a specwin window object [default: Kaiser -200dB psll]
0021 %             'olap' - segment percent overlap [default: taken from window function]
0022 %             'Nfft'  - number of samples in each fft [default: fs of input data]
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 % VERSION:    $Id: ltpda_pwelch.m,v 1.18 2008/02/20 07:39:01 hewitson Exp $
0035 %
0036 % HISTORY:    07-02-2007 M Hewitson
0037 %                Creation
0038 %
0039 % The following call returns a parameter list object that contains the
0040 % default parameter values:
0041 %
0042 % >> pl = ltpda_pwelch('Params')
0043 %
0044 % The following call returns a string that contains the routine CVS version:
0045 %
0046 % >> version = ltpda_pwelch(ao,'Version')
0047 %
0048 % The following call returns a string that contains the routine category:
0049 %
0050 % >> category = ltpda_pwelch(ao,'Category')
0051 %
0052 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0053 
0054 ALGONAME = mfilename;
0055 VERSION  = '$Id: ltpda_pwelch.m,v 1.18 2008/02/20 07:39:01 hewitson Exp $';
0056 CATEGORY = 'Signal Processing';
0057 
0058 % Initialise output
0059 bs = [];
0060 
0061 % Check if this is a call for parameters
0062 if nargin == 1
0063   in = char(varargin{1});
0064   if strcmpi(in, 'Params')
0065     varargout{1} = getDefaultPL();
0066     return
0067   elseif strcmpi(in, 'Version')
0068     varargout{1} = VERSION;
0069     return
0070   elseif strcmpi(in, 'Category')
0071     varargout{1} = CATEGORY;
0072     return
0073   end
0074 end
0075 
0076 % Collect input ao's, plist's and ao variable names
0077 in_names = {};
0078 for ii = 1:nargin
0079   in_names{end+1} = inputname(ii);
0080 end
0081 [as, pls, invars] = collect_inputs(varargin, in_names);
0082 na = length(as);
0083 
0084 % Loop over input AOs
0085 for j=1:na
0086 
0087   % Unpack AO
0088   a = as(j);
0089   d = a.data;
0090 
0091   if isa(a, 'ao')
0092     % check input data
0093     if isa(d, 'tsdata')
0094 
0095       % combine plists
0096       plo = combine(pls, getDefaultPL());
0097 
0098       % check all defaults are sensible
0099       Nfft = find(plo, 'Nfft');   % Number of points in each fft
0100       setWindow = 0;
0101       if Nfft <= 0
0102         Nfft = length(d.y);
0103         disp(sprintf('! Using default Nfft of %g', Nfft))
0104         setWindow = 1;
0105       end
0106       plo = pset(plo, 'Nfft', Nfft);
0107 
0108       Win = find(plo, 'Win');   % Window object
0109       if setWindow || length(Win.win)~=Nfft
0110         switch Win.name
0111           case {'Kaiser', 'Flattop'}
0112             Win = specwin(Win.name, Nfft, Win.psll);
0113           otherwise
0114             Win = specwin(Win.name, Nfft);
0115         end
0116         disp(sprintf('! Reset window to %s(%d)', strrep(Win.name, '_', '\_'), length(Win.win)))
0117       end
0118       plo = pset(plo, 'Win', Win);
0119 
0120       olap = find(plo, 'olap');   % Amount to overlap each fft
0121       if olap <= 0
0122         olap = Win.rov;
0123         disp(sprintf('! Using default overlap of %d /%', olap))
0124       end
0125       plo = pset(plo, 'olap', olap);
0126 
0127       % Get detrending order
0128       order = find(plo, 'Order');
0129       
0130       % scaling of spectrum
0131       scale = find(plo, 'Scale');
0132       
0133       % check parameters
0134       if length(Win.win) ~= Nfft
0135         error('### Window length should be same as Nfft.');
0136       end
0137 
0138       % Compute PSD using pwelch
0139       [f, pxx, info] = ltpda_psd(d.y, Win.win, round(olap/100*Nfft), Nfft, d.fs, scale, d.yunits, order);
0140       
0141       % create new output fsdata
0142       fs = fsdata(f, pxx, d.fs);
0143       fs = set(fs, 'name', sprintf('%s(%s)', scale, d.name));
0144       fs = set(fs, 'yunits', info.units);
0145       fs = set(fs, 'enbw', Win.nenbw*d.fs/Nfft);
0146       fs = set(fs, 'navs', info.navs);
0147 
0148       % create new output history
0149       h = history(ALGONAME, VERSION, plo, a.hist);
0150       h = set(h, 'invars', invars);
0151 
0152       % make output analysis object
0153       b = ao(fs, h);
0154 
0155       % set name
0156       if na~=length(invars) || isempty(invars{j})
0157         nj = a.name;
0158       else
0159         nj = invars{j};
0160       end
0161       b = setnh(b, 'name', sprintf('%s(%s)', scale, nj));
0162 
0163       % add to outputs
0164       bs = [bs b];
0165 
0166     else
0167       warning('### Ignoring input AO number %d (%s); it is not a time-series.', j, a.name)
0168     end
0169   end
0170 end % loop over analysis objects
0171 
0172 % set output
0173 varargout{1} = bs;
0174 
0175 
0176 %--------------------------------------------------------------------------
0177 % Get default params
0178 function plo = getDefaultPL()
0179 
0180 plo = plist();
0181 plo = append(plo, param('Nfft', -1));
0182 plo = append(plo, param('Win', specwin('Kaiser', 1, 200)));
0183 plo = append(plo, param('Olap', -1));
0184 plo = append(plo, param('Scale', 'PSD'));
0185 plo = append(plo, param('Order', 0));
0186 
0187 
0188

Generated on Fri 07-Mar-2008 15:46:43 by m2html © 2003