Home > wafo > trgauss > trangood.m

trangood

PURPOSE ^

Makes a transformation that is suitable for efficient transforms.

SYNOPSIS ^

f = trangood(ff,nmin,mini,maxi,nmax)

DESCRIPTION ^

 TRANGOOD Makes a transformation that is suitable for efficient transforms.
 
   CALL:  f = trangood(ff,nmin,mini,maxi,nmax);
 
         f    = the good transform function, [X f(X)].
         ff   = the input transform function, [X f(X)]. 
         nmin = the minimum number of points in the good transform.
                (Default  size(ff,1))
         mini = the minimum data value to transform. 
                (Default  min(ff(:,1)))
         maxi = the maximum data value to transform.
                (Default  max(ff(:,1)))
         nmax = then maximum number of points in the good transform
                (default inf)
 
  TRANGOOD interpolates ff linearly  and optionally
   extrapolate ff linearly outside the range of ff(:,1)
   with X uniformly spaced.
 
  See also   tranproc, interp1q

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function f = trangood(ff,nmin,mini,maxi,nmax)
002 %TRANGOOD Makes a transformation that is suitable for efficient transforms.
003 %
004 %  CALL:  f = trangood(ff,nmin,mini,maxi,nmax);
005 %
006 %        f    = the good transform function, [X f(X)].
007 %        ff   = the input transform function, [X f(X)]. 
008 %        nmin = the minimum number of points in the good transform.
009 %               (Default  size(ff,1))
010 %        mini = the minimum data value to transform. 
011 %               (Default  min(ff(:,1)))
012 %        maxi = the maximum data value to transform.
013 %               (Default  max(ff(:,1)))
014 %        nmax = then maximum number of points in the good transform
015 %               (default inf)
016 %
017 % TRANGOOD interpolates ff linearly  and optionally
018 %  extrapolate ff linearly outside the range of ff(:,1)
019 %  with X uniformly spaced.
020 %
021 % See also   tranproc, interp1q
022 
023 
024 % History:
025 % revised pab 07.02.2001
026 % - added nmax
027 % - replaced x = linspace(f(1,1),f(nf,1),nmin)' with x = (f(1,1):df:f(nf,1))'
028 % revised pab 07.01.2001
029 %  -fixed a bug: x = linspace(f(1,1),f(nf,1),nmin)' is safer than using
030 %    x=f(1,1)+(0:nmin-1)'/(nmin-1)*(f(nf,1)-f(1,1));  with interp1q
031 % revised pab 12.11.2000
032 %  - updated header info: A more detailed description of what TRANGOOD does.
033 %  - changed interpolation with a call to interp1q which is much faster
034 %  added nargchk and isempty(maxi),....isempty(nmin)
035 % by ???
036 
037 error(nargchk(1,5,nargin))
038 if (size(ff,2)~=2)
039   error('ff  must be a two column matrix.')
040 end
041 if (size(ff,1)<2)
042   error('ff  must have at least two rows.')
043 end
044 
045 [f,i] = sort(ff(:,1));
046 f     = [f ff(i,2)];
047 clear i;
048 df    = diff(f(:,1));
049 if ( any(df<=0)), %eps 
050   error('Duplicate x-values in  ff  not allowed.')
051 end
052 
053 nf = size(f,1);
054 if (nargin<5)|isempty(nmax),  nmax = inf; end
055 if (nargin<4)|isempty(maxi),  maxi = f(nf,1); end
056 if (nargin<3)|isempty(mini),  mini = f(1,1);  end
057 if (nargin<2)|isempty(nmin),  nmin = nf;      end
058 if (nmin<2),    nmin = 2;  end
059 if (nmax<2),    nmax = 2;  end
060 
061 if ( (nf<nmin) |(nmax<nf) | any(abs(diff(df))>10*eps*(f(nf,1)-f(1,1))) )
062   % pab 07.01.2001: Always choose the stepsize df so that 
063   % it is an exactly representable number.
064   % This is important when calculating numerical derivatives and is 
065   % accomplished by the following.
066   df = (f(nf,1)-f(1,1))/(min(nmin,nmax)-1);
067   df = donothing(df+2)-2;
068   x = (f(1,1):df:f(nf,1)).';
069   % New call pab 11.11.2000: This is much quicker
070   f = [ x interp1q(f(:,1),f(:,2),x)];     
071   %f = [ x interp1(f(:,1),f(:,2),x,'linear')];  
072 end
073 % f(:,1) is now uniformly spaced
074 df = f(2,1)-f(1,1);
075 
076 % Extrapolate linearly outside the range of ff
077 %---------------------------------------------- 
078 if (mini<f(1,1)),
079   f1 = df*(floor((mini-f(1,1))/df):1:-1)';
080   f2 = f(1,2)+f1*(f(2,2)-f(1,2))/(f(2,1)-f(1,1));
081   f  = [f1+f(1,1) f2;f];
082 end
083 n = size(f,1);
084 if (maxi>f(n,1))
085   f1 = df*(1:1:ceil((maxi-f(n,1))/df))';
086   f2 = f(n,2)+f1*(f(n,2)-f(n-1,2))/(f(n,1)-f(n-1,1));
087   f  = [f;f1+f(n,1) f2];
088 end
089 
090 return
091 function y=donothing(x)
092   y=x;
093 return
094 
095 % Old call: (Saved it just in case...)
096  x=f(1,1)+(0:nmin-1)'/(nmin-1)*(f(nf,1)-f(1,1));
097  % "y=interp1(f(:,1),f(:,2),x,'linear');" is slow:
098  % Use the fact that transforms often are "continuous" to 
099  % find y_k=f(x_k) incrementally.
100  y=zeros(nmin,1);
101  i = 1;
102  for k=1:nmin
103    xx=x(k);
104    if (xx>f(i+1,1))
105      while (xx>f(i+1,1))
106        i=i+1;
107        if (i>=nf), i=nf-1; break;  end;
108      end
109    else % xx<=f(i+1,1)
110      while (xx<=f(i,1))
111        i=i-1;
112        if (i<1),   i=1;    break;  end;
113      end
114    end
115    x0 = f(i,1); x1 = f(i+1,1);
116    y0 = f(i,2); y1 = f(i+1,2);
117    y(k) = (xx-x0)*(y1-y0)/(x1-x0)+y0;
118  end
119  f=[x y];
120  clear x y;
121 
122 
123

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003