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

Generated on Mon 02-Jul-2007 12:19:41 by m2html © 2003