Home > atmlab > covmat > covmat1d_from_cfun.m

covmat1d_from_cfun

PURPOSE ^

COVMAT1D_FROM_CFUN Correlation function based covariance matrix

SYNOPSIS ^

function S = covmat1d_from_cfun(xp,Std,cfun,varargin)

DESCRIPTION ^

 COVMAT1D_FROM_CFUN   Correlation function based covariance matrix

    This function sets up a covariance matrix from a defined correlation 
    function. The correlation function is specified by giving a functional
    form (such as exponential decreasing) and correlation lengths. The
    correlation length is throughout defined as the distance where the 
    correlation has dropped to exp(-1). For off-diagonal values, the 
    correlation length is averaged between the two involved positions.

    Correlation matrices are obtained by setting *Std* to [].

 FORMAT   S = covmat1d_from_cfun( xp, Std, cfun, Cl, [, cco, mapfun] )
        
 OUT   S      The covariance matrix
 IN    xp     The data abscissa.
       Std    Standard deviations. Given as a two vector matrix. First column
              holds position in same unit as *xp*. The second column is the 
              standard deviation at the postions of the first column. These
              values are then interpolated to *xp*, extrapolating end
              values to +-Inf (in a "nearest" manner).
              If set to a scalar, this value is applied for all *xp*.
              If set to [], unit standard deviation is assumed.
       cfun   Correlation function. Possible choices are
               'drc' : Dirac. No correlation. Any given correlation length
                       is ignored here.
               'lin' : Linearly decreasing (down to zero).
               'exp' : Exponential decreasing (exp(-dx/cl)).
               'gau' : Gaussian (normal) deceasing (exp(-(dx/cl))^2).
 OPT   Cl     Correlation lengths. Given as a column matrix as *Std*.
              Must be given for all *cfun* beside 'drc'. Extrapolation as
              for *Std*. Scalar input is allowed.
       cco    Correlation cut-off. All values below this limit are set to 0.
       mapfun Mapping function from grid unit to unit for corrleation
              lengths. For example, if correlation lengths are given in
              pressure decades, while the basic coordinate is Pa, this is
              *mapfun* handled by setting *mapfun* to @log10.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

covmat1d_from_cfun.m

SOURCE CODE ^

0001 % COVMAT1D_FROM_CFUN   Correlation function based covariance matrix
0002 %
0003 %    This function sets up a covariance matrix from a defined correlation
0004 %    function. The correlation function is specified by giving a functional
0005 %    form (such as exponential decreasing) and correlation lengths. The
0006 %    correlation length is throughout defined as the distance where the
0007 %    correlation has dropped to exp(-1). For off-diagonal values, the
0008 %    correlation length is averaged between the two involved positions.
0009 %
0010 %    Correlation matrices are obtained by setting *Std* to [].
0011 %
0012 % FORMAT   S = covmat1d_from_cfun( xp, Std, cfun, Cl, [, cco, mapfun] )
0013 %
0014 % OUT   S      The covariance matrix
0015 % IN    xp     The data abscissa.
0016 %       Std    Standard deviations. Given as a two vector matrix. First column
0017 %              holds position in same unit as *xp*. The second column is the
0018 %              standard deviation at the postions of the first column. These
0019 %              values are then interpolated to *xp*, extrapolating end
0020 %              values to +-Inf (in a "nearest" manner).
0021 %              If set to a scalar, this value is applied for all *xp*.
0022 %              If set to [], unit standard deviation is assumed.
0023 %       cfun   Correlation function. Possible choices are
0024 %               'drc' : Dirac. No correlation. Any given correlation length
0025 %                       is ignored here.
0026 %               'lin' : Linearly decreasing (down to zero).
0027 %               'exp' : Exponential decreasing (exp(-dx/cl)).
0028 %               'gau' : Gaussian (normal) deceasing (exp(-(dx/cl))^2).
0029 % OPT   Cl     Correlation lengths. Given as a column matrix as *Std*.
0030 %              Must be given for all *cfun* beside 'drc'. Extrapolation as
0031 %              for *Std*. Scalar input is allowed.
0032 %       cco    Correlation cut-off. All values below this limit are set to 0.
0033 %       mapfun Mapping function from grid unit to unit for corrleation
0034 %              lengths. For example, if correlation lengths are given in
0035 %              pressure decades, while the basic coordinate is Pa, this is
0036 %              *mapfun* handled by setting *mapfun* to @log10.
0037 
0038 % 2005-05-20   Created by Patrick Eriksson.
0039 
0040 
0041 function S = covmat1d_from_cfun(xp,Std,cfun,varargin)
0042 %
0043 [Cl,cco,mapfun] = optargs( varargin, { [], 0, [] } );
0044 
0045 %= Check input                                                              %&%
0046 %                                                                           %&%
0047 rqre_nargin( 3, nargin );                                                   %&%
0048 rqre_datatype( xp, @istensor1 );                                            %&%
0049 if ~istensor2(Std)                                                          %&%
0050   error( 'Argument *Std* must be a 2-col matrix, a scalar or be empty.' );  %&%
0051 end                                                                         %&%
0052 if ~ischar(cfun)  |  length(cfun) ~= 3                                      %&%
0053   error( 'Argument *cfun* must be a string of length 3.' );                 %&%
0054 end                                                                         %&%
0055 if nargin >= 4  &  ~istensor2(Cl)                                           %&%
0056   error( 'Argument *Cl* must be a matrix with 2 columns.' );                %&%
0057 end                                                                         %&%
0058 rqre_datatype( cco, @istensor0 );                                           %&%
0059 if cco < 0  |  cco >= 1                                                     %&%
0060   error( 'Argument *cco* must be a scalar [0,1[' );                         %&%
0061 end                                                                         %&%
0062 
0063 
0064 %- Determine standard deviations
0065 %
0066 n  = length( xp );
0067 %
0068 if isempty( Std )
0069   si = ones( n, 1 );
0070 elseif isscalar( Std )
0071   si = repmat( Std, n, 1 );
0072 else  
0073   si = interp1( Std(:,1), Std(:,2), handle_expand(Std(:,1),xp) );
0074 end
0075 %                                                                           %&%
0076 if any( isnan(si) )                                                         %&%
0077   error( 'NaN obtained when interpolating *Std* (extrapolation?)' );        %&%
0078 end                                                                         %&%
0079 
0080 
0081 %= Handle diagonal matrices separately (note return)
0082 %
0083 if strcmp( lower(cfun), 'drc' )
0084   %
0085   S = sparse( 1:n, 1:n, si.^2, n, n );
0086   return
0087   %
0088 end
0089 
0090 
0091 %= Conversion of length unit
0092 %
0093 if ~isempty(mapfun)
0094   if ~isa( mapfun, 'function_handle' )
0095     error( 'Input *mapfun* must be empty or a function handle.' );
0096   end
0097   xp = mapfun( xp );
0098   if ~isscalar( Cl )
0099     Cl(:,1) = mapfun( Cl(:,1) );
0100   end
0101 end
0102 
0103 
0104 %= Distance matrix
0105 %
0106 [X1,X2] = meshgrid( xp, xp );
0107 D       = abs( X1 - X2 );
0108 
0109 
0110 %= Correlation length matrix
0111 %
0112 if isscalar( Cl )
0113   L = Cl;
0114 else
0115   cl = interp1( Cl(:,1), Cl(:,2), handle_expand(Cl(:,1),xp) );
0116   [X1,X2] = meshgrid( cl, cl );
0117   L       = ( X1 + X2 ) / 2;
0118   %
0119   clear X1 X2
0120 end
0121 
0122 
0123 %= Create correlation matrix
0124 %
0125 switch lower(cfun)
0126 
0127   case 'lin'
0128     %
0129     S = 1 - (1-exp(-1)) * ( D./L );
0130     % Negativa values removed by cco
0131 
0132   case 'exp'
0133     %
0134     S = exp( -D./L );
0135 
0136   case 'gau'
0137     %
0138     S = exp( -(D./L).^2 );
0139 
0140   otherwise                                                                 %&%
0141     %                                                                       %&%
0142     error( sprintf('Unknown correlation function (%s)',cfun) );             %&%
0143 end
0144 
0145 
0146 %- Remove values below correlation cut-off limit, convert to sparse and
0147 %- include standard deviations
0148 %
0149 S( find( S < cco )) = 0;
0150 %
0151 [i,j,s] = find( S );
0152 S       = (si*si') .* sparse( i, j, s, n, n );
0153 
0154 return
0155 
0156 
0157 
0158 %---
0159 function xi = handle_expand(x,xi)
0160   %
0161   v1 = min( x );
0162   i1 = find( xi<v1 );
0163   %
0164   v2 = max( x );
0165   i2 = find( xi>v2 );
0166   %
0167   if ~isempty(v1) 
0168     xi(i1) = v1;
0169   end
0170   %
0171   if ~isempty(v2) 
0172     xi(i2) = v2;
0173   end
0174   %
0175 return

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005