Home > m > sigproc > frequency_domain > ltpda_lcpsd.m

ltpda_lcpsd

PURPOSE ^

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

SYNOPSIS ^

function varargout = ltpda_lcpsd(varargin)

DESCRIPTION ^

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

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

 Usage:
   >> txy = ltpda_lcpsd(as)
   >> txy = ltpda_lcpsd(as, pl)

   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: fs/4]
     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_lcpsd('Params')

 The following call returns a string that contains the routine CVS version:

 >> version = ltpda_lcpsd(ao,'Version')

 The following call returns a string that contains the routine category:

 >> category = ltpda_lcpsd(ao,'Category')

 M Hewitson 02-02-07

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

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