Home > m > sigproc > frequency_domain > ltpda_ltfe.m

ltpda_ltfe

PURPOSE ^

LTPDA_LTFE implement weighted transfer-function estimation computed on a log frequency axis.

SYNOPSIS ^

function varargout = ltpda_ltfe(varargin)

DESCRIPTION ^

 LTPDA_LTFE implement weighted transfer-function estimation computed on a log frequency axis.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 LTPDA_LTFE implement weighted transfer-function estimation computed on a 
 log frequency axis.
 
 Usage:
   >> tfs = ltpda_ltfe(as)
   >> tfs = ltpda_ltfe(as, pl)
   >> [tfs, sig2] = ltpda_ltfe(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]
     Kmin  - minimum number of averages   [default: 10]
     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_ltfe('Params')
 
 
 M Hewitson 02-02-07
 
 $Id: ltpda_ltfe.m,v 1.9 2007/09/01 10:18:06 hewitson Exp $
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_ltfe(varargin)
0002 % LTPDA_LTFE implement weighted transfer-function estimation computed on a log frequency axis.
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % LTPDA_LTFE implement weighted transfer-function estimation computed on a
0006 % log frequency axis.
0007 %
0008 % Usage:
0009 %   >> tfs = ltpda_ltfe(as)
0010 %   >> tfs = ltpda_ltfe(as, pl)
0011 %   >> [tfs, sig2] = ltpda_ltfe(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 %     Kmin  - minimum number of averages   [default: 10]
0023 %     Jdes  - number of spectral frequencies to compute [default: fs/4]
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_ltfe('Params')
0039 %
0040 %
0041 % M Hewitson 02-02-07
0042 %
0043 % $Id: ltpda_ltfe.m,v 1.9 2007/09/01 10:18:06 hewitson Exp $
0044 %
0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0046 
0047 % Check if this is a call for parameters
0048 if nargin == 1
0049   in = char(varargin{1});
0050   if strcmp(in, 'Params')
0051     varargout{1} = getDefaultPL();
0052     return
0053   end
0054 end
0055 
0056 % capture input variable names
0057 invars = {};
0058 as     = [];
0059 ps     = [];
0060 for j=1:nargin
0061   invars = [invars cellstr(inputname(j))];
0062   if isa(varargin{j}, 'ao')
0063     as = [as varargin{j}];
0064   end
0065   if isa(varargin{j}, 'plist')
0066     ps = [ps varargin{j}];
0067   end
0068 end
0069 
0070 ALGONAME = mfilename;
0071 VERSION  = '$Id: ltpda_ltfe.m,v 1.9 2007/09/01 10:18:06 hewitson Exp $';
0072 
0073 % initialise output array
0074 bs = []; 
0075 
0076 % check plist
0077 if isempty(ps)
0078   pl = getDefaultPL();
0079 else
0080   pl = combine(ps, getDefaultPL);
0081 end
0082 
0083 % Go through each input AO and build a matrix of signals
0084 na    = length(as);
0085 
0086 %----------------- Resample all AOs
0087 
0088 disp('*** Resampling AOs...');
0089 fsmax = findFsMax(as);
0090 fspl  = plist(param('fsout', fsmax));
0091 bs = [];
0092 for i=1:na
0093   a = as(i);
0094   % check this is a time-series object
0095   if ~isa(a.data, 'tsdata')
0096     error('### ltpda_ltfe requires tsdata (time-series) inputs.');
0097   end  
0098   % Check Fs
0099   if a.data.fs ~= fsmax
0100     warning(sprintf('!!! Resampling AO %s/%s to %f Hz', a.name, a.data.name, fsmax));
0101     a = resample(a, fspl);
0102   end
0103   bs = [bs a];
0104 end
0105 as = bs;
0106 clear bs;
0107 disp('*** Done.');
0108 
0109 %----------------- Truncate all vectors
0110 
0111 % Get shortest vector
0112 disp('*** Truncating all vectors...');
0113 lmin = findShortestVector(as);
0114 nsecs = lmin / fsmax;
0115 bs = [];
0116 for i=1:na
0117   a = as(i);
0118   if len(a) ~= lmin
0119     warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0120     bs = [bs  select(a, 1:lmin)];
0121   else
0122     bs = [bs a];
0123   end
0124 end
0125 as = bs;
0126 clear bs;
0127 disp('*** Done.');
0128 
0129 
0130 %----------------- Build signal Matrix
0131 N     = len(as(1)); % length of first signal
0132 iS    = zeros(na, N);
0133 for j=1:na
0134   a = as(j);
0135   iS(j,:) = a.data.x;
0136 end
0137 
0138 %----------------- check input parameters
0139 
0140 % Averaging
0141 Kdes = find(pl, 'Kdes');   % desired averages
0142 Kmin = find(pl, 'Kmin');    % minimum averages
0143 Jdes = find(pl, 'Jdes');    % num desired spectral frequencies
0144 if Jdes == -1
0145   Jdes = round(fsmax/4);
0146   disp(sprintf('! Using default Jdes of %d', Jdes))
0147   pl = set(pl, 'Jdes', Jdes);
0148 end
0149 Win   = find(pl, 'Win');   % Window object
0150 Nolap = find(pl, 'Nolap');    % desired overlap
0151 if Nolap == -1
0152   % use window function rov
0153   Nolap = Win.rov/100;
0154   disp(sprintf('! Using recommended overlap of %d', Nolap))
0155   pl = set(pl, 'Nolap', Nolap);
0156 end
0157 Order = find(pl, 'Order');    % desired overlap
0158 if isempty(Order)
0159   Order = 0;
0160   disp(sprintf('! Using default detrending order of %d', Order))
0161   pl = append(pl, param('Order', Order));
0162 end
0163 
0164 
0165 %----------------- Get frequency vector
0166 [f,r,m,L] = ltpda_compute_f(fsmax, lmin, Kdes, Kmin, Jdes, Nolap);
0167   
0168 %----------------- compute TF Estimates
0169 
0170 if nargout == 2
0171   disp('** Computing TF matrix with variance estimates')
0172   tic
0173   [Txy, ENBW, Sig] = mltfe(iS, f, r, m, L, fsmax, Win, Order);
0174   toc
0175 else
0176   disp('** Computing TF matrix')
0177   tic
0178   [Txy, ENBW] = mltfe(iS, f, r, m, L, fsmax, Win, Order);
0179   toc
0180 end
0181   
0182 %----------------- Build output Matrix of AOs
0183 disp('** Building output AO matrix')
0184 bs  = cell(na);
0185 sig = cell(na);
0186 for j=1:na  % Loop over input signal
0187   bi = as(j);
0188   for k=1:na  % Loop over output signal
0189     bo = as(k);   
0190     % create new output fsdata
0191     fs = fsdata(f, squeeze(Txy(j,k,:)), fsmax);
0192     fs = set(fs, 'name', sprintf('LTFe(%s->%s)', bi.data.name, bo.data.name));
0193     fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0194     fs = set(fs, 'enbw', squeeze(ENBW(j,k,:)));
0195     % create new output history
0196     h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0197     try
0198       h = set(h, 'invars', invars([j k]));
0199     catch
0200       warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0201     end
0202     % make output analysis object
0203     b = ao(fs, h);
0204     % set name
0205     try
0206       if isempty(invars{k})
0207         nk = bo.name;
0208       else
0209         nk = invars{k};
0210       end
0211       if isempty(invars{j})
0212         nj = bi.name;
0213       else
0214         nj = invars{j};
0215       end
0216     catch
0217       warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0218       nk = bo.name;
0219       nj = bi.name;
0220     end
0221 
0222     b = set(b, 'name', sprintf('LTFe(%s->%s)', nj, nk));
0223     % add to outputs
0224     bs(j,k) = {b};
0225     
0226     % Make variance outputs
0227     if nargout == 2
0228       % create new output fsdata
0229       fs = fsdata(f, squeeze(Sig(j,k,:)), fsmax);
0230       fs = set(fs, 'name', ['\sigma^2' sprintf('(%s->%s)', bi.data.name, bo.data.name)]);
0231       fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0232       fs = set(fs, 'enbw', squeeze(ENBW(j,k,:)));
0233       % create new output history
0234       h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0235       try
0236         h = set(h, 'invars', invars([j k]));
0237       catch
0238         warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0239       end
0240       % make output analysis object
0241       b = ao(fs, h);
0242       % set name
0243       try
0244         if isempty(invars{k})
0245           nk = bo.name;
0246         else
0247           nk = invars{k};
0248         end
0249         if isempty(invars{j})
0250           nj = bi.name;
0251         else
0252           nj = invars{j};
0253         end
0254       catch
0255         nk = bo.name;
0256         nj = bi.name;
0257       end
0258       b = set(b, 'name', ['\sigma^2' sprintf('(%s->%s)', nj, nk)]);
0259       % add to outputs
0260       sig(j,k) = {b};
0261     
0262     end
0263   end % End loop over output signal
0264 end % End loop over input signal
0265 
0266 % now we have a cell matrix of AOs but
0267 % we want a normal matrix
0268 b  = [bs{:}];
0269 bs = reshape(b, na, na);
0270 varargout{1} = bs;
0271 if nargout == 2
0272   b  = [sig{:}];
0273   s  = reshape(b, na, na);
0274   varargout{2} = s;
0275 end
0276 %--------------------------------------------------------------------------
0277 % Returns the length of the shortest vector in samples
0278 function lmin = findShortestVector(as)
0279 
0280 lmin = 1e20;
0281 for j=1:length(as)
0282   if len(as(j)) < lmin
0283     lmin = len(as(j));
0284   end
0285 end
0286 
0287 
0288 %--------------------------------------------------------------------------
0289 % Returns the max Fs of a set of AOs
0290 function fs = findFsMax(as)
0291 
0292 fs = 0;
0293 for j=1:length(as)
0294   a = as(j);
0295   if a.data.fs > fs
0296     fs = a.data.fs;
0297   end
0298 end
0299 
0300 %--------------------------------------------------------------------------
0301 % Get default params
0302 function plo = getDefaultPL()
0303 
0304 disp('* creating default plist...');
0305 plo = plist();
0306 plo = append(plo, param('Kdes', 100));
0307 plo = append(plo, param('Jdes', 1000));
0308 plo = append(plo, param('Kmin', 2));
0309 plo = append(plo, param('Win',  specwin('Kaiser', 1, 150)));
0310 plo = append(plo, param('Nolap', -1));
0311 plo = append(plo, param('Order', 0));
0312 disp('* done.');
0313 
0314 
0315 %--------------------------------------------------------------------------
0316 % Compute all TF combinations
0317 function varargout = mltfe(varargin)
0318 
0319 % MLTFE
0320 %
0321 % Txy = mltfe(X,f,r,m,L,fs,win)
0322 %
0323 % M Hewitson 19-05-07
0324 %
0325 %
0326 
0327 % Get inputs
0328 X     = varargin{1};
0329 f     = varargin{2};
0330 r     = varargin{3};
0331 m     = varargin{4};
0332 L     = varargin{5};
0333 fs    = varargin{6};
0334 win   = varargin{7};
0335 order = varargin{8};
0336 
0337 % do we need to compute variance?
0338 if nargout == 2
0339   computeVariance = 1;
0340 else
0341   computeVariance = 0;
0342 end
0343 
0344 twopi    = 2.0*pi;
0345 twopiofs = 2.0*pi/(1.0*fs);
0346 
0347 si   = size(X);
0348 nx   = si(2);
0349 nc   = si(1);
0350 nf   = length(f);
0351 Txy  = zeros(nc,nc,nf);
0352 Sig  = zeros(nc,nc,nf);
0353 ENBW = zeros(nc,nc,nf);
0354 
0355 disp_each = round(nf/100)*10;
0356 
0357 for fi=1:nf
0358   
0359   % compute DFT exponent and window
0360   l = L(fi);
0361   switch win.name
0362     case 'Kaiser'       
0363       win = specwin(win.name, l, win.psll);
0364     otherwise
0365       win = specwin(win.name, l);
0366   end    
0367       
0368   p     = 1i * twopi * m(fi)/l.*[0:l-1];
0369   C     = win.win .* exp(p);
0370   Xolap = (1-win.rov/100);
0371   
0372   if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0373     disp(sprintf('++ computing frequency %d of %d: %f Hz', fi, nf, f(fi)));
0374   end
0375   
0376   % Loop over input channels
0377   for j=1:nc
0378     % loop over output channels
0379     for k=1:nc      
0380       
0381       % Count N segments
0382       start = 1;
0383       nsegs = 0;
0384       while start + l - 1 <= nx
0385         % make next step
0386         start = floor(start + l*Xolap);
0387         nsegs = nsegs + 1;        
0388       end
0389       
0390       % do segments
0391       A   = 0.0;
0392       T   = zeros(nsegs,1);
0393       i2  = zeros(nsegs,1);
0394       
0395       start = 1;
0396       nsegs = 0;
0397       p = 1;
0398       while start + l - 1 <= nx
0399 
0400         % get segments
0401         idx = start:start+l-1;
0402         xi    = X(j,idx).';
0403         xo    = X(k,idx).';
0404         
0405         % detrend segments
0406         switch order
0407           case -1
0408             % do nothing
0409           case 0
0410             xi = xi - mean(xi);
0411             xo = xo - mean(xo);
0412           case 1
0413             xi = detrend(xi);
0414             xo = detrend(xo);
0415           otherwise
0416             xi = polydetrend(xi, order);
0417             xo = polydetrend(xo, order);
0418         end
0419 
0420         % compute DFT
0421         ti    = C*xi;
0422         to    = C*xo;
0423         i2(p) = (ti*conj(ti));
0424         T(p)  = to/ti;
0425         A     = A + T(p)*i2(p);
0426                 
0427         % make next step
0428         start = floor(start + l*Xolap);
0429         nsegs = nsegs + 1;        
0430         p     = p + 1;
0431       end
0432 
0433       sw = sum(i2);
0434       Txy(j,k,fi)  = conj(A/sw);
0435       
0436       S1      = win.ws;
0437       S12     = S1*S1;
0438       S2      = win.ws2;
0439       ENBW(j,k,fi) = fs*S2/S12;      
0440       
0441       % compute variance
0442       if computeVariance
0443         A     = 0;
0444         start = 1;
0445         p     = 1;
0446         while start + l - 1 <= nx
0447 %         for p=1:ns
0448           S  = (T(p)-Txy(j,k,fi))^2;
0449           A  = A + S * i2(p);
0450           % make next step
0451           start = floor(start + l*Xolap);
0452           p     = p + 1;
0453         end
0454         Sig(j,k,fi) = abs(A/sw);        
0455       end
0456       
0457     end % End loop over output channels
0458   end % End loop over input channels
0459 end % End loop over frequencies
0460 
0461 varargout{1} = Txy;
0462 if nargout >= 2
0463   varargout{2} = ENBW;
0464 end
0465 if nargout == 3
0466   varargout{3} = Sig;
0467 end
0468 
0469 % END

Generated on Mon 03-Sep-2007 12:12:34 by m2html © 2003