Home > m > sigproc > frequency_domain > ltpda_lcohere.m

ltpda_lcohere

PURPOSE ^

LTPDA_LCOHERE implement weighted coherence estimation computed on a log frequency axis.

SYNOPSIS ^

function varargout = ltpda_lcohere(varargin)

DESCRIPTION ^

 LTPDA_LCOHERE implement weighted coherence estimation computed on a log frequency axis.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 LTPDA_LCOHERE implement weighted cross-power estimation computed on a
 log frequency axis.

 Usage:
   >> txy = ltpda_lcohere(as)
   >> txy = ltpda_lcohere(as, pl)
   >> [txy, sig2] = ltpda_lcohere(as, pl)  % returns variance estimates

   Inputs:
     as  - array of analysis objects
     pl  - parameter list (see below)

   Outputs:
     bs  - array of analysis objects, one for each input

   Parameter list:
     Kdes  - desired number of averages   [default: 100]
     Lmin  - minimum segment length   [default: 0]
     Jdes  - number of spectral frequencies to compute [default: 1000]
     Win   - a specwin window object
             Only the design parameters of the window object are used; the
             window is recomputed for each DFT length inside the lpsd_core
             algorithm. [default: Kaiser with -200dB PSLL]
     Nolap - desired overlap [default: taken from window]
     Order - order of detrending.
             -1 - no detrending
              0 - subtract mean
              1 - subtract linear fit
              N - subtract fit of polynomial, order N

 The following call returns a parameter list object that contains the
 default parameter values:

 >> pl = ltpda_lcohere('Params')


 M Hewitson 02-02-07

 $Id: ltpda_lcohere.html,v 1.12 2008/03/31 10:27:39 hewitson Exp $

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_lcohere(varargin)
0002 % LTPDA_LCOHERE implement weighted coherence estimation computed on a log frequency axis.
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % LTPDA_LCOHERE implement weighted cross-power estimation computed on a
0006 % log frequency axis.
0007 %
0008 % Usage:
0009 %   >> txy = ltpda_lcohere(as)
0010 %   >> txy = ltpda_lcohere(as, pl)
0011 %   >> [txy, sig2] = ltpda_lcohere(as, pl)  % returns variance estimates
0012 %
0013 %   Inputs:
0014 %     as  - array of analysis objects
0015 %     pl  - parameter list (see below)
0016 %
0017 %   Outputs:
0018 %     bs  - array of analysis objects, one for each input
0019 %
0020 %   Parameter list:
0021 %     Kdes  - desired number of averages   [default: 100]
0022 %     Lmin  - minimum segment length   [default: 0]
0023 %     Jdes  - number of spectral frequencies to compute [default: 1000]
0024 %     Win   - a specwin window object
0025 %             Only the design parameters of the window object are used; the
0026 %             window is recomputed for each DFT length inside the lpsd_core
0027 %             algorithm. [default: Kaiser with -200dB PSLL]
0028 %     Nolap - desired overlap [default: taken from window]
0029 %     Order - order of detrending.
0030 %             -1 - no detrending
0031 %              0 - subtract mean
0032 %              1 - subtract linear fit
0033 %              N - subtract fit of polynomial, order N
0034 %
0035 % The following call returns a parameter list object that contains the
0036 % default parameter values:
0037 %
0038 % >> pl = ltpda_lcohere('Params')
0039 %
0040 %
0041 % M Hewitson 02-02-07
0042 %
0043 % $Id: ltpda_lcohere.html,v 1.12 2008/03/31 10:27:39 hewitson Exp $
0044 %
0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0046 
0047 ALGONAME = mfilename;
0048 VERSION  = '$Id: ltpda_lcohere.html,v 1.12 2008/03/31 10:27:39 hewitson Exp $';
0049 CATEGORY = 'Signal Processing';
0050 
0051 % Check if this is a call for parameters
0052 if nargin == 1
0053   in = char(varargin{1});
0054   if strcmp(in, 'Params')
0055     varargout{1} = getDefaultPL();
0056     return
0057   elseif strcmp(in, 'Version')
0058     varargout{1} = VERSION;
0059     return
0060   elseif strcmp(in, 'Category')
0061     varargout{1} = CATEGORY;
0062     return
0063   end
0064 end
0065 
0066 % capture input variable names
0067 % Collect input ao's, plist's and ao variable names
0068 in_names = {};
0069 for ii = 1:nargin
0070   in_names{end+1} = inputname(ii);
0071 end
0072 [as, ps, invars] = collect_inputs(varargin, in_names);
0073 
0074 % initialise output array
0075 bs = [];
0076 
0077 % check plist
0078 if isempty(ps)
0079   pl = getDefaultPL();
0080 else
0081   pl = combine(ps, getDefaultPL);
0082 end
0083 
0084 % Go through each input AO and build a matrix of signals
0085 na    = length(as);
0086 
0087 %----------------- Resample all AOs
0088 
0089 disp('*** Resampling AOs...');
0090 fsmax = findFsMax(as);
0091 fspl  = plist(param('fsout', fsmax));
0092 bs = [];
0093 for i=1:na
0094   a = as(i);
0095   % check this is a time-series object
0096   if ~isa(a.data, 'tsdata')
0097     error('### ltpda_lcohere requires tsdata (time-series) inputs.');
0098   end
0099   % Check Fs
0100   if a.data.fs ~= fsmax
0101     warning(sprintf('!!! Resampling AO %s/%s to %f Hz', a.name, a.data.name, fsmax));
0102     a = resample(a, fspl);
0103   end
0104   bs = [bs a];
0105 end
0106 as = bs;
0107 na = length(as);
0108 clear bs;
0109 disp('*** Done.');
0110 
0111 %----------------- Truncate all vectors
0112 
0113 % Get shortest vector
0114 disp('*** Truncating all vectors...');
0115 lmin = findShortestVector(as);
0116 nsecs = lmin / fsmax;
0117 bs = [];
0118 for i=1:na
0119   a = as(i);
0120   if len(a) ~= lmin
0121     warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0122     bs = [bs  select(a, 1:lmin)];
0123   else
0124     bs = [bs a];
0125   end
0126 end
0127 as = bs;
0128 clear bs;
0129 disp('*** Done.');
0130 
0131 
0132 %----------------- Build signal Matrix
0133 N     = len(as(1)); % length of first signal
0134 iS    = zeros(na, N);
0135 for j=1:na
0136   a = as(j);
0137   iS(j,:) = a.data.y;
0138 end
0139 
0140 %----------------- check input parameters
0141 
0142 % Averaging
0143 Kdes = find(pl, 'Kdes');   % desired averages
0144 Kmin = find(pl, 'Kmin');   % minimum averages
0145 Jdes = find(pl, 'Jdes');   % num desired spectral frequencies
0146 Win   = find(pl, 'Win');   % Window object
0147 Nolap = find(pl, 'Nolap');    % desired overlap
0148 if Nolap == -1
0149   % use window function rov
0150   Nolap = Win.rov/100;
0151   disp(sprintf('! Using recommended overlap of %d', Nolap))
0152   pl = pset(pl, 'Nolap', Nolap);
0153 end
0154 Order = find(pl, 'Order');    % desired overlap
0155 Lmin = find(pl, 'Lmin'); % Minimum segment length
0156 if isempty(Lmin)
0157   Lmin = 0;
0158 end
0159 
0160 % Debug
0161 debug = find(pl, 'DEBUG');
0162 
0163 %----------------- Get frequency vector
0164 [f, r, m, L, K] = ltpda_ltf_plan(lmin, fsmax, Nolap, 1, Lmin, Jdes, Kdes);
0165 % [f,r,m,L] = ltpda_compute_f(fsmax, lmin, Kdes, Kmin, Jdes, Nolap, debug);
0166 
0167 %----------------- compute CPSD Estimates
0168 
0169 Txy = mlcohere(iS, f, r, m, L, fsmax, Win, Order, Nolap*100, Lmin);
0170 
0171 
0172 %   if nargout == 2
0173 %     disp('** Computing Coherence matrix with variance estimates')
0174 %     tic
0175 %     [Txy, ENBW, Sig] = mlcohere(iS, f, r, m, L, fsmax, Win, Order);
0176 %     toc
0177 %   else
0178 %     disp('** Computing Coherence matrix')
0179 %     tic
0180 %     [Txy, ENBW]  = mlcohere(iS, f, r, m, L, fsmax, Win, Order);
0181 %     toc
0182 %   end
0183 
0184 %----------------- Build output Matrix of AOs
0185 disp('** Building output AO matrix')
0186 bs  = cell(na);
0187 for j=1:na  % Loop over input signal
0188   bi = as(j);
0189   for k=1:na  % Loop over output signal
0190     bo = as(k);
0191     % create new output fsdata
0192     fs = fsdata(f, squeeze(Txy(j,k,:)), fsmax);
0193     fs = set(fs, 'name', sprintf('LCOHERE(%s->%s)', bi.data.name, bo.data.name));
0194     fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0195     % create new output history
0196     h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0197     if na == length(invars)
0198       h = set(h, 'invars', invars([j k]));
0199       nk = invars{k};
0200       nj = invars{j};
0201     else
0202       nk = bo.name;
0203       nj = bi.name;
0204     end
0205     % make output analysis object
0206     b = ao(fs, h);
0207     % set name
0208     b = setnh(b, 'name', sprintf('LCOHERE(%s->%s)', nj, nk));
0209     % add to outputs
0210     bs(j,k) = {b};
0211   end % End loop over output signal
0212 end % End loop over input signal
0213 
0214 % now we have a cell matrix of AOs but
0215 % we want a normal matrix
0216 b  = [bs{:}];
0217 bs = reshape(b, na, na);
0218 varargout{1} = bs;
0219 
0220 %--------------------------------------------------------------------------
0221 % Returns the length of the shortest vector in samples
0222 function lmin = findShortestVector(as)
0223 
0224 lmin = 1e20;
0225 for j=1:length(as)
0226   if len(as(j)) < lmin
0227     lmin = len(as(j));
0228   end
0229 end
0230 
0231 
0232 %--------------------------------------------------------------------------
0233 % Returns the max Fs of a set of AOs
0234 function fs = findFsMax(as)
0235 
0236 fs = 0;
0237 for j=1:length(as)
0238   a = as(j);
0239   if a.data.fs > fs
0240     fs = a.data.fs;
0241   end
0242 end
0243 
0244 %--------------------------------------------------------------------------
0245 % Get default params
0246 function plo = getDefaultPL()
0247 
0248 disp('* creating default plist...');
0249 plo = plist();
0250 plo = append(plo, param('Kdes', 100));
0251 plo = append(plo, param('Jdes', 1000));
0252 plo = append(plo, param('Kmin', 2));
0253 plo = append(plo, param('Lmin', 0));
0254 plo = append(plo, param('Win',  specwin('Kaiser', 1, 200)));
0255 plo = append(plo, param('Nolap', -1));
0256 plo = append(plo, param('Order', 0));
0257 disp('* done.');
0258 
0259 
0260 %--------------------------------------------------------------------------
0261 % Compute all TF combinations
0262 function varargout = mlcohere(varargin)
0263 
0264 % MLCOHERE
0265 %
0266 % Txy = mlcohere(X,f,r,m,L,fs,win,order,olap,Lmin)
0267 %
0268 % M Hewitson 19-05-07
0269 %
0270 %
0271 
0272 % Get inputs
0273 X     = varargin{1};
0274 f     = varargin{2};
0275 r     = varargin{3};
0276 m     = varargin{4};
0277 L     = varargin{5};
0278 fs    = varargin{6};
0279 win   = varargin{7};
0280 order = varargin{8};
0281 olap  = varargin{9};
0282 Lmin  = varargin{10};
0283 
0284 
0285 % --- Prepare some variables
0286 twopi      = 2.0*pi;
0287 si         = size(X);
0288 nc         = si(1);
0289 nf         = length(f);
0290 Txy        = zeros(nc,nc,nf);
0291 disp_each  = round(nf/100)*10;
0292 minReached = 0;
0293 
0294 % ----- Loop over Frequency
0295 for fi=1:nf
0296 
0297   % compute DFT exponent and window
0298   l = L(fi);
0299   
0300   % Check if we need to update the window values
0301   % - once we reach Lmin, the window values don't change.
0302   if ~minReached
0303     switch win.name
0304       case 'Kaiser'
0305         win = specwin(win.name, l, win.psll);
0306       otherwise
0307         win = specwin(win.name, l);
0308     end
0309     if l==Lmin
0310       minReached = 1;
0311     end
0312   end
0313   
0314   % Compute DFT coefficients
0315   p     = 1i * twopi * m(fi)/l.*[0:l-1];
0316   C     = win.win .* exp(p);
0317 
0318   if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0319     disp(sprintf('++ computing frequency %04d of %04d: %f Hz', fi, nf, f(fi)));
0320   end
0321 
0322   % Loop over input channels
0323   for j=1:nc
0324     % loop over output channels
0325     for k=1:nc
0326       if k>j         
0327         % Core cross-DFT part in C-mex file
0328         xi    = X(j,:);
0329         xo    = X(k,:);
0330         % We need cross-spectrum and Power spectrum
0331         [XY, XX, YY, nsegs] = ltpda_dft(xi, xo, l, C, olap, order);
0332         Txy(j,k,fi)  = (abs(XY).^2)/(XX.*YY);
0333       else
0334         % Trivial case of equal input and output
0335         Txy(j,k,fi)  = 1;
0336       end
0337     end % End loop over output channels
0338   end % End loop over input channels
0339 end % End loop over frequencies
0340 
0341 % Fill other elements
0342 for j=1:nc
0343   for k=1:nc
0344     if k<j
0345       Txy(j,k,:) = Txy(k,j,:);
0346     end
0347   end
0348 end
0349 
0350 % Set output
0351 varargout{1} = Txy;
0352 
0353 
0354 % END

Generated on Mon 31-Mar-2008 12:20:24 by m2html © 2003