Initial import of Taylor Matlab Toolbox. This is v1.10.

Author: taylora

689_032_0002.abf

@@ -0,0 +1,4 @@
+function f_alias=f(f_true,f_sampling)
+
+f_alias=abs(f_true-f_sampling*round(f_true/f_sampling));
+

basics/alias_frequency.m

@@ -0,0 +1,11 @@
+function x = f(y)
+
+% all elements of y must be in the interval (0,1)
+
+dims=size(y);
+lower_bound=repmat(1e-15,dims);
+upper_bound=repmat(1-1e-15,dims);
+y=max(y,lower_bound);
+y=min(y,upper_bound);
+x = -log(y.^(-1)-1);
+

basics/antilogistic.m

@@ -0,0 +1,10 @@
+function rev_n=f(n)
+
+% i should be a uint8
+% this is probably hecka slow
+% if you're using this, you should think about using bitrevorder
+
+rev_n=uint8(0);
+for j=1:8
+  rev_n=bitset(rev_n,j,bitget(n,9-j));
+end

basics/bit_reverse.m

@@ -0,0 +1,149 @@
+function [x_seg,y_normed_seg,y_eb_normed_seg] = ...
+  break_at_wrap_points(x,y,y_lim,y_eb)
+
+% this function breaks an angular quantity y wherever it passes the "wrap
+% points" in y_lim.  For instance, y might be an angle in radians, and 
+% y_lim might equal [-pi +pi].  It returns x_seg and y_normed_seg, both
+% cell arrays of "segments".  If y crosses the wrap limits n_break times, 
+% then there will be n_break+1 segments, each containing a stretch of
+% values where y did not cross the wrap limits.  Any y_normed_seg{i}(j), 
+% for all i and j, will be s.t. y_lim(1) <= y_normed_seg{i}(j) <= y_lim(2).
+% This is useful is you would like to plot an angular quantity, for
+% instance.  The function is smart about interpolating near the break
+% points so that the segments end exactly at the y limits, not just near
+% them.  If the y_eb argument is provided, this error bar is also wrapped
+% appropriately, and y_eb_normed_seg{i}(j,:) is the wrapped error bar for
+% y_eb_normed_seg{i}(j).  Note however that the error bar limits can lie
+% outside of [y_lim(1),y_lim(2)].  (This is usually a good thing for
+% plotting purposes.)
+%
+% we assume x, y, are col vectors of length n
+% we assume y_lim is 2x1, with y_lim(1)<y_lim(2)
+% if present, we assume y_eb is n x 2, with 
+%  y_eb(i,1) <= y_eb(i,2) for all i
+
+% deal with args
+calc_error_bars=(nargin>=4);
+if nargin<4 
+  calc_error_bars=false;
+  y_eb_normed_seg=[];
+else
+  calc_error_bars=true;
+end
+
+% figure out where the breaks are
+y_lo=y_lim(1);
+y_hi=y_lim(2);
+y_span=y_hi-y_lo;
+y_shift=y-y_lo;
+y_phase=y_shift/y_span;  
+  % y_phase takes on integer values at wrap points
+y_n_wraps=floor(y_phase);
+  % y_n_wraps is the (signed) number of times you pass the wrap point in 
+  % going from between y_lo and y_hi to that value of y.
+i_break=find(diff(y_n_wraps)~=0);
+  % there's a break between i_break(j) and i_break(j)+1, for all j
+n_break=length(i_break);
+n_seg=n_break+1;  % a "seg" (==segment) is a run between breaks
+  % there's also a seg before the first break, and after the last
+
+% calculate the x coordinate of each break
+x_break=zeros(n_break,1);
+for i=1:n_break
+  j=i_break(i);
+  x_break(i)=interp1([y_phase(j);y_phase(j+1)],...
+                     [x(j);x(j+1)],...
+                     max(y_n_wraps(j),y_n_wraps(j+1)));
+end
+
+%
+% useful quantities for making the segments
+%
+n_wraps_seg=[y_n_wraps(i_break);y_n_wraps(end)];
+  % n_wraps_seg gives the value of y_n_wraps for each seg (y_n_wraps is 
+  % constant within a seg)
+
+% calculate the offset of the error bars from y
+if calc_error_bars
+  dy_eb=y_eb-repmat(y,[1 2]);
+end
+
+% calc the version of y restricted to the bounds
+y_normed      =y      -y_n_wraps*y_span;  % y_lo <= this < y_hi
+if calc_error_bars
+  y_eb_normed=repmat(y_normed,[1 2])+dy_eb;
+end
+
+% calculate the value of y on either side of each break
+y_normed_break_pre =y_span*(1+diff(n_wraps_seg))/2+y_lo;
+y_normed_break_post=y_span*(1-diff(n_wraps_seg))/2+y_lo;
+if isempty(y_normed_break_pre)
+  % this matters downstream
+  y_normed_break_pre =zeros(0,1);
+  y_normed_break_post=zeros(0,1);  
+end
+
+% calculate the offset of the error bars from y at each break
+if calc_error_bars
+  y_break=interp1(x,y,x_break);
+  y_eb_break=interp1(x,y_eb,x_break);
+  dy_eb_break=y_eb_break-repmat(y_break,[1 2]);
+end
+
+% calculate the error bars on either side of each break
+if calc_error_bars
+  y_eb_normed_break_pre =repmat(y_normed_break_pre ,[1 2])+dy_eb_break;
+  y_eb_normed_break_post=repmat(y_normed_break_post,[1 2])+dy_eb_break;
+end
+
+% make the segments
+x_seg=cell(n_seg,1);
+y_normed_seg=cell(n_seg,1);
+for i=1:n_seg
+  if i==1 && i==n_seg
+    % i.e., if there are no breaks
+    x_seg{i}=x;
+    y_normed_seg{i}=y_normed;
+    if calc_error_bars
+      y_eb_normed_seg{i}=y_eb_normed;
+    end
+  elseif i==1 
+    x_seg{i}=[x(1:i_break(1)) ; ...
+              x_break(1)];
+    y_normed_seg{i}=[y_normed(1:i_break(1)) ; ...
+                     y_normed_break_pre(1)];
+    if calc_error_bars               
+      dy_eb_normed_seg_this=...
+        [dy_eb(1:i_break(1),:) ; ...
+         dy_eb_break(1,:)    ];
+      y_eb_normed_seg{i}=repmat(y_normed_seg{i},[1 2])+...
+                         dy_eb_normed_seg_this; 
+    end
+  elseif i==n_seg
+    x_seg{i}=[x_break(end) ; ...
+              x(i_break(end)+1:end)];
+    y_normed_seg{i}=[y_normed_break_post(end) ; ...
+                     y_normed(i_break(end)+1:end)];
+    if calc_error_bars               
+      dy_eb_normed_seg_this=...
+        [dy_eb_break(end,:) ; ...
+         dy_eb(i_break(end)+1:end,:) ]; ...
+      y_eb_normed_seg{i}=repmat(y_normed_seg{i},[1 2])+...
+                         dy_eb_normed_seg_this; 
+    end
+  else
+    x_seg{i}=[x_break(i-1);...
+              x(i_break(i-1)+1:i_break(i));...
+              x_break(i)];
+    y_normed_seg{i}=[y_normed_break_post(i-1); ...
+                     y_normed(i_break(i-1)+1:i_break(i)); ...
+                     y_normed_break_pre(i)];
+    if calc_error_bars               
+      dy_eb_normed_seg_this=[dy_eb_break(i-1,:); ...
+                             dy_eb(i_break(i-1)+1:i_break(i),:); ...
+                             dy_eb_break(i,:)];
+      y_eb_normed_seg{i}=repmat(y_normed_seg{i},[1 2])+...
+                         dy_eb_normed_seg_this;
+    end
+  end
+end

basics/break_at_wrap_points.m

@@ -0,0 +1,2 @@
+clear;
+close all;

basics/cclf.m

@@ -0,0 +1,13 @@
+function result = f(data)
+
+% data must be a vector
+% result is square root of the "S squared" stat over the "x bar" stat
+% I don't really know much about what you can say about how this relates to
+% sqrt(Var[X]])/E[X].  Is it biased?  Are there other estimators of this
+% quantity that strictly dominate it in terms of squared-error risk?  I have
+% no idea... 
+
+n=length(data);
+x_bar=mean(data);
+std_dev_est=sqrt(sum((data-x_bar).^2)/(n-1));
+result=std_dev_est/x_bar;

basics/coeff_of_var_est.m

@@ -0,0 +1,26 @@
+function cmap=bipolar(varargin)
+
+if nargin<1
+  n=256;
+else
+  n=varargin{1};
+end
+
+% n is the number of colors
+clr_pos=[0 1 0];
+clr_neg=[1 0 0];
+
+clr_pos_lab=srgb2lab(clr_pos);
+clr_neg_lab=srgb2lab(clr_neg);
+
+n_half=ceil(n/2);
+scale=linspace(0,1,n_half)';
+cmap_pos_lab=repmat(scale,[1 3]).*repmat(clr_pos_lab,[n_half 1]);
+cmap_neg_lab=repmat(scale,[1 3]).*repmat(clr_neg_lab,[n_half 1]);
+if mod(n,2)==0
+  cmap_lab=[flipud(cmap_neg_lab);cmap_pos_lab];
+else
+  cmap_lab=[flipud(cmap_neg_lab); 0 0 0; cmap_pos_lab];
+end
+cmap=lab2srgb(cmap_lab);
+cmap=max(min(cmap,1),0);

basics/colormaps/bipolar.m

@@ -0,0 +1,17 @@
+function color_map=f(L,x)
+
+% L is a scalar, x is a col vector
+
+% this is the RGB-to-XYZ matrix
+M=[0.412453 0.357580 0.180423;
+   0.212671 0.715160 0.072169;
+   0.019334 0.119193 0.950227];
+
+% this is the code
+n_samples=length(x);
+Y=l2y(L);
+r=1-x;
+g=ungamma_srgb( (Y-M(2,1)*gamma_srgb(1-x)) / M(2,2) );
+b=repmat(0,[n_samples 1]);
+color_map=[r g b];
+

basics/colormaps/blue_off_edge.m

@@ -0,0 +1,17 @@
+function sRGB=f(L,x)
+
+% L is a scalar, x is a col vector
+
+% this is the RGB-to-XYZ matrix
+M=[0.412453 0.357580 0.180423;
+   0.212671 0.715160 0.072169;
+   0.019334 0.119193 0.950227];
+
+% this is the code
+n_samples=length(x);
+Y=l2y(L);
+r=x;
+g=ungamma_srgb( (Y-M(2,1)*gamma_srgb(x)-M(2,3)) / M(2,2) );
+b=repmat(1,[n_samples 1]);
+sRGB=[r g b];
+

basics/colormaps/blue_on_edge.m

@@ -0,0 +1,18 @@
+function cmap=blue_to_yellow(n)
+
+% deal with args
+if nargin<1 || isempty(n)
+  n=256;
+end
+
+% n is the number of colors
+clr_pos=[1 1 0];  % yellow
+clr_neg=[0 0 1];  % blue
+
+clr_one_lab=srgb2lab(clr_pos);
+clr_zero_lab=srgb2lab(clr_neg);
+
+scale=linspace(0,1,n)';
+cmap_lab=interp1([0;1],[clr_zero_lab;clr_one_lab],scale);
+cmap=lab2srgb(cmap_lab);
+cmap=max(min(cmap,1),0);

basics/colormaps/blue_to_yellow.m

@@ -0,0 +1,33 @@
+function clr=f(x)
+
+%  x is a col vector, on [0,1]
+
+% calculates a color based on x, which goes from blue to cyan to green to
+% yellow to red
+
+% the four cases to be handled
+x_frac=mod(6*x,1);
+x_floor=floor(6*x)+1;
+special=(x==1);
+x_frac(special)=1;
+x_floor(special)=6;
+
+% make weights
+n_x=length(x);
+w=zeros(n_x,7);
+for i=1:n_x
+  w(i,x_floor(i)  )=1-x_frac(i);
+  w(i,x_floor(i)+1)=  x_frac(i);
+end
+
+% the colors
+C=[0 0 0;...
+   0 0 1;...
+   0 1 1;...
+   0 1 0;...
+   1 1 0;...
+   1 0 0;...
+   1 1 1];
+
+% do the matrix mult
+clr=w*C;

basics/colormaps/bspectrumw_of_x.m

@@ -0,0 +1,17 @@
+function cmap=f(n_colors)
+
+% generate a non-smooth colormap on a very fine grid
+% calculate the path length, and from that get the phase for each
+% point on the fine grid
+n_samples=4096;
+x=linspace(0,1,n_samples)';
+clr=bspectrumw_of_x(x);
+clr_lab=srgb2lab(clr);
+ds=dist_lab(clr_lab(1:end-1,:),clr_lab(2:end,:));
+s=[0 ; cumsum(ds)];
+phase=s/s(end);  % normalized path length
+
+% make a colormap with inter-color spacings equal
+phase_samples=linspace(0,1,n_colors)';
+x_samples=interp1(phase,x,phase_samples,'linear');
+cmap=bspectrumw_of_x(x_samples);

basics/colormaps/bspectrumw_smooth.m

@@ -0,0 +1,5 @@
+function clr=f(x)
+
+%  x is a col vector, on [0,1]
+
+clr=x*[1 1 1];

basics/colormaps/bw_of_x.m

@@ -0,0 +1,17 @@
+function cmap=f(n_colors)
+
+% generate a non-smooth colormap on a very fine grid
+% calculate the path length, and from that get the phase for each
+% point on the fine grid
+n_samples=4096;
+x=linspace(0,1,n_samples)';
+clr=bw_of_x(x);
+clr_lab=srgb2lab(clr);
+ds=dist_lab(clr_lab(1:end-1,:),clr_lab(2:end,:));
+s=[0 ; cumsum(ds)];
+phase=s/s(end);  % normalized path length
+
+% make a colormap with inter-color spacings equal
+phase_samples=linspace(0,1,n_colors)';
+x_samples=interp1(phase,x,phase_samples,'linear');
+cmap=bw_of_x(x_samples);

basics/colormaps/bw_smooth.m

@@ -0,0 +1,36 @@
+function im=coh2l75_border(z,z_thresh)
+
+% z is mxn of complex, assumed to have mag(z)<=1
+% im is mxnx3, RGB image
+
+if nargin<2 || isempty(z_thresh)
+  z_thresh=0;
+end
+
+n_clr=1000;
+cmap_outer_rgb=l75_border(n_clr);
+cmap_outer_lab=srgb2lab(cmap_outer_rgb);
+z_mag=abs(z);
+z_phase=wrap(angle(z));
+% deal with NaN's, infs, etc.
+messed=~isfinite(z_mag);
+z_mag(messed)=0;
+z_phase(messed)=0;
+% deal with values that might be (hopefully only slightly) out of bounds
+z_mag=min(max(z_mag,0),1);
+im_l=75*z_mag;
+im_l(z_mag<z_thresh)=0;
+ind=round((n_clr-1)*(((z_phase/pi)+1)/2))+1;
+im_lab_r_max=reshape(cmap_outer_lab(ind,:),[size(z) 3]);
+im_ab_r_max=im_lab_r_max(:,:,2:3);
+im_ab=repmat(z_mag,[1 1 2]).*im_ab_r_max;
+%im_ab=im_ab_r_max;
+im_lab=zeros([size(z) 3]);
+im_lab(:,:,1)=im_l;
+im_lab(:,:,2:3)=im_ab;
+
+% convert Lab to sRGB
+im_lab_rows=reshape(im_lab,[size(z,1)*size(z,2) 3]);
+im_rows=lab2srgb(im_lab_rows);
+im=reshape(im_rows,[size(z) 3]);
+im=min(max(im,0),1);  % some things might be a little bit out-of-bounds