Home > m > sigproc > frequency_domain > ltpda_pwelch.m

ltpda_pwelch

PURPOSE ^

LTPDA_PWELCH makes power spectral density estimates of the time-series

SYNOPSIS ^

function varargout = ltpda_pwelch(varargin)

DESCRIPTION ^

 LTPDA_PWELCH makes power spectral density estimates of the time-series
 objects in the input analysis objects.
 
 Usage:  b = ltpda_pwelch(a1,a2,a3,...,pl)
 
   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)
 
 M Hewitson 07-02-07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Fri 08-Jun-2007 16:09:11 by m2html © 2003