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
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