Home > m > sigproc > frequency_domain > ltpda_ltfe.m

ltpda_ltfe

PURPOSE ^

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

SYNOPSIS ^

function varargout = ltpda_ltfe(varargin)

DESCRIPTION ^

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

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

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