Home > m > sigproc > frequency_domain > ltpda_ltfe.m

ltpda_ltfe

PURPOSE ^

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

SYNOPSIS ^

function varargout = ltpda_ltfe(varargin)

DESCRIPTION ^

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

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

 Usage:
   >> tfs = ltpda_ltfe(as)
   >> tfs = ltpda_ltfe(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 ltpda_dft
             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')

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

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

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

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

 M Hewitson 02-02-07

 $Id: ltpda_ltfe.html,v 1.14 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_ltfe(varargin)
0002 % LTPDA_LTFE implement transfer-function estimation computed on a log frequency axis.
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % LTPDA_LTFE implement 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 %
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 ltpda_dft
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_ltfe('Params')
0038 %
0039 % The following call returns a string that contains the routine CVS version:
0040 %
0041 % >> version = ltpda_ltfe(ao,'Version')
0042 %
0043 % The following call returns a string that contains the routine category:
0044 %
0045 % >> category = ltpda_ltfe(ao,'Category')
0046 %
0047 % M Hewitson 02-02-07
0048 %
0049 % $Id: ltpda_ltfe.html,v 1.14 2008/03/31 10:27:39 hewitson Exp $
0050 %
0051 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0052 
0053 ALGONAME = mfilename;
0054 VERSION  = '$Id: ltpda_ltfe.html,v 1.14 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_ltfe 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 na = length(as);
0113 clear bs;
0114 disp('*** Done.');
0115 
0116 %----------------- Truncate all vectors
0117 
0118 % Get shortest vector
0119 disp('*** Truncating all vectors...');
0120 lmin = findShortestVector(as);
0121 nsecs = lmin / fsmax;
0122 bs = [];
0123 for i=1:na
0124   a = as(i);
0125   if len(a) ~= lmin
0126     warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0127     bs = [bs  select(a, 1:lmin)];
0128   else
0129     bs = [bs a];
0130   end
0131 end
0132 as = bs;
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 if find(pl, 'OLD SCHEDULER')
0177   [f,r,m,L,rr,rrr] = ltpda_compute_f(fsmax, lmin, Kdes, Kmin, Jdes, Nolap, 0);
0178 else
0179   [f, r, m, L, K] = ltpda_ltf_plan(lmin, fsmax, Nolap, 1, Lmin, Jdes, Kdes);
0180 end
0181 
0182 %----------------- compute TF Estimates
0183 
0184 Txy = mltfe(iS, f, r, m, L, fsmax, Win, Order, Nolap*100, Lmin);
0185 
0186 %----------------- Build output Matrix of AOs
0187 disp('** Building output AO matrix')
0188 bs  = cell(na);
0189 for j=1:na  % Loop over input signal
0190   bi = as(j);
0191   for k=1:na  % Loop over output signal
0192     bo = as(k);
0193     % create new output fsdata
0194     fs = fsdata(f, squeeze(Txy(j,k,:)), fsmax);
0195     fs = set(fs, 'name', sprintf('LTFE(%s->%s)', bi.data.name, bo.data.name));
0196     fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0197     % create new output history
0198     h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0199     try
0200       h = set(h, 'invars', invars([j k]));
0201     catch
0202       warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0203     end
0204     % make output analysis object
0205     b = ao(fs, h);
0206     % set name
0207     try
0208       if isempty(invars{k})
0209         nk = bo.name;
0210       else
0211         nk = invars{k};
0212       end
0213       if isempty(invars{j})
0214         nj = bi.name;
0215       else
0216         nj = invars{j};
0217       end
0218     catch
0219       warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0220       nk = bo.name;
0221       nj = bi.name;
0222     end
0223 
0224     b = setnh(b, 'name', sprintf('LTFE(%s->%s)', nj, nk));
0225     % add to outputs
0226     bs(j,k) = {b};
0227 
0228   end % End loop over output signal
0229 end % End loop over input signal
0230 
0231 % now we have a cell matrix of AOs but
0232 % we want a normal matrix
0233 b  = [bs{:}];
0234 bs = reshape(b, na, na);
0235 varargout{1} = bs;
0236 
0237 %--------------------------------------------------------------------------
0238 % Returns the length of the shortest vector in samples
0239 function lmin = findShortestVector(as)
0240 
0241 lmin = 1e20;
0242 for j=1:length(as)
0243   if len(as(j)) < lmin
0244     lmin = len(as(j));
0245   end
0246 end
0247 
0248 
0249 %--------------------------------------------------------------------------
0250 % Returns the max Fs of a set of AOs
0251 function fs = findFsMax(as)
0252 
0253 fs = 0;
0254 for j=1:length(as)
0255   a = as(j);
0256   if a.data.fs > fs
0257     fs = a.data.fs;
0258   end
0259 end
0260 
0261 %--------------------------------------------------------------------------
0262 % Get default params
0263 function plo = getDefaultPL()
0264 
0265 plo = plist();
0266 plo = append(plo, param('Kdes', 100));
0267 plo = append(plo, param('Jdes', 1000));
0268 plo = append(plo, param('Kmin', 2));
0269 plo = append(plo, param('Lmin', 0));
0270 plo = append(plo, param('Win',  specwin('Kaiser', 1, 150)));
0271 plo = append(plo, param('Nolap', -1));
0272 plo = append(plo, param('Order', 0));
0273 
0274 
0275 %--------------------------------------------------------------------------
0276 % Compute all TF combinations
0277 function varargout = mltfe(varargin)
0278 
0279 % MLTFE
0280 %
0281 % Txy = mltfe(X,f,r,m,L,fs,win)
0282 %
0283 % M Hewitson 19-05-07
0284 %
0285 %
0286 
0287 % Get inputs
0288 X     = varargin{1};
0289 f     = varargin{2};
0290 r     = varargin{3};
0291 m     = varargin{4};
0292 L     = varargin{5};
0293 fs    = varargin{6};
0294 win   = varargin{7};
0295 order = varargin{8};
0296 olap  = varargin{9};
0297 Lmin  = varargin{10};
0298 
0299 
0300 % --- Prepare some variables
0301 twopi      = 2.0*pi;
0302 si         = size(X);
0303 nc         = si(1);
0304 nf         = length(f);
0305 Txy        = zeros(nc,nc,nf);
0306 disp_each  = round(nf/100)*10;
0307 minReached = 0;
0308 
0309 % ----- Loop over Frequency
0310 for fi=1:nf
0311 
0312   % compute DFT exponent and window
0313   l = L(fi);
0314   
0315   % Check if we need to update the window values
0316   % - once we reach Lmin, the window values don't change.
0317   if ~minReached
0318     switch win.name
0319       case 'Kaiser'
0320         win = specwin(win.name, l, win.psll);
0321       otherwise
0322         win = specwin(win.name, l);
0323     end
0324     if l==Lmin
0325       minReached = 1;
0326     end
0327   end
0328   
0329   % Compute DFT coefficients
0330   p     = 1i * twopi * m(fi)/l.*[0:l-1];
0331   C     = win.win .* exp(p);
0332 
0333   if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0334     disp(sprintf('++ computing frequency %d of %d: %f Hz', fi, nf, f(fi)));
0335   end
0336 
0337   % Loop over input channels
0338   for j=1:nc
0339     % loop over output channels
0340     for k=1:nc
0341       if k~=j         
0342         % Core cross-DFT part in C-mex file
0343         xi    = X(j,:);
0344         xo    = X(k,:);
0345         % We need cross-spectrum and Power spectrum
0346         [XY, XX, YY, nsegs] = ltpda_dft(xi, xo, l, C, olap, order);
0347         Txy(j,k,fi)  = XY/XX;
0348       else
0349         % Trivial case of equal input and output
0350         Txy(j,k,fi)  = 1;
0351       end
0352     end % End loop over output channels
0353   end % End loop over input channels
0354 end % End loop over frequencies
0355 
0356 % Set output
0357 varargout{1} = Txy;
0358 
0359 % END

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