Mercurial > hg > ltpda

% test csd for ao/noisegen2D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: Output a model for mdc1 cross spectral density
%
% CALL:
%       [TF,CSD] = get_test_csd_ao_noisegen2D(pl)
%
% PARAMETERS: 'fs'      -  Sampling frequency
%             'f1'      -  the start frequency
%             'f2'      -  the stop frequency
%             'nf'      -  number of evaluation points
%             'scale'   - spacing of frequencies: 'lin' or 'log'
%         or
%             'f'       -  a vector of frequency values or an AO
%                          whereby the x-axis is taken for the frequency values
%
% OUTPUTS:
%              CSD     -  a 2x2 matrix of analysis objects containing CSD
%                       models
%              TF      -  a 2x2 matrix of analysis objects containing TF
%                       models
%
% M-FILE INFO: Get information about this methods by calling
%              >> ao.getInfo('get_test_csd_ao_noisegen2D')
%
%              Get information about a specified set-plist by calling:
%              >> ao.getInfo('get_test_csd_ao_noisegen2D', 'none')
%
% VERSION:     $Id: get_test_csd_ao_noisegen2D.m,v 1.4 2011/03/24 19:37:53 ingo Exp $
%
%
% HISTORY:     22-12-2008 L Ferraioli
%                 Creation
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = get_test_csd_ao_noisegen2D(varargin)

  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end

  import utils.const.*
  utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);

  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end

  % Collect all AOs and plists
  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);

  % Get frequency vector
  f = find(pl, 'f');
  if isa(f, 'ao') && (isa(f.data, 'fsdata') || isa(f.data, 'xydata'))
    f = f.data.getX;
  end

  % Get sampling frequency
  fs = find(pl, 'fs');

  if isempty(f)
    % Compute from frequency range
    f1 = find(pl, 'f1');
    f2 = find(pl, 'f2');
    ndata = find(pl, 'nf');
    scale = find(pl, 'scale');

    switch scale
      case 'lin'
        f   = linspace(f1, f2, ndata);
      case 'log'
        f = logspace(log10(f1), log10(f2), ndata);
    end
  end


  % Models response calculation
  tf11 = calcTF11(f,fs);
  tf12 = calcTF12(f,fs);
  tf21 = calcTF21(f,fs);
  tf22 = calcTF22(f,fs);

  % CSD calculation
  csd11 = tf11.*conj(tf11)+tf12.*conj(tf12);
  csd12 = tf11.*conj(tf21)+tf12.*conj(tf22);
  csd22 = tf22.*conj(tf22)+tf21.*conj(tf21);
  csd21 = conj(csd12);
  csd11.setName('CSD11');
  csd12.setName('CSD12');
  csd21.setName('CSD21');
  csd22.setName('CSD22');

  % Output data
  if nargout == 1
    varargout{1} = [csd11 csd12;csd21 csd22];
  elseif nargout == 2
    varargout{1} = [csd11 csd12;csd21 csd22];
    varargout{2} = [tf11 tf12;tf21 tf22];
  end

end

%--------------------------------------------------------------------------
% Get Info Object
%--------------------------------------------------------------------------
function ii = getInfo(varargin)
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pls   = [];
  elseif nargin == 1&& ~isempty(varargin{1}) && ischar(varargin{1})
    sets{1} = varargin{1};
    pls = getDefaultPlist(sets{1});
  else
    sets = {'List', 'Range'};
    pls = [];
    for kk=1:numel(sets)
      pls = [pls getDefaultPlist(sets{kk})];
    end
  end
  % Build info object
  ii = minfo(mfilename, 'ao', '', utils.const.categories.mdc01, '$Id: get_test_csd_ao_noisegen2D.m,v 1.4 2011/03/24 19:37:53 ingo Exp $', sets, pls);
  ii.setModifier(false);
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------
function plo = getDefaultPlist(set)
  switch set
    case 'List'
      plo = plist('fs', 10, 'f', []);
    case 'Range'
      plo = plist('fs', 10, 'f1', -1, 'f2', -1, 'nf', 1000, 'scale', 'log');
    otherwise
      error('### Unknown set [%s] for a default list.', set);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOCAL FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tf11 = calcTF11(f,fs)


  % Model Stefano TF11 coefficients
  dRes11 = [2.44554138162509e-011 - 1.79482547894083e-011i;
  2.44554138162509e-011 + 1.79482547894083e-011i;
  2.66402334803101e-009 +  1.1025122049153e-009i;
  2.66402334803101e-009 -  1.1025122049153e-009i;
  -7.3560293387644e-009;
  -1.82811618589835e-009 - 1.21803627800855e-009i;
  -1.82811618589835e-009 + 1.21803627800855e-009i;
  1.16258677367555e-009;
  1.65216557639319e-016;
  -1.78092396888606e-016;
  -2.80420398962379e-017;
  9.21305973049041e-013 - 8.24686706827269e-014i;
  9.21305973049041e-013 + 8.24686706827269e-014i;
  5.10730060739905e-010 - 3.76571756625722e-011i;
  5.10730060739905e-010 + 3.76571756625722e-011i;
  3.45893698149735e-009;
  3.98139182134446e-014 - 8.25503935419059e-014i;
  3.98139182134446e-014 + 8.25503935419059e-014i;
  -1.40595719147164e-011];

  dPoles11 = [0.843464045655194 - 0.0959986292915475i;
  0.843464045655194 + 0.0959986292915475i;
  0.953187595424927 - 0.0190043625473383i;
  0.953187595424927 + 0.0190043625473383i;
  0.967176277937188;
  0.995012027005247 - 0.00268322602801729i;
  0.995012027005247 + 0.00268322602801729i;
  0.996564761885673;
  0.999999366165445;
  0.999981722418555;
  0.999921882627659;
  0.999624431675213 - 0.000813407848742761i;
  0.999624431675213 + 0.000813407848742761i;
  0.997312006278751 - 0.00265611346834941i;
  0.997312006278751 + 0.00265611346834941i;
  0.990516544257531;
  0.477796923118318 - 0.311064085401834i;
  0.477796923118318 + 0.311064085401834i;
  0];

  dDTerms11 = 0;

  % response calculation
  pfparams.type = 'disc';
  pfparams.freq = f;
  pfparams.fs = fs;
  pfparams.res = dRes11;
  pfparams.pol = dPoles11;
  pfparams.dterm = dDTerms11;
  pfr = utils.math.pfresp(pfparams);
  mtf11 = pfr.resp;

  % building AOs
  tf11 = ao(plist('xvals', f, 'yvals', mtf11, 'fs', fs, 'type', 'fsdata'));
  % name for this object
  tf11.setName('TF11');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tf12 = calcTF12(f,fs)

  % Model Stefano TF12
  dRes12 = [1.44258422208796e-017 + 7.07359428613009e-019i;
  1.44258422208796e-017 - 7.07359428613009e-019i;
  -3.4918408053655e-021 - 1.05662874569329e-021i;
  -3.4918408053655e-021 + 1.05662874569329e-021i;
  -7.61773292876976e-021;
  4.84357724603939e-020 + 2.38824204294595e-019i;
  4.84357724603939e-020 - 2.38824204294595e-019i;
  -4.07088520945753e-020 - 2.31474543846105e-019i;
  -4.07088520945753e-020 + 2.31474543846105e-019i;
  8.73316588658882e-023;
  -5.21840635377469e-020;
  1.8461911504859e-023;
  5.20105247464461e-020;
  -4.68960092394415e-022;
  -1.44261407664171e-017 +  6.8922564526833e-019i;
  -1.44261407664171e-017 -  6.8922564526833e-019i;
  3.13688133935426e-022];

  dPoles12 = [0.477546340377332 -     0.310830571032376i;
            0.477546340377332 +     0.310830571032376i;
             0.99790715414307 -    0.0028490561287024i;
             0.99790715414307 +    0.0028490561287024i;
            0.998014205354671                         ;
            0.999585354543332 -  0.000780408757425194i;
            0.999585354543332 +  0.000780408757425194i;
             0.99966003029931 -  0.000830944038363768i;
             0.99966003029931 +  0.000830944038363768i;
            0.999962770401331                         ;
            0.999981881865521                         ;
            0.999999365763457                         ;
            0.999981706320212                         ;
             0.99992421574188                         ;
            0.477898460791003 +     0.311001926610074i;
            0.477898460791003 -     0.311001926610074i;
                            0];
  dDTerms12 = 0;

  % response calculation
  pfparams.type = 'disc';
  pfparams.freq = f;
  pfparams.fs = fs;
  pfparams.res = dRes12;
  pfparams.pol = dPoles12;
  pfparams.dterm = dDTerms12;
  pfr = utils.math.pfresp(pfparams);
  mtf12 = pfr.resp;

  % building AOs
  tf12 = ao(plist('xvals', f, 'yvals', mtf12, 'fs', fs, 'type', 'fsdata'));
  % name for this object
  tf12.setName('TF12');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tf21 = calcTF21(f,fs)

  % Model Stefano Tf21
  dRes21 = [-1.80035241582968e-016 + 1.99543917791863e-015i;
       -1.80035241582968e-016 - 1.99543917791863e-015i;
       -1.85590889333759e-013 - 1.23844418827409e-014i;
       -1.85590889333759e-013 + 1.23844418827409e-014i;
        5.03656596876842e-013                         ;
       -2.62470963499904e-013 + 2.30024232938878e-012i;
       -2.62470963499904e-013 - 2.30024232938878e-012i;
       -9.83780507870955e-018                         ;
        3.40426735130194e-021                         ;
        9.78322351492755e-018                         ;
       -1.65010934542937e-020                         ;
        2.60918565203438e-015 +  1.0546609464659e-015i;
        2.60918565203438e-015 -  1.0546609464659e-015i;
        3.18105585405455e-014 + 2.48839990780042e-013i;
        3.18105585405455e-014 - 2.48839990780042e-013i;
        3.23021641947666e-013                         ;
        4.81265000078114e-016 - 3.18269170053848e-017i;
        4.81265000078114e-016 + 3.18269170053848e-017i;
        5.16260024128201e-018];

  dPoles21 = [0.872004077421604 -     0.110344282822693i;
            0.872004077421604 +     0.110344282822693i;
            0.956884129232757 -    0.0225532091775074i;
            0.956884129232757 +    0.0225532091775074i;
            0.966514825697177                         ;
            0.995140550419744 -  0.000755127639524413i;
            0.995140550419744 +  0.000755127639524413i;
            0.999981802194393                         ;
             0.99999936576546                         ;
            0.999981722418555                         ;
            0.999921882627659                         ;
            0.999624431675213 +  0.000813407848742761i;
            0.999624431675213 -  0.000813407848742761i;
            0.997312006278751 +   0.00265611346834941i;
            0.997312006278751 -   0.00265611346834941i;
            0.990516544257531                         ;
            0.477796923118318 +     0.311064085401834i;
            0.477796923118318 -     0.311064085401834i;
                            0];

  dDTerms21 = 0;

  % response calculation
  pfparams.type = 'disc';
  pfparams.freq = f;
  pfparams.fs = fs;
  pfparams.res = dRes21;
  pfparams.pol = dPoles21;
  pfparams.dterm = dDTerms21;
  pfr = utils.math.pfresp(pfparams);
  mtf21 = pfr.resp;

  % building AOs
  tf21 = ao(plist('xvals', f, 'yvals', mtf21, 'fs', fs, 'type', 'fsdata'));
  % name for this object
  tf21.setName('TF21');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tf22 = calcTF22(f,fs)

  % Model Stefano Tf22
  dRes22 = [1.1284521501259e-014;
       -3.72133611555879e-014 - 2.08232683444075e-014i;
       -3.72133611555879e-014 + 2.08232683444075e-014i;
        9.84930639106637e-014 - 1.46640810672565e-013i;
        9.84930639106637e-014 + 1.46640810672565e-013i;
        2.51323684013671e-014                         ;
       -5.64078525288305e-014                         ;
       -1.62476406586366e-014                         ;
       -1.08424815979566e-011 + 8.32328079357669e-012i;
       -1.08424815979566e-011 - 8.32328079357669e-012i;
        2.41831559776112e-011];

  dPoles22 = [0.988511243978897;
            0.997305870640646 +   0.00211760900132725i;
            0.997305870640646 -   0.00211760900132725i;
            0.999626453270255 +    0.0008125673525946i;
            0.999626453270255 -    0.0008125673525946i;
            0.999999366366222                         ;
            0.999981706320212                         ;
             0.99992421574188                         ;
            0.477898460791003 -     0.311001926610074i;
            0.477898460791003 +     0.311001926610074i;
                            0];

  dDTerms22 = 0;

  % response calculation
  pfparams.type = 'disc';
  pfparams.freq = f;
  pfparams.fs = fs;
  pfparams.res = dRes22;
  pfparams.pol = dPoles22;
  pfparams.dterm = dDTerms22;
  pfr = utils.math.pfresp(pfparams);
  mtf22 = pfr.resp;

  % building AOs
  tf22 = ao(plist('xvals', f, 'yvals', mtf22, 'fs', fs, 'type', 'fsdata'));
  % name for this object
  tf22.setName('TF22');

end
author	Daniele Nicolodi <nicolodi@science.unitn.it>
date	Tue, 06 Dec 2011 18:42:11 +0100
parents
children