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

Generated on Wed 04-Jul-2007 19:03:10 by m2html © 2003