Home > m > sigproc > frequency_domain > ltpda_cpsd.m

ltpda_cpsd

PURPOSE ^

LTPDA_CPSD makes cross-spectral density estimates of the time-series

SYNOPSIS ^

function varargout = ltpda_cpsd(varargin)

DESCRIPTION ^

 LTPDA_CPSD makes cross-spectral density estimates of the time-series
 objects in the input analysis objects. CPSDs are computed
 using MATLAB's cpsd (>> help cpsd).
 
 
 Usage:  b = ltpda_cpsd(a1,a2,a3,...,pl)
 
   b    - output analysis objects
   aN   - input analysis objects (at least two)
   pl   - input parameter list
 
 The function makes CPSD estimates between a1 and all other input AOs. 
 Therefore, if the input argument list contains N analysis objects, the 
 output a will contain N-1 CPSD estimates.
 
 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_cpsd(varargin)
0002 
0003 % LTPDA_CPSD makes cross-spectral density estimates of the time-series
0004 % objects in the input analysis objects. CPSDs are computed
0005 % using MATLAB's cpsd (>> help cpsd).
0006 %
0007 %
0008 % Usage:  b = ltpda_cpsd(a1,a2,a3,...,pl)
0009 %
0010 %   b    - output analysis objects
0011 %   aN   - input analysis objects (at least two)
0012 %   pl   - input parameter list
0013 %
0014 % The function makes CPSD estimates between a1 and all other input AOs.
0015 % Therefore, if the input argument list contains N analysis objects, the
0016 % output a will contain N-1 CPSD estimates.
0017 %
0018 % If the last input argument is a parameter list (plist) it is used. The
0019 % following parameters are recognised.
0020 %
0021 % Parameters:
0022 %   'Win'   - a specwin window object [default: Kaiser -200dB psll]
0023 %   'Nolap' - segment overlap (number of samples) [default: taken from window function]
0024 %   'Nfft'  - number of samples in each fft [default: fs of input data]
0025 %   'Debug' - debug level for terminal output (0-1)
0026 %
0027 % M Hewitson 07-02-07
0028 %
0029 
0030 ALGONAME = mfilename;
0031 VERSION  = '$Id: ltpda_cpsd.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $';
0032 
0033 
0034 % Check if this is a call for parameters
0035 if nargin == 1
0036   in = char(varargin{1});
0037   if strcmp(in, 'Params')
0038     varargout{1} = getDefaultPL();
0039     return
0040   end
0041 end
0042 
0043 % capture input variable names
0044 invars = {};
0045 as     = [];
0046 ps     = [];
0047 for j=1:nargin
0048   invars = [invars cellstr(inputname(j))];
0049   if isa(varargin{j}, 'ao')
0050     as = [as varargin{j}];
0051   end
0052   if isa(varargin{j}, 'plist')
0053     ps = [ps varargin{j}];
0054   end
0055 end
0056 
0057 
0058 % check plist
0059 if isempty(ps)
0060   pl = getDefaultPL();
0061 else
0062   pl = combine(ps, getDefaultPL);
0063 end
0064 
0065 varargout{1} = ltpda_xspec(as, pl, 'CPSD', ALGONAME, VERSION, invars);
0066 
0067 
0068 %--------------------------------------------------------------------------
0069 % Get default params
0070 function plo = getDefaultPL()
0071 
0072 disp('* creating default plist...');
0073 plo = plist();
0074 plo = append(plo, param('Nfft', -0.5));
0075 plo = append(plo, param('Win', specwin('Kaiser', 1, 100)));
0076 plo = append(plo, param('Nolap', -1));
0077 disp('* done.');
0078 
0079 
0080 
0081 %
0082 % % Initialise output
0083 % bs = [];
0084 % if nargin < 2
0085 %   error('### LTPDA_CPSD needs at least two AOs to make a transfer function.');
0086 % end
0087 %
0088 % % Read inputs
0089 % invars = {};
0090 % as   = [];
0091 % pl   = [];
0092 % for j=1:nargin
0093 %   if isa(varargin{j}, 'ao')
0094 %     as = [as varargin{j}];
0095 %     invars = [invars cellstr(inputname(j))];
0096 %   end
0097 %   if isa(varargin{j}, 'plist')
0098 %     pl = varargin{j};
0099 %   end
0100 % end
0101 % if isempty(pl)
0102 %   pl = plist();
0103 % end
0104 % na = length(as);
0105 %
0106 % % Output AO
0107 % oa = as(1);
0108 % od = oa.data;
0109 % if ~isa(od, 'tsdata')
0110 %   error('### the first analysis object must contain time-series data.');
0111 % end
0112 %
0113 % % name for this object
0114 % if isempty(invars{1})
0115 %   n1 = oa.name;
0116 % else
0117 %   n1 = invars{1};
0118 % end
0119 %
0120 % % Loop over input AOs
0121 % for j=2:na
0122 %
0123 %   % Unpack AO
0124 %   a = as(j);
0125 %   d = a.data;
0126 %
0127 %   % check input data
0128 %   if isa(d, 'tsdata')
0129 %
0130 %     plo = plist();
0131 %
0132 %     % Unpack inputs
0133 %     Nfft = find(pl, 'Nfft');   % Number of points in each fft
0134 %     if isempty(Nfft)
0135 %       Nfft = length(d.x);
0136 %       disp(sprintf('! Using default sample rate of %g', Nfft))
0137 %     end
0138 %     plo = append(plo, param('Nfft', Nfft));
0139 %
0140 %     Win = find(pl, 'Win');   % Window object
0141 %     if isempty(Win)
0142 %       Win = specwin('Kaiser', Nfft, 200);
0143 %       disp(sprintf('! Using default Window of %s', strrep(Win.name, '_', '\_')))
0144 %     end
0145 %     plo = append(plo, param('Win', Win));
0146 %
0147 %     Nolap = find(pl, 'Nolap');   % Amount to overlap each fft
0148 %     if isempty(Nolap)
0149 %       Nolap = round(Win.rov*Nfft/100);
0150 %       disp(sprintf('! Using default overlap of %d samples', Nolap))
0151 %     end
0152 %     plo = append(plo, param('Nolap', Nolap));
0153 %
0154 %     Debug = find(pl, 'Debug');   % Debug level
0155 %     if isempty(Debug)
0156 %       Debug = 0;
0157 %       disp(sprintf('! Using default debug level of %d', Debug))
0158 %     end
0159 %     plo = append(plo, param('Debug', Debug));
0160 %
0161 %     % display input information
0162 %     switch Debug
0163 %       case 0
0164 %
0165 %       case 1
0166 %         disp(char(plo));
0167 %       otherwise
0168 %         error('### Unknown debug level');
0169 %     end
0170 %
0171 %     % Check this input data makes sense
0172 %     %
0173 %     % Note: we should resample here. That's something for later.
0174 %     %
0175 %     if od.fs ~= d.fs
0176 %       error('### all analysis objects must contain time-series with the same sample rate for now.');
0177 %     end
0178 %
0179 %     % Compute TF using tfestimate
0180 %     [pxy,f] = cpsd(d.x,od.x,Win.win,Nolap,Nfft,d.fs);
0181 %
0182 %     % create new output fsdata
0183 %     fs = fsdata(f, pxy, d.fs);
0184 %     fs = set(fs, 'name', sprintf('CPSD(%s,%s)', d.name, od.name));
0185 %     fs = set(fs, 'enbw', Win.nenbw);
0186 %
0187 %     % create new output history
0188 %     h = history(ALGONAME, VERSION, plo, [oa.hist a.hist]);
0189 %     h = set(h, 'invars', invars);
0190 %
0191 %     % make output analysis object
0192 %     b = ao(fs, h);
0193 %
0194 %     % set name
0195 %     if isempty(invars{j})
0196 %       nj = a.name;
0197 %     else
0198 %       nj = invars{j};
0199 %     end
0200 %     b = set(b, 'name', sprintf('CPSD(%s,%s)', n1, nj));
0201 %
0202 %     % add to outputs
0203 %     bs = [bs b];
0204 %
0205 %   else
0206 %     warning('### Ignoring input AO number %d (%s); it is not a time-series.', j, a.name)
0207 %   end
0208 % end % loop over analysis objects
0209 %
0210

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