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]
             'Nolap' - segment overlap (number of samples) [default: taken from window function]
             'Nfft'  - number of samples in each fft [default: fs of input data]
             'Debug' - debug level for terminal output (0-1)

 VERSION:    $Id: ltpda_pwelch.m,v 1.12 2007/07/16 12:52:21 ingo 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')

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

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 %             'Nolap' - segment overlap (number of samples) [default: taken from window function]
0022 %             'Nfft'  - number of samples in each fft [default: fs of input data]
0023 %             'Debug' - debug level for terminal output (0-1)
0024 %
0025 % VERSION:    $Id: ltpda_pwelch.m,v 1.12 2007/07/16 12:52:21 ingo Exp $
0026 %
0027 % HISTORY:    07-02-2007 M Hewitson
0028 %                Creation
0029 %
0030 % The following call returns a parameter list object that contains the
0031 % default parameter values:
0032 %
0033 % >> pl = ltpda_pwelch('Params')
0034 %
0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0036 
0037 ALGONAME = mfilename;
0038 VERSION  = '$Id: ltpda_pwelch.m,v 1.12 2007/07/16 12:52:21 ingo Exp $';
0039 
0040 % Initialise output
0041 bs = [];
0042 
0043 % Check if this is a call for parameters
0044 if nargin == 1
0045   in = char(varargin{1});
0046   if strcmp(in, 'Params')
0047     varargout{1} = getDefaultPL();
0048     return
0049   end
0050 end
0051 
0052 % Read inputs
0053 invars = {};
0054 as   = [];
0055 pls  = [];
0056 for j=1:nargin
0057   if isa(varargin{j}, 'ao')
0058     as = [as varargin{j}];
0059     invars = [invars cellstr(inputname(j))];
0060   end
0061   if isa(varargin{j}, 'plist')
0062     pls = [pls varargin{j}];
0063   end
0064 end
0065 if isempty(pls)
0066   pls = plist();
0067 end
0068 na = length(as);
0069 
0070 % Loop over input AOs
0071 for j=1:na
0072 
0073   % Unpack AO
0074   a = as(j);
0075   d = a.data;
0076 
0077   if isa(a, 'ao')
0078     % check input data
0079     if isa(d, 'tsdata')
0080 
0081       % combine plists
0082       plo = combine(pls, getDefaultPL());
0083 
0084       % check all defaults are sensible
0085       Nfft = find(plo, 'Nfft');   % Number of points in each fft
0086       setWindow = 0;
0087       if Nfft <= 0
0088         Nfft = length(d.x);
0089         disp(sprintf('! Using default Nfft of %g', Nfft))
0090         setWindow = 1;
0091       end
0092       plo = set(plo, 'Nfft', Nfft);
0093 
0094       Win = find(plo, 'Win');   % Window object
0095       if setWindow || length(Win.win)~=Nfft
0096         switch Win.name
0097           case {'Kaiser', 'Flattop'}
0098             Win = specwin(Win.name, Nfft, Win.psll);
0099           otherwise
0100             Win = specwin(Win.name, Nfft);
0101         end
0102         disp(sprintf('! Reset window to %s(%d)', strrep(Win.name, '_', '\_'), length(Win.win)))
0103       end
0104       plo = set(plo, 'Win', Win);
0105 
0106       Nolap = find(plo, 'Nolap');   % Amount to overlap each fft
0107       if Nolap <= 0
0108         Nolap = round(Win.rov*Nfft/100);
0109         disp(sprintf('! Using default overlap of %d samples', Nolap))
0110       end
0111       plo = set(plo, 'Nolap', Nolap);
0112 
0113       % check parameters
0114       if length(Win.win) ~= Nfft
0115         error('### Window length should be same as Nfft.');
0116       end
0117 
0118       % Compute PSD using pwelch
0119       [f, pxx, info] = ltpda_psd(d.x, Win.win, Nolap, Nfft, d.fs, 'ASD', d.yunits);
0120 
0121       % create new output fsdata
0122       fs = fsdata(f, pxx, d.fs);
0123       fs = set(fs, 'name', sprintf('PSD(%s)', d.name));
0124       fs = set(fs, 'yunits', info.units);
0125       fs = set(fs, 'enbw', Win.nenbw*d.fs/Nfft);
0126       fs = set(fs, 'navs', info.navs);
0127 
0128       % create new output history
0129       h = history(ALGONAME, VERSION, plo, a.hist);
0130       h = set(h, 'invars', invars);
0131 
0132       % make output analysis object
0133       b = ao(fs, h);
0134 
0135       % set name
0136       if na~=length(invars) || isempty(invars{j})
0137         nj = a.name;
0138       else
0139         nj = invars{j};
0140       end
0141       b = set(b, 'name', sprintf('PSD(%s)', nj));
0142 
0143       % add to outputs
0144       bs = [bs b];
0145 
0146     else
0147       warning('### Ignoring input AO number %d (%s); it is not a time-series.', j, a.name)
0148     end
0149   end
0150 end % loop over analysis objects
0151 
0152 % set output
0153 varargout{1} = bs;
0154 
0155 
0156 %--------------------------------------------------------------------------
0157 % Get default params
0158 function plo = getDefaultPL()
0159 
0160 disp('* creating default plist...');
0161 plo = plist();
0162 plo = append(plo, param('Nfft', -1));
0163 plo = append(plo, param('Win', specwin('Kaiser', 1, 200)));
0164 plo = append(plo, param('Nolap', -1));
0165 disp('* done.');
0166 
0167 
0168

Generated on Mon 03-Sep-2007 12:12:34 by m2html © 2003