 
 
  ******* ./Wavelab850/Biorthogonal/ *******
  Biorthogonal:Contents v850 -- Bi-Orthogonal Wavelet Transform tools
  DownDyadHi_PBS -- Hi-Pass Downsampling operator (periodized,symmetric)
  DownDyadLo_PBS -- Lo-Pass Downsampling operator (periodized,symmetric)
  DownDyad_SBS -- Symmetric Downsampling operator
  dyadpartition -- determine dyadic partition in wavelet transform of 
  extend -- perform various kinds of symmetric extension
  FWT2_PBS -- 2-dimensional wavelet transform (periodized, biorthonormal)
  FWT2_PBS -- 2-dimensional wavelet transform (periodized, biorthonormal)
  FWT2_SBS -- 2-dimensional wavelet transform
  FWT_PB -- Forward Wavelet Transform (periodized, biorthogonal)
  FWT_PBS -- Forward Wavelet Transform (periodized, biorthogonal, symmetric)
  FWT_SBS -- Forward Wavelet Transform (symmetric extension, biorthogonal, symmetric)
  IWT2_PBS -- Inverse 2d Wavelet Transform (periodized, biorthogonal)
  IWT2_PBS -- Inverse 2d Wavelet Transform (periodized, biorthogonal)
  IWT2_SBS -- Inverse 2d Wavelet Transform
  IWT_PB -- Inverse Wavelet Transform (periodized, biorthogonal)
  IWT_PBS -- Inverse Wavelet Transform (periodized, biorthogonal,symmetric)
  iwt_po -- Inverse Wavelet Transform (symmetric extension, biorthogonal, symmetric)
  MakeBSFilter -- Generate Biorthonormal QMF Filter Pairs
  MakeWavelet -- Make periodized orthogonal wavelet
  MirrorSymmFilt -- apply (-1)^t modulation to symmetric filter
  PlotPBSMultiRes -- Multi-Resolution Display of 1-d Wavelet Transform
  symm_aconv -- Symmetric Convolution Tool for Two-Scale Transform
  symm_iconv -- Symmetric Convolution Tool for Two-Scale Transform
  UpDyadHi_PBS -- Hi-Pass Upsampling operator; periodized
  UpDyadLo_PBS -- Lo-Pass Upsampling operator; periodized
  UpDyad_SBS -- Symmetric Upsampling operator
 
 
  ******* ./Wavelab850/Books/WaveTour/WTCh02/ *******
  WTCh02Demo -- Demo Browser for chapter 2.
 
 
  ******* ./Wavelab850/Books/WaveTour/WTCh04/ *******
  ImageAWTnew -- Image of Continuous Wavelet Transform
  MakeChirps -- Make artificial chirp signal
  MakeWindow -- Make artificial Window
  WindowFT -- Window Fourier Transform
  WTCh04Demo -- Demo Browser for chapter 4.
 
 
  ******* ./Wavelab850/Books/WaveTour/WTCh05/ *******
  WTCh05Demo -- Demo Browser for chapter 5.
 
 
  ******* ./Wavelab850/Books/WaveTour/WTCh06/ *******
  aconv2 -- 2-D seperable convolution
  DownShift -- Shift matrix down by one, move bottom line to the top.
  MakeBrownian -- Create Fractional Brownian Signal
  RightShift -- Shift right matrix, the last column is moved to the first.
  WTCh06Demo -- Demo Browser for chapter 6.
 
 
  ******* ./Wavelab850/Books/WaveTour/WTCh07/ *******
  PlotWaveCoeff -- Spike-plot display of wavelet coefficients
  newPlotSpikes -- Plot 1-d signal as baseline with series of spikes
  PlotWaveCoeff -- Spike-plot display of wavelet coefficients
  WTCh07Demo -- Demo Browser for chapter 7.
 
 
  ******* ./Wavelab850/Books/WaveTour/WTCh08/ *******
  WTCh08Demo -- Demo Browser for chapter 8.
 
 
  ******* ./Wavelab850/Books/WaveTour/WTCh09/ *******
  WTCh09Demo -- Demo Browser for chapter 9.
 
 
  ******* ./Wavelab850/Books/WaveTour/WTCh10/ *******
  congauss -- Convolution of a signal with a gaussian filter.
  DiscrDenoise -- Denoising of inverse filtered signal
  DisStatTree -- Build tree with entropy numbers
  'HTE_Compr' -- Thresholding Estimator (white noise).
  'HTE_Discr' -- Thresholding Estimator (colored noise).
  MyPlotBasisTree -- Display basis tree with decorated branch lengths
  wiener -- Wiener filter for blurred and noisy image.
  wienerdec -- Wiener filter for deconvolution.
  Wien_Filt_GWN -- Wiener filtering of the gaussian smoothed signal.
  -- filter f
  WTCh10Demo -- Demo Browser for chapter 10.
 
 
  ******* ./Wavelab850/Browsers/WaveTour/ *******
  WTBrowser -- Mallat's "WaveLet Tour of Signal Processing" Browser
  WTBrowser -- Mallat's "WaveLet Tour of Signal Processing" Browser
 
 
  ******* ./Wavelab850/Browsers/One-D/ *******
  clr_global -- Called by WLBrowser
  Browsers/One-D:Contents v850 -- 1-d Signal Browser
  interaction, such as interactive wavelet thresholding -- changing a
  def_bell -- Called by WLBrowser
  def_data -- Called by WLBrowser
  def_entropy -- Called by WLBrowser
  def_filter -- Called by WLBrowser
  def_nlevel -- Called by WLBrowser
  def_noise -- Called by WLBrowser
  def_nonlin -- Called by WLBrowser
  def_shrinker -- Called by WLBrowser
  def_siglen -- Called by WLBrowser
  def_signal -- Called by WLBrowser
  def_tut -- Called by WLBrowser
  do_action -- Called by WLBrowser
  do_button -- Called by WLBrowser
  do_cpcompress -- Called by WLBrowser
  do_detail -- Called by WLBrowser
  do_global -- Called by WLBrowser
  do_scalegram -- Called by WLBrowser
  do_scalogram -- Called by WLBrowser
  do_shrinkage -- Called by WLBrowser
  do_spectrum -- Called by WLBrowser
  do_transform -- Called by WLBrowser
  do_wpcompress -- Called by WLBrowser
  do_wtcompress -- Called by WLBrowser
  live_cpcompress -- Called by WLBrowser
  live_thresh -- Called by WLBrowser
  live_wpcompress -- Called by WLBrowser
  live_wtcompress -- Called by WLBrowser
  NewShrink -- Called by WLBrowser
  NewZoom -- Called by WLBrowser
  plot_new_data -- Called by WLBrowser
  reset_slider -- Called by WLBrowser
  setfigure -- Called by WLBrowser
  sonfiy -- Called by WLBrowser
  update_plot -- Called by WLBrowser
  WLBrowserIntro -- Called by WLBrowser
  -- Data for Actions menu #4 --- %
  -- Data for Data menu # 1 --- %
  -- Data for Params menu #5 --- %
  -- Data for Signals menu #2 --- %
  -- Data for Xforms menu #3 --- %
  WLBrowser -- WaveLab Browser
  wl_spectrum -- Called by WLBrowser
 
 
  ******* ./Wavelab850/Continuous/ *******
  AWT -- Analytical Wavelet Transform
  BuildSkelMapFast -- Chain together Ridges of Wavelet Transform
  skellen vector of length nchain -- length of skellists
  skelptr vector of length nchain -- pointers to head of chain
  BuildSkelMap -- Chain together Ridges of Wavelet Transform
  skellen vector of length nchain -- length of skellists
  CalcCWTNorms -- Calculate Norm at each Scale of CWT
  CalcCWTPars -- Calculate Underlying Parameters of CWT Structure
  CalcCWTScale -- Calculate Scales and TickMarks for CWT Display
  CalcFracSpectrum -- Calculate Spectrum of Local Scaling Exponents
  CalcGenFracDimen -- Calculate Generalized Fractal dimensions
  CalcMomentGenFun -- Calculate Moment Generating Function
  CalcThermoPartition -- Build Thermodynamic Partition Function
  BuildSkelMapFast Build Map of Skeleton of Wavelet Transform -- Faster
  Continuous:Contents v802 -- Continuous Wavelet Transform tools
  different than in the usual orthogonal wavelet transforms -- they are
  Transforms of 1-d signals, and tasks associated with CWT -- ridge and
  CWT -- Continuous Wavelet Transform
  ExtractRidge -- Pull One Ridge Continuous Wavelet Transform
  skellen vector of length nchain -- length of skellists
  skelptr vector of length nchain -- pointers to heads of chains
  ImageCWT -- Image of Continuous Wavelet Transform
  ImageRWT -- Image of Continuous Wavelet Transform
  ImageWTMM -- Maxima of Continuous Wavelet Transform
  iRWT -- Inverse Real Wavelet Transform
  MM_RWT -- Modulus Maxima of a Real Wavelet Transform
  PlotGenFracDimen -- Plot Generalized Fractal Dimension Function
  PlotMomentGenFun -- Display Moment Generating Function
  PlotMultiSpectrum -- Display Multifractal Spectrum
  PlotRidges -- Log-Log Plot of Ridges in Continuous Wavelet Transform
  skellen vector of length nchain -- length of skellists
  skelptr vector of length nchain -- pointers to heads of chains
  PlotSkelMap -- Display Skeleton of Continuous Wavelet Transform
  skellen vector of length nchain -- length of skellists
  skelptr vector of length nchain -- pointers to heads of chains
  PlotThermoPartition -- Plot Thermodynamic Partition Function
  PruneSkelMap -- Remove Thin Bones in Skeleton of Continuous Wavelet Transform
  skellen vector of length nchain -- length of skellists
  skl vector of length nchain -- length of skellists
  skp vector of length nchain -- pointers to heads of chains
  Ridge_AWT -- Ridges of an Analytic Wavelet Transform
  RWT -- Real Wavelet Transform
  CWTscale -- Complement to Continuous Wavelet Transform
  skellen vector of length nchain -- length of skellists
  SkelMap -- Chain together Ridges of Wavelet Transform
  WTMM -- Wavelet Transform Maximum Modulus
 
 
  ******* ./Wavelab850/Datasets/ *******
  BrowseImages -- Show Images with Point & Click Interface
  Datasets:Contents v850 -- Datasets, Documentation, and Readers
  ImageFig -- Part of BrowseImages demo
  Make2dSignal -- Make artificial 2d signal
  MakeBrownian -- Create Fractional Brownian Signal
  makediag -- Make a diagonal pattern
  MakeFractal -- Create Random or Deterministic Self-Similar signal
  MakeImage -- Make artificial 2d signal
  sig = MakeSignal(Name,n)
% MakeSignal -- Make artificial signal
% Usage
% sig = MakeSignal(Name,n)
% Inputs
% Name string: 'HeaviSine', 'Bumps', 'Blocks',
% 'Doppler', 'Ramp', 'Cusp', 'Sing', 'HiSine',
% 'LoSine', 'LinChirp', 'TwoChirp', 'QuadChirp',
% 'MishMash', 'WernerSorrows' (Heisenberg),
% 'Leopold' (Kronecker), 'Piece-Regular' (Piece-Wise Smooth),
% 'Riemann','HypChirps','LinChirps', 'Chirps', 'Gabor'
% 'sineoneoverx','Cusp2','SmoothCusp','Gaussian'
% 'Piece-Polynomial' (Piece-Wise 3rd degree polynomial)
% n desired signal length
% Outputs
% sig 1-d signal
%
% References
% Various articles of D.L. Donoho and I.M. Johnstone
%
 if nargin > 1,
 t = (1:n) ./n;
 end
 if strcmp(Name,'HeaviSine'),
 sig = 4.*sin(4*pi.*t);
 sig = sig - sign(t - .3) - sign(.72 - t);
 elseif strcmp(Name,'Bumps'),
 pos = [ .1 .13 .15 .23 .25 .40 .44 .65 .76 .78 .81];
 hgt = [ 4 5 3 4 5 4.2 2.1 4.3 3.1 5.1 4.2];
 wth = [.005 .005 .006 .01 .01 .03 .01 .01 .005 .008 .005];
 sig = zeros(size(t));
 for j =1:length(pos)
 sig = sig + hgt(j)./( 1 + abs((t - pos(j))./wth(j))).^4;
 end 
 elseif strcmp(Name,'Blocks'),
 pos = [ .1 .13 .15 .23 .25 .40 .44 .65 .76 .78 .81];
 hgt = [4 (-5) 3 (-4) 5 (-4.2) 2.1 4.3 (-3.1) 2.1 (-4.2)];
 sig = zeros(size(t));
 for j=1:length(pos)
 sig = sig + (1 + sign(t-pos(j))).*(hgt(j)/2) ;
 end
 elseif strcmp(Name,'Doppler'),
 sig = sqrt(t.*(1-t)).*sin((2*pi*1.05) ./(t+.05));
 elseif strcmp(Name,'Ramp'),
 sig = t - (t >= .37);
 elseif strcmp(Name,'Cusp'),
 sig = sqrt(abs(t - .37));
 elseif strcmp(Name,'Sing'),
 k = floor(n * .37);
 sig = 1 ./abs(t - (k+.5)/n);
 elseif strcmp(Name,'HiSine'),
 sig = sin( pi * (n * .6902) .* t);
 elseif strcmp(Name,'LoSine'),
 sig = sin( pi * (n * .3333) .* t);
 elseif strcmp(Name,'LinChirp'),
 sig = sin(pi .* t .* ((n .* .500) .* t));
 elseif strcmp(Name,'TwoChirp'),
 sig = sin(pi .* t .* (n .* t)) + sin((pi/3) .* t .* (n .* t));
 elseif strcmp(Name,'QuadChirp'),
 sig = sin( (pi/3) .* t .* (n .* t.^2));
 elseif strcmp(Name,'MishMash'), % QuadChirp + LinChirp + HiSine
 sig = sin( (pi/3) .* t .* (n .* t.^2)) ;
 sig = sig + sin( pi * (n * .6902) .* t);
 sig = sig + sin(pi .* t .* (n .* .125 .* t));
 elseif strcmp(Name,'WernerSorrows'),
 sig = sin( pi .* t .* (n/2 .* t.^2)) ;
 sig = sig + sin( pi * (n * .6902) .* t);
 sig = sig + sin(pi .* t .* (n .* t));
 pos = [ .1 .13 .15 .23 .25 .40 .44 .65 .76 .78 .81];
 hgt = [ 4 5 3 4 5 4.2 2.1 4.3 3.1 5.1 4.2];
 wth = [.005 .005 .006 .01 .01 .03 .01 .01 .005 .008 .005];
 for j =1:length(pos)
 sig = sig + hgt(j)./( 1 + abs((t - pos(j))./wth(j))).^4;
 end 
 elseif strcmp(Name,'Leopold'),
 sig = (t == floor(.37 * n)/n); % Kronecker
 elseif strcmp(Name,'Riemann'),
 sqn = round(sqrt(n));
 sig = t .* 0; % Riemann's Non-differentiable Function
 sig((1:sqn).^2) = 1. ./ (1:sqn);
 sig = real(ifft(sig));
 elseif strcmp(Name,'HypChirps'), % Hyperbolic Chirps of Mallat's book
 alpha = 15*n*pi/1024;
 beta = 5*n*pi/1024;
 t = (1.001:1:n+.001)./n; 
 f1 = zeros(1,n);
 f2 = zeros(1,n); 
 f1 = sin(alpha./(.8-t)).*(0.1<t).*(t<0.68);
 f2 = sin(beta./(.8-t)).*(0.1<t).*(t<0.75);
 M = round(0.65*n);
 P = floor(M/4);
 enveloppe = ones(1,M); % the rising cutoff function 
 enveloppe(1:P) = (1+sin(-pi/2+((1:P)-ones(1,P))./(P-1)*pi))/2;
 enveloppe(M-P+1:M) = reverse(enveloppe(1:P));
 env = zeros(1,n);
 env(ceil(n/10):M+ceil(n/10)-1) = enveloppe(1:M);
 sig = (f1+f2).*env;
 elseif strcmp(Name,'LinChirps'), % Linear Chirps of Mallat's book
 b = 100*n*pi/1024;
 a = 250*n*pi/1024;
 t = (1:n)./n; 
 A1 = sqrt((t-1/n).*(1-t));
 sig = A1.*(cos((a*(t).^2)) + cos((b*t+a*(t).^2)));
 elseif strcmp(Name,'Chirps'), % Mixture of Chirps of Mallat's book
 t = (1:n)./n.*10.*pi; 
 f1 = cos(t.^2*n/1024);
 a = 30*n/1024;
 t = (1:n)./n.*pi; 
 f2 = cos(a.*(t.^3));
 f2 = reverse(f2);
 ix = (-n:n)./n.*20;
 g = exp(-ix.^2*4*n/1024);
 i1 = (n/2+1:n/2+n);
 i2 = (n/8+1:n/8+n);
 j = (1:n)/n;
 f3 = g(i1).*cos(50.*pi.*j*n/1024);
 f4 = g(i2).*cos(350.*pi.*j*n/1024);
 sig = f1+f2+f3+f4;
 enveloppe = ones(1,n); % the rising cutoff function 
 enveloppe(1:n/8) = (1+sin(-pi/2+((1:n/8)-ones(1,n/8))./(n/8-1)*pi))/2;
 enveloppe(7*n/8+1:n) = reverse(enveloppe(1:n/8));
 sig = sig.*enveloppe;
 elseif strcmp(Name,'Gabor'), % two modulated Gabor functions in 
 % Mallat's book
 N = 512; 
 t = (-N:N)*5/N;
 j = (1:N)./N;
 g = exp(-t.^2*20);
 i1 = (2*N/4+1:2*N/4+N);
 i2 = (N/4+1:N/4+N);
 sig1 = 3*g(i1).*exp(i*N/16.*pi.*j);
 sig2 = 3*g(i2).*exp(i*N/4.*pi.*j);
 sig = sig1+sig2;
 elseif strcmp(Name,'sineoneoverx'), % sin(1/x) in Mallat's book
 N = 1024;
 i1 = (-N+1:N);
 i1(N) = 1/100;
 i1 = i1./(N-1);
 sig = sin(1.5./(i1));
 sig = sig(513:1536);
 elseif strcmp(Name,'Cusp2'),
 N = 64;
 i1 = (1:N)./N;
 x = (1-sqrt(i1)) + i1/2 -.5;
 M = 8*N;
 sig = zeros(1,M);
 sig(M-1.5.*N+1:M-.5*N) = x;
 sig(M-2.5*N+2:M-1.5.*N+1) = reverse(x);
 sig(3*N+1:3*N + N) = .5*ones(1,N);
 elseif strcmp(Name,'SmoothCusp'),
 sig = MakeSignal('Cusp2');
 N = 64;
 M = 8*N;
 t = (1:M)/M;
 sigma = 0.01;
 g = exp(-.5.*(abs(t-.5)./sigma).^2)./sigma./sqrt(2*pi);
 g = fftshift(g);
 sig2 = iconv(g',sig)'/M; 
 elseif strcmp(Name,'Piece-Regular'),
 sig1=-15*MakeSignal('Bumps',n);
 t = (1:fix(n/12)) ./fix(n/12);
 sig2=-exp(4*t);
 t = (1:fix(n/7)) ./fix(n/7);
 sig5=exp(4*t)-exp(4);
 t = (1:fix(n/3)) ./fix(n/3);
 sigma=6/40;
 sig6=-70*exp(-((t-1/2).*(t-1/2))/(2*sigma^2));
 sig(1:fix(n/7))= sig6(1:fix(n/7));
 sig((fix(n/7)+1):fix(n/5))=0.5*sig6((fix(n/7)+1):fix(n/5));
 sig((fix(n/5)+1):fix(n/3))=sig6((fix(n/5)+1):fix(n/3));
 sig((fix(n/3)+1):fix(n/2))=sig1((fix(n/3)+1):fix(n/2));
 sig((fix(n/2)+1):(fix(n/2)+fix(n/12)))=sig2;
 sig((fix(n/2)+2*fix(n/12)):-1:(fix(n/2)+fix(n/12)+1))=sig2;
sig(fix(n/2)+2*fix(n/12)+fix(n/20)+1:(fix(n/2)+2*fix(n/12)+3*fix(n/20)))=...
-ones(1,fix(n/2)+2*fix(n/12)+3*fix(n/20)-fix(n/2)-2*fix(n/12)-fix(n/20))*25;
 k=fix(n/2)+2*fix(n/12)+3*fix(n/20);
 sig((k+1):(k+fix(n/7)))=sig5;
 diff=n-5*fix(n/5);
 sig(5*fix(n/5)+1:n)=sig(diff:-1:1);
 % zero-mean
 bias=sum(sig)/n;
 sig=bias-sig;
 elseif strcmp(Name,'Piece-Polynomial'),
 t = (1:fix(n/5)) ./fix(n/5);
 sig1=20*(t.^3+t.^2+4);
 sig3=40*(2.*t.^3+t) + 100;
 sig2=10.*t.^3 + 45;
 sig4=16*t.^2+8.*t+16;
 sig5=20*(t+4);
 sig6(1:fix(n/10))=ones(1,fix(n/10));
 sig6=sig6*20;
 sig(1:fix(n/5))=sig1;
 sig(2*fix(n/5):-1:(fix(n/5)+1))=sig2;
 sig((2*fix(n/5)+1):3*fix(n/5))=sig3;
 sig((3*fix(n/5)+1):4*fix(n/5))=sig4;
 sig((4*fix(n/5)+1):5*fix(n/5))=sig5(fix(n/5):-1:1);
 diff=n-5*fix(n/5);
 sig(5*fix(n/5)+1:n)=sig(diff:-1:1);
 %sig((fix(n/20)+1):(fix(n/20)+fix(n/10)))=-ones(1,fix(n/10))*20;
 sig((fix(n/20)+1):(fix(n/20)+fix(n/10)))=ones(1,fix(n/10))*10;
 sig((n-fix(n/10)+1):(n+fix(n/20)-fix(n/10)))=ones(1,fix(n/20))*150;
 % zero-mean
 bias=sum(sig)/n;
 sig=sig-bias;
 elseif strcmp(Name,'Gaussian'),
 sig=GWN(n,beta);
 g=zeros(1,n);
 lim=alpha*n;
 mult=pi/(2*alpha*n);
 g(1:lim)=(cos(mult*(1:lim))).^2;
 g((n/2+1):n)=g((n/2):-1:1);
 g = rnshift(g,n/2);
 g=g/norm(g);
 sig=iconv(g,sig);
 else
 disp(sprintf('MakeSignal: I don*t recognize <<%s>>',Name))
 disp('Allowable Names are:')
 disp('HeaviSine'),
 disp('Bumps'),
 disp('Blocks'),
 disp('Doppler'),
 disp('Ramp'),
 disp('Cusp'),
 disp('Crease'),
 disp('Sing'),
 disp('HiSine'),
 disp('LoSine'),
 disp('LinChirp'),
 disp('TwoChirp'),
 disp('QuadChirp'),
 disp('MishMash'),
 disp('WernerSorrows'),
 disp('Leopold'),
 disp('Sing'),
 disp('HiSine'),
 disp('LoSine'),
 disp('LinChirp'),
 disp('TwoChirp'),
 disp('QuadChirp'),
 disp('MishMash'),
 disp('WernerSorrows'),
 disp('Leopold'),
 disp('Riemann'),
 disp('HypChirps'),
 disp('LinChirps'),
 disp('Chirps'),
 disp('sineoneoverx'),
 disp('Cusp2'),
 disp('SmoothCusp'),
 disp('Gabor'),
 disp('Piece-Regular');
 disp('Piece-Polynomial');
 disp('Gaussian');
 end
 
%
% Originally made by David L. Donoho.
% Function has been enhanced.
 
  ReadImage -- Read Image file in 8-bit raw binary format
  ReadSignal -- Read 1-d signal from file in ascii format
 
 
  ******* ./Wavelab850/DeNoising/ *******
  CohWave -- Denoising of 1-d signal with wavelet thresholding.
  DeNoising:Contents v850 -- Tools for DeNoising Signals with White Gaussian Noise
  CPDeNoise -- De-Noising in an adaptively chosen CP basis
  HardThresh -- Apply Hard Threshold 
  HT2 -- Hard Threshold Applied to Image Wavelet Coefficients.
  HT -- Hard Threshold Applied to Wavelet Coefficients.
  HybridThresh -- Adaptive Threshold Selection Using Principle of SURE
  IdealWavDenoise -- Simulation of an Ideal Thresholding Applied to
  InvShrink -- Shrinkage with Exponential Factor applied
  MinMaxThresh -- Minimax Thresholding
  MultiHybrid -- Apply Shrinkage to Wavelet Coefficients
  MultiMAD -- Apply Shrinkage with level-dependent Noise level estimation
  MultiSURE -- Apply Shrinkage to Wavelet Coefficients
  MultiVisu -- Apply sqrt(2log(n)) Thresholding to Wavelet Coefficients
  NormNoise -- Estimates noise level, Normalize signal to noise level 1
  SoftThresh -- Apply Soft Threshold 
  ST2 -- Sofd Threshold Applied to Wavelet Coefficients.
  ST -- Soft Threshold Applied to Wavelet Coefficients.
  SUREThresh -- Adaptive Threshold Selection Using Principle of SURE
  ThreshCP -- Cosine packet decomposition of a noisy signal,
  ThreshWave2 -- Denoising of 2-d image with wavelet thresholding.
  ThreshWave -- Denoising of 1-d signal with wavelet thresholding.
  ThreshWP -- Wavepacket decomposition of a noisy signal,
  TIDenoiseHard2 -- Translation-Invariant Hard Wavelet Thresholding of
  TIDenoiseSoft2 -- Translation-Invariant Soft Wavelet Thresholding of
  ValSUREThresh -- Adaptive Threshold Selection Using Principle of SURE
  VisuThresh -- Visually calibrated Adaptive Smoothing
  WaveShrink -- Soft Threshold Shrinkage Applied to Wavelet Coefficients
  Wiener -- Wiener filter for images.
  Wiener -- Wiener filter.
  WPDeNoise -- De-Noising in an adaptively chosen WP basis
 
 
  ******* ./Wavelab850/Documentation/ *******
  ADDINGNEWFEATURES -- To Add New Features to WaveLab
  MyOwnCleanup -- clears all globals created by the demo
  MyOwnDemo -- starts the Demonstration, invokes Choices
  MyOwnFig -- called from Choices
  MyOwnInit -- sets up data structures
  MyOwnIntro -- help file, explains contents of directories
  one script and then used in later scripts -- these variables
  BUGREPORT -- Report Problem with WaveLab
  Documentation:Contents v850 -- Directory of Documentation Tools
  COPYING -- WaveLab Copying Permissions
  2-d Packet Table There are no 2-d packet tables in WaveLab -- they would
  DATASTRUCTURES -- basic data structures in WaveLab
  FEEDBACK -- Encouraging User Feedback about WaveLab
  0. Install the Software -- see INSTALLATION.m
  GETTINGSTARTED -- ideas for getting started with WaveLab
  INSTALLATION -- Installation of WaveLab
  LIMITATIONS -- WaveLab known limitations
  PAYMENT -- No Charge for WaveLab Software
  Books -- in alphabetical order
  READING -- Sources for further reading about wavelets
  REGISTRATION -- WaveLab Registration
  SUPPORT -- WaveLab Support
  THANKS -- Thanks to contributors
  WARRANTY -- No Warranty on WaveLab software
 
 
  ******* ./Wavelab850/FastAlgorithms/ *******
  FastAlgorithms:Contents v850 -- Tools for Fast Matrix Multiplication in Wavelet Basis
  FNWT_PO -- Nonstandard Wavelet Transform (Periodized, Orthogonal)
  FSFT_PO -- Convert Matrix to Standard Form (Periodized,Orthogonal)
  INWT_PO -- Inverse Nonstandard Wavelet Transform (Periodized, Orthogonal)
  ISFT_PO -- Convert Standard Form to Matrix Form (Periodized,Orthogonal)
  Mat2SparseNSForm -- Put Matrix into Sparse Nonstandard Form
  Mat2SparseSForm -- Put Matrix into Sparse Standard Form
  ndyad -- Index dyad of nonstandard wavelet transform
  NonStdWaveMult -- Nonstandard Wavelet Matrix Multiplication
  StdWaveMult -- Standard Wavelet Matrix Multiplication
  StretchMatrix -- Stretch matrix into BCR nonstandard form
 
 
  ******* ./Wavelab850/Fractals/ *******
  FracPartition -- Fractal Partition Function based on wavelet modulus maxima
  FracScalExp -- Calculate Moment Generating Function
  FracSingSpect -- Calculate Spectrum of Local Scaling Exponents
  Fractals:Contents v850 -- Fractal Analysis
  FracPartition -- Fractal Partition Function based on wavelet modulus maxima
  FracScalExp -- Calculate Moment Generating Function
  FracSingSpect -- Calculate Spectrum of Local Scaling Exponents
 
 
  ******* ./Wavelab850/Interpolating/ *******
  AIDyadDown -- Average-Interpolating DownSampling Operator
  AIDyadUp -- Average-Interpolating UpSampling Operator
  AIRefine2d -- 2-d Refinement based on Average-Interpolation
  AIRefine -- Average-Interpolating Refinement Scheme
  Interpolating:Contents v850 -- Interpolating Refinements and Wavelet Transforms
  DDDyadDown -- Deslauriers-Dubuc Interpolating DownSampling Operator
  DDDyadUp -- Interpolating UpSampling Operator
  DDRefine -- Deslauriers-Dubuc Refinement Scheme
  FWT_AI -- Average-Interpolating wavelet transform
  FWT_DD -- Interpolating Wavelet Transform 
  HCoarsen2d -- 2-d Coarsening operator based on Haar Pyramid
  HRefine2d -- 2-d Refinement based on Haar Pyramid
  IWT_AI -- Inverse transform, average-interpolating wavelets
  IWT_DD -- Inverse transform Deslauriers-Dubuc wavelets
  MakeAIBdryFilter -- Edge filters for Average-Interpolating Wavelets
  MakeAIFilter -- Filters for Average-Interpolating Wavelets
  MakeDDBdryFilter -- Edge filters, Interpolating (Deslauriers-Dubuc) Refinement
  MakeDDFilter -- Filter for interpolating (Deslauriers-Dubuc) Refinement
  PlotAIMultiRes -- Average-Interpolating MRA
  PlotDDMultiRes -- Interpolating MRA
 
 
  ******* ./Wavelab850/Invariant/ *******
  BestShift -- Best Circulant Shift of Signal
  CalcStatPath -- build path-tree filled with entropy numbers
  'Entropy' -- Coifman-Wickerhauser
  'Log' -- sum log(abs(th_i))
  'l^p' -- sum th_i^p, 0 < p < 2, p = par
  'N(eps)' -- #>= eps, eps = par
  'Risk' -- sum min(th_i^2,eps^2), eps=par
  'SURE' -- SURE(Thresholding), thresh = par
  Invariant:Contents v850 -- Invariant Wavelet Transform tools
  DisplayDWT -- Display Dyadic Wavelet Transform by Scale 
  DisplayDWT -- Display Dyadic Wavelet Transform by Scale 
  FWT_TI -- 2-D translation invariant forward wavelet transform
  FWT_ATrou -- Fast Dyadic Wavelet Transform (periodized, orthogonal)
  FWT_Stat -- stationary wavelet transform
  FWT_TI -- translation invariant forward wavelet transform
  IWT2_TI -- Invert 2-d translation invariant wavelet transform
  IWT_ATrou -- Inverse Dyadic Wavelet Transform
  IWT_Stat -- Invert Stationary Wavelet Transform
  IWT_TI -- Invert translation invariant wavelet transform
  MakeATrouFilter -- Generate Biorthonormal Quadratic Spline Filter Pair
  MakeDyadFilter -- Generate Biorthonormal Quadratic Spline Filter Pair
  PackShiftCoeff -- stuff shifted coefficients into packet table
  PlotColumns -- Plot each Column of input matrix
  PlotStatTable -- Plot Stationary Wavelet Transform
  Stat2TI -- Convert Stationary Wavelet Transform to Translation-Invariant Transform 
  TI2Stat -- Convert Translation-Invariant Transform to Stationary Wavelet Transform
  UnpackShiftCoeff -- get wavelet coeff. of circulant shift
 
 
  ******* ./Wavelab850/Median/ *******
  Median:Contents v.850 -- Median Interpolating Pyramid Transform Tools
  FMIPT -- Forward Median Interpolationg Pyramid Transform
  IMIPT -- Inverse Median Transform
  MedITriadDown -- Med Triad Down
  MedITriadUp -- Med Triad Up
  PlotMRA_MT -- Multi-Resolution Display of 1-d Median Transform
  PlotMTCoeff -- Spike-plot display of Median Transform coefficients
  quadmedian -- find the median value of a quadratic polynomial on an interval
  QuadMedRef -- Quadratic Interpolation given median information
  triad -- Index entire j-th triad of 1-d wavelet xform
 
 
  ******* ./Wavelab850/Median/HigherDegree/ *******
  BlockMedian -- Compute median of a polynomial over an interval
  Median:Contents v.850 -- Median Interpolating Pyramid Transform Tools
  MedianInterp -- Perform Polynomial Interpolation of Block Medians
 
 
  ******* ./Wavelab850/Meyer/ *******
  CoarseMeyerCoeff -- Resume coefficients, coarse level
  CoarseMeyerProj -- Invert Meyer Transform, coarse level C
  CombineCoeff -- Combine local trig. coeff's into wavelet coeff's
  different than in the usual orthogonal wavelet transforms -- they are
  frequency domain -- rather than filtering and decimation in the time
  Meyer:Contents v850 -- Periodic Meyer Wavelet Transform Tools
  dct_iii -- Discrete Cosine Transform of Type III
  dct_ii -- Discrete Cosine Transform of Type II
  DetailMeyerCoeff -- Detail coefficients, level j, 3 <= j <= J-2
  DetailMeyerProj -- Invert Meyer Transform, level j, 3 <= j <= J-2
  dst_iii -- Discrete Sine Transform of Type III
  dst_ii -- Discrete Sine Transform of Type II
  dst_i -- Discrete Sine Transform, Type I
  ExtendProj -- Extend a projection to all of the integers -n/2+1 -> n/2
  FineMeyerCoeff -- Meyer Coefficients, finest level
  FineMeyerProj -- Invert Meyer Transform, finest level J-1
  FoldMeyer -- Fold a vector onto itself using a specified window
  FTWT2_YN -- 2-dimensional Tensor Product Meyer Wavelet Transform
  FWT2_YM -- 2d Forward Wavelet Transform (periodized Meyer Wavelet)
  FWT_YM -- Forward Wavelet Transform (periodized Meyer Wavelet)
  ITWT2_YM -- Inverse 2d Tensor Meyer Wavelet Transform
  IWT2_YM -- 2d Inverse Wavelet Transform (periodized Meyer Wavelet)
  IWT_YM -- Inverse Wavelet Transform (periodized Meyer Wavelet)
  PlotYMMultiRes -- Plot MRA for Meyer Wavelet
  QuasiDCT -- Nearly Discrete Cosine Transform of Type I.
  QuasiDST -- Nearly Discrete Sine Transform of Type I.
  SeparateCoeff -- Separate wavelet coeff. into local trig. coeff.
  UnfoldMeyer -- Unfold a vector using a specified window.
  WindowMeyer -- auxiliary window function for Meyer wavelets.
 
 
  ******* ./Wavelab850/Orthogonal/ *******
  aconv -- Convolution Tool for Two-Scale Transform
  CDJV -- Cohen, Daubechies, Jawerth and Vial, 1992.
  CDJVDyadDown -- Downsampling with Boundary Correction
  CDJV -- Cohen, Daubechies, Jawerth and Vial, 1992.
  CDJVDyadUp -- UpSampling with CDJV Boundary Correction
  Orthogonal:Contents v850 -- Orthogonal Wavelet Transform Tools
  ContourMultiRes -- Multi-Resolution Mesh Display of 1-d Wavelet Transform
  DisplayMultiRes -- Multi-Resolution Mesh Display of 1-d Wavelet Transform
  DisplayWaveCoeff -- Display wavelet coefficients
  DownDyadHi -- Hi-Pass Downsampling operator (periodized)
  DownDyadLo -- Lo-Pass Downsampling operator (periodized)
  dyad2ix -- Convert wavelet indexing into linear indexing
  dyadlength -- Find length and dyadic length of array
  dyad -- Index entire j-th dyad of 1-d wavelet xform
  FTWT2_PO -- 2-d tensor wavelet transform (periodized, orthogonal)
  FWT2_PO -- 2-d MRA wavelet transform (periodized, orthogonal)
  FWT_CDJV -- Forward Wavelet Transform (boundary-corrected)
  FWT_IO -- Forward Wavelet Transform (boundary-corrected)
  FWT_PO -- Forward Wavelet Transform (periodized, orthogonal)
  iconv -- Convolution Tool for Two-Scale Transform
  IWT2_PO -- Inverse 2-d tensor wavelet transform (periodized, orthogonal)
  IWT2_PO -- Inverse 2-d MRA wavelet transform (periodized, orthogonal)
  IWT_CDJV -- Inverse Wavelet Transform (boundary corrected)
  IWT_IO -- Inverse Wavelet Transform (boundary corrected)
  IWT_PO -- Inverse Wavelet Transform (periodized, orthogonal)
  lshift -- Circular left shift of 1-d signal
  Make2dWavelet -- Make 2-d Wavelet
  MakeCDJVFilter -- Set up filters for CDJV Wavelet Transform
  MakeOBFilter -- Makes Orthogonal Boundary conjugate mirror filters 
  MakeONFilter -- Generate Orthonormal QMF Filter for Wavelet Transform
  MakeWavelet -- Make periodized orthogonal wavelet
  MirrorFilt -- Apply (-1)^t modulation
  PlotMultiRes -- Multi-Resolution Display of 1-d Wavelet Transform
  PlotSpikes -- Plot 1-d signal as baseline with series of spikes
  PlotWaveCoeff -- Spike-plot display of wavelet coefficients
  quad2ix -- Translate 2-d wavelet indexing into linear indexing
  quadlength -- Find length and dyadic length of square matrix
  reverse -- Reverse order of elements in 1-d signal
  rshift -- Circular right shift of 1-d signal
  UpDyadHi -- Hi-Pass Upsampling operator; periodized
  UpDyadLo -- Lo-Pass Upsampling operator; periodized
  UpSample -- Upsampling operator
 
 
  ******* ./Wavelab850/Packets/One-D/ *******
  CPTour -- Cosine packet decomposition & Best Basis analysis
  AutoWPTour -- Wavelet packet decomposition & Best Basis analysis
  BestBasis -- Coifman-Wickerhauser Best-Basis Algorithm
  CalcStatTree -- Build tree with entropy numbers
  'Entropy' -- Coifman-Wickerhauser
  'Log' -- sum log(abs(th_i))
  'l^p' -- sum |th_i|^p, 0 < p < 2, p = par
  'N(eps)' -- #>= eps, eps = par
  'Risk' -- sum min(th_i^2,eps^2), eps=par
  'Sum' -- sum th_i
  'SURE' -- SURE(Thresholding), thresh = par
  CalcTreeHeight -- Measure the total height of a stat-tree
  CalcWPLocation -- Calculate location of wavelet packet entry
  CompareStdBases -- Calculate entropy of some standard bases
  Packets/One-D:Contents v850 -- One-d Wavelet- and Cosine- Packet Tools
  CPAnalysis -- Dyadic table of all Cosine Packet Coefficients
  CPSynthesis -- Synthesize signal from cosine packet table
  CPTour -- Cosine packet decomposition & Best Basis analysis
  dct_iv -- Type (IV) Discrete Cosine Xform
  DrawHeisenberg -- Draw Heisenberg rectangle
  edgefold -- Perform folding projection with (+,-) polarity at EDGES
  edgeunfold -- Undo folding projection with (+,-) polarity at EDGES
  fold -- Folding projection with (+,-) polarity
  FPT_CP -- Fast transform into specific cosine packet basis
  % terminal node -- analyze by dct_iv
  FPT_WP -- Analyze image into specific wavelet packet basis
  ImageGaborPhase -- Time-Frequency Display with congruent rectangles
  ImagePacketPhase -- Partition phase space by rectangular blocks
  ImagePhasePlane -- Partition phase space by rectangular blocks
  ImagePhaseVarious -- Compare four phase plane displays
  IPT_CP -- Synthesize signal from cosine packet coefficients
  IPT_WP -- Synthesize signal from wavelet packet coefficients
  MakeBasis -- Create a basis tree for a standard basis
  MakeCosinePacket -- Make Cosine Packet signal
  MakeONBell -- Make Bell for Orthonormal Local Cosine Analysis
  MakeWaveletPacket -- Make periodized orthogonal wavelet packet
  node -- Tree indexing function
  PackBasisCoeff -- Stuff basis coefficients into packet table
  packet -- Packet table indexing
  PlotBasisTree -- Display basis tree with decorated branch lengths
  PlotCoeffComparison -- Compare coefficients in various bases
  PlotCompressNumbers -- Plot Compression Numbers
  PlotGaborPhase -- Partition phase space with congruent rectangles
  PlotPacketTable -- Display entries in wavelet, cosine packet tables
  PlotPhasePlane -- Plot Time-Frequency content of signal
  PlotPhaseTiling -- partition phase space by rectangular blocks
  PlotPhaseVarious -- Show several different phase plane images
  PlotWavePhase -- Partition phase space by wavelet tiling
  unfold -- Undo folding projection with (+,-) polarity
  UnpackBasisCoeff -- Extract basis coefficients from packet table
  UnpackStatValue -- Unpack the value for given basis from the stat tree
  WPAnalysis -- Dyadic table of all Wavelet Packet coefficients
  WPSynthesis -- Synthesize signal from wavelet packet coefficients
  WPTour -- Wavelet packet decomposition & Best Basis analysis
  WrapAround -- Circular map onto 1:n
  WTWPTour -- Wavelet packet decomposition & Best Basis analysis
 
 
  ******* ./Wavelab850/Packets/Two-D/ *******
  Best2dBasis -- Coifman-Wickerhauser Best-2d-Basis Algorithm
  BuildFBITree -- Build the 2-d Basis Tree for fingerprint compression
  Calc2dPktTable -- Put Wavelet/Cosine Packet Coefficients into 2-d Table 
  Calc2dSQTree -- Put Packet Table Statistics into Quad Tree 
  'Entropy' -- Coifman-Wickerhauser
  'Log' -- sum log |th_i| 
  'l^p' -- sum |th_i|^p, 0 < p < 2, p = par
  'N(eps)' -- #>= eps, eps = par
  'Risk' -- sum min(th_i^2,eps^2), eps=par
  'Sum' -- sum th_i
  'SURE' -- SURE(Thresholding), thresh = par
  Calc2dStatTree -- Put Wavelet/Cosine Packet Statistics into Quad Tree
  'Entropy' -- Coifman-Wickerhauser
  'Log' -- sum log(abs(th_i))
  'l^p' -- sum |th_i|^p, 0 < p < 2, p = EntPar
  'N(eps)' -- #>= eps, eps = EntPar
  'Risk' -- sum min(th_i^2,eps^2), eps = EntPar
  'Sum' -- sum th_i
  'SURE' -- SURE(Thresholding), thresh = EntPar
  Calc2dTreeHeight -- Measure the total height of a stat-quad-tree
  CalcEntropy -- Calculate entropy number of array
  'Entropy' -- Coifman-Wickerhauser
  'Log' -- sum log(|th_i|)
  'l^p' -- sum |th_i|^p, 0 < p < 2, p = par
  'N(eps)' -- #>= eps, eps = par
  'Risk' -- sum min(th_i^2,eps^2), eps=par
  'Sum' -- sum th_i
  'SURE' -- SURE(Thresholding), thresh = par
  Packets/Two-D:Contents v850 -- 2-d Wavelet- and Cosine- Packet Tools
  CP2dTour -- Cosine Packet Analysis in Adaptively Chosen Basis
  CP2dTour -- Cosine Packet Analysis in Adaptively Chosen Basis
  dct2_iv -- 2-dimensional discrete cosine transform (type IV)
  DownQuad -- Split 2-d image into 4 subbands
  FPT2_CP -- Analyze image into 2-d cosine packet coefficients
  scale 0, fold around edges -- to be added in a later version
  % terminal node -- analyze by dct2_iv
  FPT2_WP -- Analyze image into 2-d wavelet packet coefficients
  IPT2_CP -- Synthesize image from 2-d cosine packet coefficients
  scale 0, unfold around edges -- in a later version
  % terminal node -- reconstruct by dct2_iv
  IPT2_WP -- Synthesize image from 2-d wavelet packet coefficients
  % reconstruct -- apply UpQuad operator
  Make2dCosinePacket -- Create 2d Cosine packet
  Make2dWaveletPacket -- Create 2d Wavelet Packet
  Plot2dBasisTree -- Show quad tree of best basis
  Plot2dPartition -- Show partition of Time/Freq plane by best basis
  qnode -- Quad tree indexing
  qpkt2ix -- Translate quad-packet index to linear index
  quadbounds -- x,y bounds of quadlet
  Unpack2dBasisCoeff -- Extract basis coefficients from 2-d packet table
  UpQuad -- Merge four subbands into 2d image
  WP2dTour -- 2d Wavelet Packet Analysis in Adaptively Chosen Basis
  WP2dTour -- 2d Wavelet Packet Analysis in Adaptively Chosen Basis
 
 
  ******* ./Wavelab850/Papers/Adapt/ *******
  AdaptCleanup
% AdaptCleanup -- Clean up Adapt global variables
% Usage
% AdaptCleanup
%

clear global ADAPTFIGNUM Spline_Data xheavi 
clear global Blocks Spline_Pnty yDoppler 
clear global Bumps yblocks 
clear global Doppler qmf ybumps 
clear global HeaviSine t yheavi 
clear global L xDoppler rho 
clear global xblocks N xbumps 
 
 
  AdaptFig(num)
% AdaptFig -- callback for display of individual figures in AdaptDemo
global ADAPTFIGNUM
global N rho qmf L
if 1 <= num & num <=15 ,
 nihil = ~exist('ADAPTFIGNUM');
 if ~nihil,
 empty = isempty(ADAPTFIGNUM);
 else
 empty = 1;
 end
 if ~empty,
 nil = ADAPTFIGNUM == 0;
 else
 nil = 1;
 end
 if nil,
 AdaptInit
 end
 figname = FigNameString('adfig',num);
 clc; home; help(figname);
 eval(figname);
 RegisterPlot(figname,'{R}');
 drawnow;
 ADAPTFIGNUM = num;
elseif num==19, % change default settings
 N = input('sample size for all figures, default 2048, must be power of 2?');
 rho = input('signal-to-noise ratio, default 7?');
 fam = input('family name of QMF, in quotes, default Symmlet?');
 ind = input('index of QMF family, default 8?');
 qmf = MakeONFilter(fam,ind);
 L = input('index of coarsest resolution level, default 6');
 ADAPTFIGNUM = 0;
elseif num==20, %run all figures
 AdaptInit
 for num=1:15,
 figname = FigNameString('adfig',num);
 eval(figname);
 RegisterPlot(figname,'{R}');
 clc; home; help(figname);
 drawnow;
 ADAPTFIGNUM = num;
 end
else
 disp('Bizarre argument in AdaptFig')
 disp(sprintf('num =%g',num))
end


  AdaptInit -- Initialize all global variables for AdaptDemo
global N rho
global t 
global L qmf
global xblocks xbumps xheavi xDoppler
global yblocks ybumps yheavi yDoppler
global Blocks Bumps HeaviSine Doppler
%
 if isempty(N)
 N = 2048; 
 rho = 7;
 qmf = MakeONFilter('Symmlet',8);
 L = 6;
 end
 t = (0:(N-1))./N;
%
 HeaviSine = MakeSignal('HeaviSine',N);
 Bumps = MakeSignal('Bumps',N);
 Blocks = MakeSignal('Blocks',N);
 Doppler = MakeSignal('Doppler',N);
%
 rand('seed',0);
 [xblocks,yblocks] = NoiseMaker(Blocks,rho);
 [xbumps,ybumps] = NoiseMaker(Bumps,rho);
 [xheavi,yheavi] = NoiseMaker(HeaviSine,rho);
 [xDoppler,yDoppler] = NoiseMaker(Doppler,rho);
%
disp('Initializing Parameters:')
disp(sprintf('N = %g',N))
disp(sprintf('rho = %g',rho))
disp(sprintf('L = %g',L))

 
 
  AdaptIntro -- Info about AdaptDemo
%
% The files in this directory can reproduce the figures in the paper
% ``Adapting to Unknown Smoothness via Wavelet Shrinkage''
% by D.L. Donoho and I.M. Johnstone.
%
% This paper is available by www from 
% http://www-stat.stanford.edu/~donoho/Reports/1993/ausws.ps.Z
%
% The paper proposes soft thresholding of wavelet coefficients
% with an adaptively-chosen threshold. The threshold is chosen
% using the SURE (Stein's unbiased estimate of risk) for the
% Soft Thresholding Estimator. The paper shows that this has
% various optimality properties for mean-squared error estimation.
%

 
 
  adfig01 -- Adapt Figure 01: Four Spatially Inhomogenoeus Signals
%
% These four signals are caricatures of signals arising
% in scientific signal processing.
%
% Blocks is a piecewise constant signal, a caricature of a scan line
% through a 2-d image depicting an object with several edges.
%
% Bumps is a very peaky signal, a caricature of an NMR Spectrum.
% 
% HeaviSine is a sinusoid with two clicks.
%
% Doppler is an attenuating sinusoid with time-varying frequency.
%
% These objects are interesting to us because they portray significant
% spatial inhomogeneity.
%
global t 
global xblocks xbumps xheavi xDoppler
global yblocks ybumps yheavi yDoppler
global Blocks Bumps HeaviSine Doppler
%
 %clf;
 versaplot(221,t,xblocks, [],' 1 (a) Blocks ' ,[],[])
 versaplot(222,t,xbumps, [],' 1 (b) Bumps ' ,[],[])
 versaplot(223,t,xheavi, [],' 1 (c) HeaviSine ',[],[])
 versaplot(224,t,xDoppler,[],' 1 (d) Doppler ' ,[],[])

 
 
  adfig02 -- Adapt Figure 02: Noisy Versions of Four Signals
%
% The four objects of Figure 1 with white noise superposed.
% The noise is normally distributed with variance 1.
%
% 

global N L t 
global yblocks ybumps yheavi yDoppler
%
%clf;
versaplot(221,t,yblocks, [],' 2 (a) Noisy Blocks ' ,[],[])
versaplot(222,t,ybumps, [],' 2 (b) Noisy Bumps ' ,[],[])
versaplot(223,t,yheavi, [],' 2 (c) Noisy HeaviSine ',[],[])
versaplot(224,t,yDoppler,[],' 2 (d) Noisy Doppler ' ,[],[])

 
 
  adfig03 -- Adapt Figure 03: Wavelet Shrinkage of the Four Noisy Signals
%
% Here we aapply a specific thresholding rule to the
% four noisy signals depicted in Figure 2.
%
% The procedure for DeNoising:
% 1. Transform to Wavelet Domain, Using Nearly Symmetric Wavelet
% with 8 vanishing moments.
% 2. Apply a soft thresholding nonlinearity, with threshold selected
% by the Stein's Unbiased Risk Estimate (SURE) in the interval
% [0,sqrt{2 log(n)}]
% 3. Transform back to the signal domain.
% 
% The reconstructions suppress the noise, while preserving the sharp structure
% in the neighborhood of the highly-variable spatial components.
%
global L qmf
global yblocks ybumps yheavi yDoppler
global t
%

 [xhat,xw] = WaveShrink(yblocks,'Hybrid',L,qmf);
 versaplot(221,t,xhat,[],' 3 (a) SUREShrink[Blocks]',[],[])
%
 [xhat,xw] = WaveShrink(ybumps,'Hybrid',L,qmf);
 versaplot(222,t,xhat,[],' 3 (b) SUREShrink[Bumps]',[],[])
%
 [xhat,xw] = WaveShrink(yheavi,'Hybrid',L,qmf);
 versaplot(223,t,xhat,[],' 3 (c) SUREShrink[HeaviSine]',[],[])
%
 [xhat,xw] = WaveShrink(yDoppler,'Hybrid',L,qmf);
 versaplot(224,t,xhat,[],' 3 (d) SUREShrink[Doppler]',[],[])

 
 
 
  adfig04 -- Adapt Figure 04: LPJS reconstruction
%
% The James-Stein Estimate is applied to Dyadic Fourier
% Coronae. This shrinker can be described as follows:
% segment the Fourier Domain into blocks of size 2^j, 
% estimate the signal strength in each block, and apply that
% scalar shrinkage to each block which is appropriate for
% the estimated signal strength. The first block begins at 
% L = 5.
%
% Despite the promising connection with classical multivariate
% shrinkage, the estimator has a disappointing visual appearance.
%
global yblocks ybumps yheavi yDoppler
global L t
%
%clf;
%
[xhat,xw] = CoronaJS(yblocks,L);
versaplot(221,t,xhat,[],' 4 (a) LPJS[Blocks]',[],[])
[xhat,xw] = CoronaJS(ybumps,L);
versaplot(222,t,xhat,[],' 4 (b) LPJS[Bumps]',[],[])
[xhat,xw] = CoronaJS(yheavi,L);
versaplot(223,t,xhat,[],' 4 (c) LPJS[HeaviSine]',[],[])
[xhat,xw] = CoronaJS(yDoppler,L);
versaplot(224,t,xhat,[],' 4 (d) LPJS[Doppler]',[],[])

 
 
 
  adfig05 -- Adapt Figure 05: A Display of all four wavelets, j=6, k=32;
%
% The display shows four mother wavelets from the Haar,
% original Daubechies, Coiflet, and Nearly Symmetric Families.
%
 xw = zeros(1,2048);
 xw(dyad2ix(6,32)) = 1;
 HaarQMF = MakeONFilter('Haar');
 HaarWave = IWT_PO(xw,5,HaarQMF);
 %HaarWave = MakeWavelet(6,32,'Haar',1,'Mother',2048);
 D4QMF = MakeONFilter('Daubechies',4);
 DaubWave = IWT_PO(xw,5,D4QMF);
 %DaubWave = MakeWavelet(6,32,'Daubechies',4,'Mother',2048);
 C3QMF = MakeONFilter('Coiflet',3);
 CoifWave = IWT_PO(xw,5,C3QMF);
 %CoifWave = MakeWavelet(6,32,'Coiflet',3,'Mother',2048);
 S8QMF = MakeONFilter('Symmlet',8);
 SymmWave = IWT_PO(xw,5,S8QMF);
 %SymmWave = MakeWavelet(6,32,'Symmlet',8,'Mother',2048);
%
% 2 * 2 array of plots
%
 %clf;
 tt = (.5:2047.5)/2048;
 versaplot(221,tt,HaarWave,[],'5 (a) Haar Wavelet', [.4 .6 -.3 +.3],[])
 versaplot(222,tt,DaubWave,[],'5 (b) Daubechies D4 Wavelet',[.4 .6 -.3 +.3],[])
 versaplot(223,tt,CoifWave,[],'5 (c) Daubechies C3 Coiflet',[.4 .6 -.3 +.3],[])
 versaplot(224,tt,SymmWave,[],'5 (d) Daubechies S8 Symmlet',[.4 .6 -.3 +.3],[])
 
 
  adfig06 -- Adapt Figure 06: MultiResolution Analysis of the Four Functions
% 
% Mallat's MRA groups together all the terms in the wavelet 
% decomposition which occupy the same scale (with perhaps
% different locations).
%
% Here we display the MRA of the four functions from Figure 1,
% using the Daubechies Nearly Symmetric S8 Wavelet.
%
global Blocks Bumps HeaviSine Doppler
S8QMF = MakeONFilter('Symmlet',8);
L = 6; %clf;
% 
 subplot(221);
 xw = FWT_PO(Blocks,L,S8QMF);
 IMJPlotMultiRes(xw,L,.2,S8QMF);
 title(' 6 (a) Blocks ')
%
 subplot(222);
 xw = FWT_PO(Bumps,L,S8QMF);
 IMJPlotMultiRes(xw,L,.2,S8QMF);
 title(' 6 (b) Bumps ')
%
 subplot(223)
 xw = FWT_PO(HeaviSine,L,S8QMF);
 IMJPlotMultiRes(xw,L,.2,S8QMF);
 title(' 6 (c) HeaviSine ')
%
 subplot(224)
 xw = FWT_PO(Doppler,L,S8QMF);
 IMJPlotMultiRes(xw,L,.2,S8QMF);
 title(' 6 (d) Doppler ')
 
 
 
  adfig07 -- Adapt Figure 07: MultiResolution Analysis of Blocks and HeaviSine
% Using Haar & Db4
%
% Mallat's MRA groups together all the terms in the wavelet 
% decomposition which occupy the same scale (with perhaps
% different locations).
%
% Here we display the MRA of the four functions from Figure 1,
% using the Daubechies Nearly Symmetric S8 Wavelet. Each curve
% represents the contributions arising from a single scale.
%

global L
global Blocks Bumps HeaviSine Doppler
%
HaarQMF = MakeONFilter('Haar');
%
 %clf; 
 subplot(221);
 xw = FWT_PO(Blocks,L,HaarQMF);
 IMJPlotMultiRes(xw,L,.2,HaarQMF);
 title(' 7 (a) Haar[Blocks] ')
%
 subplot(222);
 xw = FWT_PO(HeaviSine,L,HaarQMF);
 IMJPlotMultiRes(xw,L,.2,HaarQMF);
 title(' 7 (b) Haar[HeaviSine] ')
%
Db4QMF = MakeONFilter('Daubechies',4);
%
 subplot(223)
 xw = FWT_PO(Blocks,L,Db4QMF);
 IMJPlotMultiRes(xw,L,.2,Db4QMF);
 title(' 7 (c) Db4[Blocks] ')
%
 subplot(224)
 xw = FWT_PO(HeaviSine,L,Db4QMF);
 IMJPlotMultiRes(xw,L,.2,Db4QMF);
 title(' 7 (d) Db4[HeaviSine] ')

 
 
  adfig08 -- Adapt Figure 08: S8 Wavelet Coefficients of the Four Functions
%
% Here we display the wavelet coefficients of the four functions,
% stratified by scale and location. The significant coefficients
% occur in locations of significant spatial variability. For
% locally smooth functions, the coefficients tend to be small
% at fine scales.
%
global L
global Blocks Bumps HeaviSine Doppler
%
S8QMF = MakeONFilter('Symmlet',8);
%clf;
%
 subplot(221)
 wx = FWT_PO(Blocks,L,S8QMF);
 IMJPlotWaveCoeff(wx,L,0.);
 title(' 8 (a) Blocks ')
%
 subplot(222)
 wx = FWT_PO(Bumps,L,S8QMF);
 IMJPlotWaveCoeff(wx,L,0.);
 title(' 8 (b) Bumps ')
%
 subplot(223)
 wx = FWT_PO(HeaviSine,L,S8QMF);
 IMJPlotWaveCoeff(wx,L,0.);
 title(' 8 (c) HeaviSine ')
%
 subplot(224)
 wx = FWT_PO(Doppler,L,S8QMF);
 IMJPlotWaveCoeff(wx,L,0.);
 title(' 8 (d) Doppler ')
%axis
 
 
  adfig09 -- Adapt Figure 09: Haar Coefficients of the Four Functions
%
% Here we display the wavelet coefficients of the four functions,
% stratified by scale and location. The significant coefficients
% occur in locations of significant spatial variability. For
% locally smooth functions, the coefficients tend to be small
% at fine scales. The coefficients do not decay with scale as 
% rapidly as the coefficients in an S8 wavelet expansion. This
% is particularly evident by comparing (Panel d) of this figure
% with (Panel d) of the previous figure.
%

global Blocks Bumps HeaviSine Doppler
HaarQMF = MakeONFilter('Haar');
HL = 3;
%clf;
% 
 subplot(221)
 xw = FWT_PO(Blocks,HL,HaarQMF);
 IMJPlotWaveCoeff(xw,HL,0.);
 title(' 9 (a) Haar[Blocks] ')
%
 subplot(222)
 xw = FWT_PO(Bumps,HL,HaarQMF);
 IMJPlotWaveCoeff(xw,L,0.);
 title(' 9 (b) Haar[Bumps] ')
%
 subplot(223)
 xw = FWT_PO(HeaviSine,HL,HaarQMF);
 IMJPlotWaveCoeff(xw,HL,0.);
 title(' 9 (c) Haar[HeaviSine] ')
%
 subplot(224)
 xw = FWT_PO(Doppler,HL,HaarQMF);
 IMJPlotWaveCoeff(xw,HL,0.);
 title(' 9 (d) Haar[Doppler] ')
%
 
 
  adfig10 -- Adapt Figure 10: Wavelet Shrinkage of object yBlocks in Haar Basis
%
% (Panel a) depicts the noisy object yBlocks, its Haar transform (Panel c),
% wavelet shrinkage reconstruction using the Haar wavelet (Panel b), and 
% the Haar Transform of the reconstruction (Panel d).
%
% The viewer is supposed to notice that in the Haar domain, the
% noise is spread out among all coefficients, while the signal is
% concentrated in only a few coefficients. Hence thresholding mostly
% affects the noise without disturbing the signal.
%
global t
global yblocks
%
HQMF = MakeONFilter('Haar');
HL = 3; 
[xhat,xw] = WaveShrink(yblocks,'Hybrid',HL,HQMF);
yw = FWT_PO(yblocks,HL,HQMF);
%clf;
ax = [ 0 1 (-10) 30];
%
 versaplot(221,t,yblocks,[],' 10 (a) Noisy Data y' ,ax,[]);
 versaplot(222,t,xhat ,[],' 10 (b) HaarShrink Reconstruction',ax,[]);
 subplot(223); IMJPlotWaveCoeff(yw,HL,.025); title(' 10 (c) Haar[y]')
 subplot(224); IMJPlotWaveCoeff(xw,HL,.025); title(' 10 (d) Haar[Reconstruction]')
 
 
 
  adfig11 -- Adapt Figure 11: Simulate the SURE threshold estimation procedure
%
% (Panel a) A 128-long vector, where only the first 16 
% entries are nonzero, contaminated with Normal(0,1) noise.
% (Panel b) Reconstruction with soft thresholding chosen
% by minimizing the SURE criterion. 
% (Panel c) The SURE criterion for soft thresholding, as a
% function of the threshold. (Panel d) The actual loss of the 
% thresholding rule, as a function of thethreshold.

%
% 1. Create a normal vector with sparse means
%
 theta = zeros(1,128);
 theta(1:16) = 4 .* ones(1,16);
 rand('seed',0);
 x = theta + WhiteNoise(theta);
 %clf;
 n = 128;
%
% 2. Calculate the SURE curve
%
 a = sort(abs(x)).^2 ;
 b = cumsum(a);
 n = length(x);
 c = linspace(n-1,0,n);
 s = b+c.*a;
 risk = (n - ( 2 .* (1:n )) + s)/n;
 [guess,ibest] = min(risk);
 lambda = sqrt( a(ibest) );
 res = (abs(x) - lambda);
 res = (res + abs(res))/2;
 xhat = sign(x).*res;
%
% 3. Calculate the Loss curve
%
 loss = zeros(size(x));
 lam = sqrt(a);
 for i=1:128
 loss(i) = sum( (SoftThresh(x,lam(i)) - theta).^2 ) ./ 128;
 end
%
% 4. Display
%
 versaplot(221,1:n, x ,[],'11(a) raw data',[0 128 -4 6],[])
 versaplot(222,1:n, xhat,[],'11(b) estimate',[0 128 -4 6],[])
 ax = [ 0 max(sqrt(a)) 0 max(risk) ];
 versaplot(223,lam, risk,[],'11(c) estimated risk versus lambda',ax,[])
 versaplot(224,sqrt(a), loss,[],'11(d) loss versus lambda',ax,[])

 
 
  adfig12 -- Adapt Figure 12: Simulating the behavior of various threshold rules
%
% Here we compare the behavior of several schemes for threshold
% choice, when the noise is Normal(0,1). Several Different degrees
% of sparsity are studied.
%
% (Panel a) Choice by minimizing pure SURE criterion.
% (Panel b) Choice by fixed sqrt{2 log(n)} criterion.
% (Panel c) Choice by Hybrid of (a) and (b).
% (Panel d) Choice by Hybrid on random half-samples.
%
NProb=1024; M=25;
rand('seed',0);
%
 EvalSure
 EvalFix
 EvalHybrid
 EvalRigor
%
ax = [0 .3 0 1];
%clf;
versaplot(221,SureRecords(:,2),SureRecords(:,3).^2,'o',...
 '12 (a) SURE Threshold',ax,[]);
xlabel('Epsilon')
ylabel('root MSE')
versaplot(222,FixRecords(:,2),FixRecords(:,3).^2,'.',...
 '12 (b) Fixed Threshold',ax,[]);
xlabel('Epsilon')
ylabel('root MSE')
versaplot(223,HybridRecords(:,2),HybridRecords(:,3).^2,'*',...
 '12 (c) Hybrid(*) Threshold',ax,[]);
xlabel('Epsilon')
ylabel('root MSE')
versaplot(224,PlusRecords(:,2),PlusRecords(:,3).^2,'+',...
 '12 (d) Hybrid(+) Threshold',ax,[]);
xlabel('Epsilon')
ylabel('root MSE')
 
 

  adfig13 -- Adapt Figure 13: WaveJS Reconstructions from Noisy Data; S8 Wavelet
%
% A classical James-Stein shrinker applied in the wavelet domain
% to different resolution levels independently of one another.
%
% The results are noisier-looking than the results of the 
% nonlinear soft thresholding of wavelet coefficients.
%
global L qmf 
global yblocks ybumps yheavi yDoppler
global t
%
 %clf;
 [xhat,xw] = WaveJS(yblocks,L,qmf);
 versaplot(221,t,xhat,[],' 13 (a) WaveJS[Blocks]',[],[])
%
 [xhat,xw] = WaveJS(ybumps,L,qmf);
 versaplot(222,t,xhat,[],' 13 (b) WaveJS[Bumps]',[],[])
%
 [xhat,xw] = WaveJS(yheavi,L,qmf);
 versaplot(223,t,xhat,[],' 13 (c) WaveJS[HeaviSine]',[],[])
%
 [xhat,xw] = WaveJS(yDoppler,L,qmf);
 versaplot(224,t,xhat,[],' 13 (d) WaveJS[Doppler]',[],[])
 
 
 
  adfig14 -- Adapt Figure 14: AutoTrunc: Fourier Truncated
% Reconstructions from Noisy Data;
%
% Here we depict a reconstruction using truncated Fourier expansions,
% in which only the low-order fourier components of the noisy signal
% are retained in reconstruction. The cut-off point for reconstruction
% in the frequency domain is automatically chosen; it is the empirical 
% minimizer of the Stein Unbiased Estimate of Risk among truncation rules.
%
global yblocks ybumps yheavi yDoppler
global t
%
 %clf;
 xhat = AutoTrunc(yblocks,1.);
 versaplot(221,t,xhat,[],' 14 (a) AutoTrunc[Blocks]',[],[])
%
 xhat = AutoTrunc(ybumps,1.);
 versaplot(222,t,xhat,[],' 14 (b) AutoTrunc[Bumps]',[],[])
%
 xhat = AutoTrunc(yheavi,1.);
 versaplot(223,t,xhat,[],' 14 (c) AutoTrunc[HeaviSine]',[],[])
%
 xhat = AutoTrunc(yDoppler,1.);
 versaplot(224,t,xhat,[],' 14 (d) AutoTrunc[Doppler]',[],[])
 
 
 
  adfig15 -- Adapt Figure 15: AutoSpline: Periodic Thin-Plate spline Reconstructions
% with automatic choice of penalty parameter;
%
%
% Here we depict a reconstruction using penalized smoothing splines,
% in which periodic thin-plate splines are used to smooth the data,
% with roughness penalty proportional to the energy stored in the
% second derivative. The constant of proportionality (or Lagrangian
% tuning parameter) is the empirical minimizer of the Stein Unbiased
% Estimate of Risk.
%
global yblocks ybumps yheavi yDoppler
global t
%
SplineInit
%clf;
[xhat, c] = SplineUChoose(yblocks,2.);
versaplot(221,t,xhat,[],' 15 (a) AutoSpline[Blocks]',[],[])
[xhat,c] = SplineUChoose(ybumps,2.);
versaplot(222,t,xhat,[],' 15 (b) AutoSpline[Bumps]',[],[])
[xhat,c] = SplineUChoose(yheavi,2.);
versaplot(223,t,xhat,[],' 15 (c) AutoSpline[HeaviSine]',[],[])
[xhat,c] = SplineUChoose(yDoppler,2.);
versaplot(224,t,xhat,[],' 15 (d) AutoSpline[Doppler]',[],[])
 
 
 
  s = AutoTrunc(y,sig)
% AutoTrunc -- Perform Truncation
% Usage
% s = AutoTrunc(y,sig)
%
v = fft(y) ./ sqrt(length(y));
n = length(y);
hv = v(1:(n/2));
sv = steintrunc(hv,sig);
vv = [sv 0 conj(sv((n/2):-1:2))];
s = sqrt(length(y)) .* real(ifft(vv));
 
 
  [xhat,yh] = CoronaJS(y,L)
% CoronaJS -- James-Stein applied to ``Littlewood-Paley'' Fourier Coronae
% Usage
% [xhat,xhathat] = CoronaJS(y,L)
% Inputs
% y array of dyadic length 2^J
% NORMALIZED TO NOISE LEVEL 1 (See NOISENORM)
% L Low-Frequency cutoff for shrinkage (e.g. L=4)
% SHOULD BE L << J
% Outputs
% xhat estimate, obtained by James-Stein Shrinkage
% on Dyadic Coronae
% xhhat Fourier Transform of estimate
%
 [n,J] = dyadlength(y);
 yh = fft(y); i=sqrt(-1);
 for j=L:(J-1) ,
 lo = 2^(j-1)+1; hi = 2^j;
 sel = lo:hi;
 v = [real(yh(sel)) imag(yh(sel))];
 vsh = JamesStein(v./sqrt(n/2)).*sqrt(n/2);
 k = 2^(j-1);
 yh(sel) = vsh(1:k) + i*vsh((k+1):(2*k));
 end
 yh((n/2)+1) = 0;
 yh(((n/2)+2):n) = conj(yh((n/2):-1:2));
 xhat = real(ifft(yh));
 
  EvalFix -- Performance of minimax shrinkage for varying sparsity
% parameters C and eps
if ~exist('NProb'),
 NProb = input('Size of Problem');
end
if ~exist('M'),
 M = input('Number of Repetitions');
end
disp('in EvalFix ...')
disp(sprintf('doing %g Monte Carlo repetitions at N=%g',M,N))
disp('this could take a while ...')
FixRecords = zeros( M*13*3 , 3);
MedFix = zeros( 13*3, 3);
v = zeros(1,M); rcount = 1; scount = 1;
for C = [3 5 7],
 for eps = [.005 .01 .02 .04 .06 .08 .10 .12 .14 .16 .18 .20 .25],
 x = zeros(1,NProb);
 k = floor(eps*NProb);
 x(1:k) = zeros(1,k)+C;
 for m=1:M,
 y = x + WhiteNoise(x);
 xhat = MinMaxThresh(y);
 loss = rms(xhat-x); % root mean square!
 v(m) = loss;
 FixRecords(rcount,:) = [C eps loss]; rcount = rcount +1;
 if rem(rcount,10)==0, fprintf('.'), end 
 end
 MedFix(scount,:) = [C eps rms(v)]; scount = scount +1;
 end
end
plot(FixRecords(:,2),FixRecords(:,3).^2,'.')
title(' (b) Fixed Threshold')
xlabel('Epsilon')
ylabel('MSE')
 
 
  EvalHybrid -- Performance of Hybrid(*) shrinkage for varying sparsity 
% parameters C and eps
if ~exist('NProb'),
 NProb = input('Size of Problem');
end
if ~exist('M'),
 M = input('Number of Repetitions');
end
disp('in EvalHybrid')
disp(sprintf('doing %g Monte Carlo repetitions at N=%g',M,N))
disp('this could take a while ...')
HybridRecords = zeros( M*13*3 , 3);
MedHybr = zeros( 13*3, 3);
v = zeros(1,M); rcount = 1; scount = 1;
for C = [3 5 7],
 for eps = [.005 .01 .02 .04 .06 .08 .10 .12 .14 .16 .18 .20 .25],
 x = zeros(1,NProb);
 k = floor(eps*NProb);
 x(1:k) = zeros(1,k)+C;
 for m=1:M,
 y = x + WhiteNoise(x);
 xhat = HybridThresh(y);
 loss = rms(xhat-x); % root mean square!
 v(m) = loss;
 HybridRecords(rcount,:) = [C eps loss]; rcount = rcount +1;
 if rem(rcount,10)==0, fprintf('.'), end 
 end
 MedHybr(scount,:) = [C eps rms(v)]; scount = scount +1;
 end
end
plot(HybridRecords(:,2),HybridRecords(:,3).^2,'*')
hold
title('(c) Hybrid Thresholds')
xlabel('Epsilon')
ylabel('MSE')
hold off
%prtsc
 
 
 
 
  EvalRigor -- Performance of Hybrid(+) shrinkage for varying sparsity 
% parameters C and eps
if ~exist('NProb'),
 NProb = input('Size of Problem');
end
if ~exist('M'),
 M = input('Number of Repetitions');
end
disp('in EvalRigor ...')
disp(sprintf('doing %g Monte Carlo repetitions at N=%g',M,N))
disp('this could take a while ...')
PlusRecords = zeros( M*13*3 , 3);
MedPlus = zeros( 13*3, 3);
v = zeros(1,M); rcount = 1; scount = 1;
for C = [3 5 7],
 for eps = [.005 .01 .02 .04 .06 .08 .10 .12 .14 .16 .18 .20 .25],
 x = zeros(1,NProb);
 k = floor(eps*NProb);
 x(1:k) = zeros(1,k)+C;
 for m=1:M,
 y = x + WhiteNoise(x);
 xhat = RigorThresh(y);
 loss = rms(xhat-x); % root mean square!
 v(m) = loss;
 PlusRecords(rcount,:) = [C eps loss]; rcount = rcount +1;
 if rem(rcount,10)==0, fprintf('.'), end 
 end
 MedPlus(scount,:) = [C eps rms(v)]; scount = scount +1;
 end
end
plot(PlusRecords(:,2),PlusRecords(:,3).^2,'*')
hold
title('(c) Hybrid(+) Thresholds')
xlabel('Epsilon')
ylabel('MSE')
hold off
%prtsc
 
 
 
  EvalSure -- Performance of SURE shrinkage for varying sparsity 
% parameters C and eps
if ~exist('NProb'),
 NProb = input('Size of Problem');
end
if ~exist('M'),
 M = input('Number of Repetitions');
end
disp('in EvalSURE...')
disp(sprintf('doing %g Monte Carlo repetitions at N=%g',M,N))
disp('this could take a while ...')
SureRecords = zeros(M*13*3,3);
MedSure = zeros(13*3,3);
v = zeros(1,M); rcount=1; scount=1;
for C = [3 5 7],
 for eps = [.005 .01 .02 .04 .06 .08 .10 .12 .14 .16 .18 .20 .25],
 x = zeros(1,NProb);
 k = floor(eps*NProb);
 x(1:k) = zeros(1,k)+C;
 for m=1:M,
 y = x + WhiteNoise(x);
 xhat = SUREThresh(y);
 loss = rms(xhat-x); % root mean square!
 v(m) = loss;
 SureRecords(rcount,:) = [C eps loss]; rcount = rcount+1;
 if rem(rcount,10)==0, fprintf('.'), end 
 end
 MedSure(scount,:) = [C eps rms(v)]; scount = scount+1;
 end
end
plot(SureRecords(:,2),SureRecords(:,3).^2,'+')
title('(a) SURE Threshold')
xlabel('Epsilon')
ylabel('MSE')
 
 
 
 
  IMJPlotMultiRes(wavecoef,L,scal,QMF_Filter)
% IMJPlotMultiRes -- Mallat-Style MRA plot -- IJ version with changed axis labelling
%
if scal==0, 
 scal = 1. / max(abs(wavecoef)); 
end
[n,J] = dyadlength(wavecoef); 
LockAxes([0 1 (L-2) J]);
t = (.5:(n-.5))/n;
%
 for j=(J-1):-1:L
 w = zeros(1,n);
 w(dyad(j)) = wavecoef(dyad(j));
 x = IWT_PO(w,j,QMF_Filter);
 plot(t, j + scal.*x);
 end
 w = zeros(1,n);
 w(1:2^(L)) = wavecoef(1:2^(L));
 x = IWT_PO(w,L,QMF_Filter);
 plot(t,(L-1) + scal.*x);
%xlabel('Spatial Var.'); ylabel('Dyadic Resolution Level');
%title('Multi-Resolution Decomposition');
UnlockAxes;


% Revision History
% 10/1/05 AM Name of the variable QMF is changed to
% QMF_Filter

 
  xhat = JamesStein(y)
% JamesStein -- James-Stein Shrinker
% Usage:
% xhat = JamesStein(y)
% Inputs:
% y Noisy Data, Standard Deviation of Noise = 1.
% Outputs:
% xhat James-Stein Estimate of x
% Notes:
% dimension(x) > 2 is required.
%
 ss = sum(y.^2);
 n = length(y);
 shrnk = max((ss - (n-2))/ss,0);
 xhat = shrnk .* y;

 
 
  urisk = opttrunc(y)
% opttrunc -- Perform Truncation
% Usage
% urisk = opttrunc(y)
%

tails = reverse(cumsum(reverse(abs(y).^2)));
n = length(y);
penalty = (-n):2:(n-2);
urisk = tails + penalty;
 
 
  xhat = RigorThresh(y)
% RigorThresh -- Adaptive Threshold Selection using SURE
% Usage
% xhat = RigorThresh(y)
% Inputs
% y Noisy Data with Std. Deviation = 1
% Outputs
% xhat Estimate of mean vector
% Notes
% Using Principle of Stein's Unbiased Risk Estimate
% Variant with Rigorous Proof
%
% See: Donoho, D.L. and Johnstone, I.M. ``Adapting
% to Unknown Smoothness by Wavelet Shrinkage''

 [n,J] = dyadlength(y);
 magic = sqrt(2*log(n));

 % pretest statistic
 eta = (norm(y).^2 - n)/n;
 crit = J^(1.5)/sqrt(n);

 if eta < crit,

 % failed pretest -- use sqrt(2 log(n))
 xhat = SoftThresh(y,magic);

 else
 
 % Construct random subsets of size n/2 each
 xhat = rand(size(y)); v = sort(xhat);
 I = xhat > v(n/2); Ip = ~I;

 % Apply SURE to each half
 T = min(ValSUREThresh(y(I)), magic);
 Tp = min(ValSUREThresh(y(Ip)),magic);
 xhat = y;
 
 % Threshold from one half is applied to other
 xhat(I) = SoftThresh(y(I) ,Tp);
 xhat(Ip) = SoftThresh(y(Ip),T);

 end

 
 
  r = rms(y)
% rms -- Return Root Mean Square of Signal
% Usage
% r = rms(y)
% Inputs
% y signal
% Outputs
% r sqrt( mean ( y .^2 ))
%
r = sqrt( mean ( y .^2 ));
 
 

  SplineInit -- Script to initialize Spline search 
% Must run before SplineUChoose, SplineURisk
global Spline_Pnty
global Spline_Data
 
 

  [xh,c] = SplineUChoose(y,k)
% SplineUChoose -- Estimate Optimal Bandwidth for Spline Kernel
% Usage
% [xhat,risk] = SplineUChoose(y,k)
% Inputs
% y Noisy Data Normalized to Noise Level 1
% k Degree of Spline Penalty .
% Outputs
% xhat Estimate based on y with optimal lambda
% risk Estimate of Risk at this x
% Warning
% must run script SplineInit before running this
% See Also:
% SplineInit, SplineURisk
%

 global Spline_Pnty
 global Spline_Data
%
 f = fft(y);
 n = length(y);
 Spline_Data = abs(f/sqrt(n)).^2;
 Spline_Pnty = (2 .* sin(pi .* (0:(n-1)) ./ n)) .^ (2*k);
 Spline_Lmda = fminbnd('SplineURisk',0.,n.^2);
 xhh = f ./ (1 + Spline_Lmda.*Spline_Pnty);
 xh = real(ifft(xhh));
 %clf; subplot(211);
 %plot(1. ./ (1 + Spline_Lmda.*Spline_Pnty));
 %plot(xh);
 c = SplineURisk(Spline_Lmda);

 
 

  URisk = SplineURisk(lam)
% SplineURisk -- Unbiased Estimate of Risk of a given Spline Filter
% Usage
% URisk = SplineURisk(lam)
% Inputs
% lam Real Number in range [0,n]
% Global Variables
% Spline_Pnty positive variable of length n.
% Spline_Data positive variable of length n.
% [These Variables are initialized by SplineUChoose]
% Output
% Estimated Risk of Spline Est implicitly defined by
% 1/(1 + lam * 2^2k sin^2k(pi*j/n))
% Warning
% SplineInit and SplineUKern must be run first!
%
global Spline_Pnty
global Spline_Data
%
compl = (lam .* Spline_Pnty)./(1. + lam .* Spline_Pnty);
URisk = sum(compl.^2 .*Spline_Data) - 2 .* sum(compl) ...
 + length(Spline_Data);
 
 
  x = steintrunc(y,sig)
% steintrunc -- Perform Truncation
% Usage
% x = steintrunc(y,sig)
%

u = opttrunc(y ./sig);
t = min(u);
indmin = find(u==t);
x = y .* ((1:length(y)) <= indmin);
 
 
  [xh,wcoef] = WaveJS(y,L,qmf)
% WaveJS -- James-Stein Shrinkage Applied to Wavelet dyads
% Usage
% [xh,xwh] = WaveJS(y,L,qmf)
% Inputs
% y array of dyadic length 2^J
% NORMALIZED TO NOISE LEVEL 1 (See NOISENORM)
% L Low-Frequency cutoff for shrinkage (e.g. L=4)
% SHOULD BE L << J
% qmf Quadrature Mirror Filters for Wavelet Transform
% Outputs
% xh estimate, obtained by James-Stein Shrinkage
% on Wavelet Coefficients
% xwh Wavelet Transform of estimate
% See Also
% FWT_PO, IWT_PO, JamesStein
%

 [n,J] = dyadlength(y); 

 wcoef = FWT_PO(y,L,qmf) ;

 for j=(J-1):-1:L
 wcoef(dyad(j)) = JamesStein(wcoef(dyad(j))) ;
 end

 xh = IWT_PO(wcoef,L,qmf);
 

 

 
 
 
 
  ******* ./Wavelab850/Papers/Asymp/ *******
  asfig01 -- Asymptopia Figure 01: De-Noising of NMR Signal
  asfig02 -- Asymptopia Figure 02: Four Spatially Inhomogenoeus Signals
  asfig03 -- Asymptopia Figure 03: Noisy Versions of Four Signals
  asfig04 -- Asymptopia Figure 04: Wavelet Shrinkage Damping of Noisy Data
  asfig05 -- Asymptopia Figure 05: AutoSpline Reconstructions from Noisy Data
  asfig06 -- Asymptopia Figure 06: AutoTrunc Truncated Fourier Reconstructions
  asfig07 -- Asymptopia Figure 07: Unconditional Bases and Mutilation of Coefficients
  asfig08 -- Asymptopia Figure 08: Smoothing Good-Gaskins data
  asfig09 -- Asymptopia Figure 09: Comparison to Good-Gaskins Smoother
  asfig10 -- Asymptopia Figure 10: Compare various thresholds; ESCA data
  disp('File ESCAsim.mat ~exist -- calculating')
  disp('File ESCAsim.mat exists -- loading')
  AsympCleanup -- remove all global variables
  AsympFig -- callback implementing functions of AsympDemo
  AsympInit -- setup all global variables
  AsympIntro -- Info for AsympDemo
  fishrand -- generate Poisson random vector, with given intensity
  reverse -- Reverse order of elements in 1-d signal
 
 
  ******* ./Wavelab850/Papers/Blocky/ *******
  aifig1 -- Average-Interpolating Refinement of a sequence of averages.
%
% The figure shows (as a step function) the block averages 
% of the function sin(t), and (as dashed and dotted curves)
% several average-interpolating refinements to a resolution
% 5 scales finer.
%
% At the resolution of the plotting device, all of the refinements 
% are equal to the true underlying function, and are much closer to 
% that function than the piecewise constant function. 
%
global Fil2 Fil4 Fil6 Fil8
global EdgeFil2 EdgeFil4 EdgeFil6 EdgeFil8
%
t = (0:63) ./16 - 2;
ave = 16 .*(cos(t) - cos(t +1/16.));
Blocks = ave;
MaxOrd=8;
dt = t(2)-t(1);
vb = [ Blocks ; Blocks ]; vb = vb(:);
tb = [t ; (t+dt) ]; tb = tb(:);
%clf;
hold on
subplot(1,1,1)
plot(tb,vb)
title('1. Interpolating Averages of sin(t)')
%axis([0 1 0 1]); 
xlabel('t')
%
a = Blocks;
for j=1:5, 
 a = AIRefine(a,2,Fil2,EdgeFil2); 
end;
t = 4 * (0:(length(a)-1)) ./ length(a) - 2;
plot(t,a,'--')
%
a = Blocks;
for j=1:5, 
 a = AIRefine(a,4,Fil4,EdgeFil4); 
end;
plot(t,a,'-.')
%
if(MaxOrd > 6),
a = Blocks;
for j=1:5, 
 a = AIRefine(a,8,Fil8,EdgeFil8); 
end;
plot(t,a,':')
end;
UnlockAxes






 
 
  aifig2 -- Fundamental Solutions of Average-Interpolation
%
% Results of average-interpolating the Kronecker
% sequence. (Panel a) Average-Interpolation of Degree D=2.
% (Panel b) D=4. (Panel c) D=6. (Panel d) D=8.
%
%
global Fil2 Fil4 Fil6 Fil8
global EdgeFil2 EdgeFil4 EdgeFil6 EdgeFil8
%
nref = 5;
kronecker = zeros(1,16);
kronecker(9) = 1;
%
% D=2
%
 a = kronecker;
 for j=1:nref, 
 a = AIRefine(a,2,Fil2,EdgeFil2); 
 end;
 t = 16*(0:(length(a)-1)) ./ length(a) -8; 
 versaplot(221,t,a,[],'2(a) AI Fundamental Soln., D=2',[],[])
 a2 =a;
%
% D=4
%
 a = kronecker;
 for j=1:nref, 
 a = AIRefine(a,4,Fil4,EdgeFil4); 
 end;
 versaplot(222,t,a,[],'2(b) AI Fundamental Soln., D=4',[],[])
 a4=a;
%
% D=6
%
 a = kronecker;
 for j=1:nref, 
 a = AIRefine(a,6,Fil6,EdgeFil6); 
 end;
 versaplot(223,t,a,[],'2(c) AI Fundamental Soln., D=6',[],[])
 a6=a;
%
% D=8
%
 kronecker = zeros(1,32);
 kronecker(17) = 1;
 a = kronecker;
 for j=1:nref, 
 a = AIRefine(a,8,Fil8,EdgeFil8); 
 end;
 t = (0 : (length(a)-1)) .* 32 ./ (length(a)) - 16;
 inx = (t >= -8) & (t < 8); a8=a;
 versaplot(224,t(inx),a(inx),[],'2(d) AI Fundamental Soln., D=8',[],[])
 a8=a;
%
UnlockAxes;
 
 
  aifig3 -- AI Wavelets
%
% The result of average-interpolating a sequence which
% consists of (-1,1) in a sea of zeros.
%
% 
global Fil2 Fil4 Fil6 Fil8
global EdgeFil2 EdgeFil4 EdgeFil6 EdgeFil8
%
nref = 4;
kronecker = zeros(1,32);
kronecker(17) = -1; kronecker(18) = 1;
%
 a = kronecker;
 for j=1:nref, 
 a = AIRefine(a,2,Fil2,EdgeFil2); 
 end;
 t = 16*(0:(length(a)-1)) ./ length(a) -8; 
 versaplot(221,t,a,[],'3(a) AI Wavelet, D=2',[],[]) ; 
 a2 =a;
%
% D=4
%
 a = kronecker;
 for j=1:nref, 
 a = AIRefine(a,4,Fil4,EdgeFil4); 
 end;
 versaplot(222,t,a,[],'3(b) AI Wavelet, D=4',[],[]); 
 a4=a;
%
% D=6
%
 a = kronecker;
 for j=1:nref, 
 a = AIRefine(a,6,Fil6,EdgeFil6); 
 end;
 versaplot(223,t,a,[],'3(c) AI Wavelet, D=6',[],[]) ; 
 a6=a;
%
% D=8
%
 kronecker = zeros(1,64);
 kronecker(33) = -1; kronecker(34) = 1;
 a = kronecker;
 for j=1:nref, 
 a = AIRefine(a,8,Fil8,EdgeFil8); 
 end;
 t = (0 : (length(a)-1)) .* 32 ./ (length(a)) - 16;
 inx = (t >= -8) & (t < 8); a8=a;
 versaplot(224,t(inx),a(inx),[],'3(d) AI Wavelet, D=8',[],[])
 
 
  aifig4 -- The Dual Functions
%
% These functions are dual to the wavelets of the previous
% figure, in that integration of a function f against these
% functions gives the coefficients in the wavelet expansion
% using the corresponding AI wavelets.
%
global Fil2 Fil4 Fil6 Fil8
global EdgeFil2 EdgeFil4 EdgeFil6 EdgeFil8
%
 nref = 5;
v = zeros(4,32);
kronecker = zeros(1,32);
kronecker(18) = 1;
% Order D=2
D = 2;
a = 2.* kronecker;
for k=(-D/2):(D/2),
 a(17+2*k) = a(17+2*k) - Fil2(1 + D + 2*k);
 a(18+2*k) = a(18+2*k) - Fil2(1 + D + 2*k);
end
% D, a 
v(1,:) = a; 
for j=1:nref, 
 a = [a ; a]; 
 a = reshape(a,1,prod(size(a)));
end;
t = 16*(0:(length(a)-1)) ./ length(a) -8; 
inx = (t >= -5) & (t < 5);
subplot(221);
plot(t(inx),a(inx)) ; 
title('4(a) Dual AI Wavelet D=2')

D = 4;
a = 2 .* kronecker;
for k=(-D/2):(D/2),
 a(17+2*k) = a(17+2*k) - Fil4(1 + D + 2*k);
 a(18+2*k) = a(18+2*k) - Fil4(1 + D + 2*k);
end
% D, a 
v(2,:) = a; 
for j=1:nref, 
 a = [a ; a]; 
 a = reshape(a,1,prod(size(a)));
end;
subplot(222);
plot(t(inx),a(inx)) ; 
title('4(b) Dual, D=4')

D = 6;
a = 2 .* kronecker;
for k=(-D/2):(D/2),
 a(17+2*k) = a(17+2*k) - Fil6(1 + D + 2*k);
 a(18+2*k) = a(18+2*k) - Fil6(1 + D + 2*k);
end
% D, a
v(3,:) = a; 
for j=1:nref, 
 a = [a ; a]; 
 a = reshape(a,1,prod(size(a)));
end;
subplot(223);
plot(t(inx),a(inx)) ; 
title('4(c) Dual, D=6')

D = 8;
a = 2.* kronecker;
for k=(-D/2):(D/2),
 a(17+2*k) = a(17+2*k) - Fil8(1 + D + 2*k);
 a(18+2*k) = a(18+2*k) - Fil8(1 + D + 2*k);
end
%D, a
v(4,:) = a; 
for j=1:nref, 
 a = [a ; a]; 
 a = reshape(a,1,prod(size(a)));
end;
t = 16*(0:(length(a)-1)) ./ length(a) -8;
subplot(224);
plot(t(inx),a(inx)) ; 
title('4(d) Dual, D=8')
 
 

 
  aifig5 -- Show Boundary Wavelets
%
% This figure displays some boundary wavelets and
% corresponding wavelets at the interior of the 
% interval. The wavelets at the ``heart'' of the 
% interval are the same as the non-boundary adjusted 
% wavelets. The wavelets at the ``ends'' are modified
% to maintain polynomial exactness and related 
% properties.
%
global Fil2 Fil4 Fil6 Fil8
global EdgeFil2 EdgeFil4 EdgeFil6 EdgeFil8
%
% D=2 Boundary Wavelets
%
 
 subplot(221); 
 LockAxes([0 1 -20 20]);
 for j= [1 2 5 7 8],
 a = zeros(1,8); 
 a(j) = 1; 
 for k=1:5, 
 a = AIRefine(a,2,Fil2,EdgeFil2); 
 end;
 t = (0:255) ./256; 
 if j==1, pch = '-.'; end
 if j==2, pch = '--'; end
 if j==5, pch = '-'; end
 if j==7, pch ='--'; end
 if j==8, pch = '-.'; end
 plot(t,8 .*a, pch); 
 end
 title('5(a) Interval Wavelets phi_{3,k}, D=2'); 
 UnlockAxes;
%
% D=4
%
subplot(222); 
LockAxes([0 1 -20 20]);
for j= [1 2 9 15 16],
a = zeros(1,16); 
a(j) = 1; 
for k=1:5, 
 a = AIRefine(a,4,Fil4,EdgeFil4); 
end;
t4 = (0:(length(a)-1)) ./length(a); 
if j==1, pch = '-.'; end
if j==2, pch = '--'; end
if j==9, pch = '-'; end
if j==15, pch ='--'; end
if j==16, pch = '-.'; end
plot(t4,8 .*a, pch); 
end
title('5 (b) Interval Wavelets phi_{4,k}, D=4');
UnlockAxes;
%
% D=6
%
subplot(223); 
LockAxes([0 1 -20 20]); 
for j= [1 2 5 7 8],
a = zeros(1,16); 
a(2*j) = 1; a(2*j-1) = -1; 
for k=1:4, 
 a = AIRefine(a,2,Fil2,EdgeFil2); 
end;
t = (0:255) ./256; 
if j==1, pch = '-.'; end
if j==2, pch = '--'; end
if j==5, pch = '-'; end
if j==7, pch ='--'; end
if j==8, pch = '-.'; end
plot(t,8 .*a, pch); 
end
title('5 (c) Interval Wavelets psi_{3,k}, D=2');
UnlockAxes;
%
% D=8
%
subplot(224); 
LockAxes([0 1 -20 20]); 
for j= [1 2 9 15 16],
a = zeros(1,32); 
a(2*j) = 1; a(2*j-1) = -1; 
for k=1:4, 
 a = AIRefine(a,4,Fil4,EdgeFil4); 
end;
if j==1, pch = '-.'; end
if j==2, pch = '--'; end
if j==9, pch = '-'; end
if j==15, pch ='--'; end
if j==16, pch = '-.'; end
plot(t4,8 .*a, pch); 
end
title('5 (d) Interval Wavelets psi_{4,k}, D=4');
UnlockAxes
 
 
 
  aifig6 -- Depict 2-d Average-Interpolation
%
% (Panel a) A set of 2-d block averages at a fine scale,
% (Panel b) a set of 2-d block averages two scales coarser,
% (Panel c) Refinement of the coarse-scale averages (D=2)
% (Panel c) Refinement of the coarse-scale averages (D=4)
%
% Refinement of the coarse-scale information nearly reproduces
% the original fine-scale data.
%
global Fil2 Fil4 Fil6 Fil8
global EdgeFil2 EdgeFil4 EdgeFil6 EdgeFil8
%
n = 64; lev=2;
[x,y] = ...
meshgrid((-.501+1/(2*n)):(1/n):(.501-1/(2*n)),(-.5+1/(2*n)):(1/n):(.5-1/(2*n)));
y = flipud(y);
r = sqrt(x .^2 + y .^2);
subplot(221)
sr = (sin(5 .*(r)) ./ (5 .*(r))).^2; clear r x y;
mesh(sr); title('6(a) Averages of (sin(5r)/(5r))^2, Fine scale, j=6')
csr = HCoarsen2d(sr,lev);
clear sr;
Hsr = HRefine2d(csr,lev);
subplot(222);
mesh(Hsr); title('6(b) Averages at Coarse scale, j=4')
clear Hsr; 
Rsr2 = AIRefine2d(csr,lev,2,Fil2,EdgeFil2);
subplot(223);
mesh(Rsr2); title('6(c) D=2 Interpolation from Coarse to Fine') 
clear Rsr2;
Rsr4 = AIRefine2d(csr,lev,4,Fil4,EdgeFil4);
subplot(224);
mesh(Rsr4); title('6(d) D=4 Interpolation from Coarse to Fine');

% Revision History
% 09/29/99 MRD meshdom replaced by meshgrid/flipud
 
 
  aifig7 -- Noisy Differentiation
%
% Here we illustrate the use of wavelets for noisy differentiation.
% (Panel a) displays object Bumps, (panel b) its noisy cumulative, 
% and (Panel c) a naive inversion based on differencing. 
%
% The naive inversion is very noisy and displays only a hint of the 
% structure of the underlying function.
%
% We also display a wavelet-shrinkage reconstruction, in which noise is
% suppressed while the structure remains. 
%
% The wavelet-shrinkage reconstruction goes as follows
%
% 1. Transform the naive reconstruction into the wavelet domain.
% 2. Apply level-dependent thresholds which scale inversely with 
% resolution.
% 3. Return to the wavelet domain.
%
% Coiflets with 3 vanishing moments are used.
%

 n = 2048;
 t = (.5:(n-.5)) ./ n;
 Bumps = MakeSignal('Bumps',n);
%
% Build Noisy Primitive
%
 zBumps = cumsum(Bumps);
 x = zBumps + rand(size(zBumps));
%
% Naive Reconstruction
%
 z = diff([0 x]);
%
% Wavelet Shrinkage
%
 QSymm6 = MakeONFilter('Symmlet',6);
 w = FWT_PO(z, 5, QSymm6);
 ws = InvShrink(w,5,4,1);
 zrec = IWT_PO(ws,5,QSymm6);
%
% Display Results

 subplot(221); plot(t,Bumps); title('7 (a) Object Bumps')
 subplot(222); plot(t,x); title('7 (b) Noisy Primitive')
 subplot(223); plot(t,z); title('7 (c) Naive Reconstruction via Differences')
 subplot(224); plot(t,zrec); title('7 (d) Wavelet De-Noising Reconstruction')
 
 
 
  BlockyCleanup -- Remove global datastructures for BlockyDemo
clear global Fil2 Fil4 Fil6 Fil8
clear global EdgeFil2 EdgeFil4 EdgeFil6 EdgeFil8
clear labels callbacks header
 
 
  BlockyDemo -- Show Blocky Figs with Point & Click Interface
% To run, type BlockyDemo
header = 'Blocky Figures';
labels = ['Change Settings' ; 'Run All Scripts' ; ...
 'Run One Script ' ; 'Clean Up '];
callbacks = ['BlockyFig(19)'; 'BlockyFig(20)'; ...
 'BlockyFig(21)'; 'BlockyCleanup'];
global BLOCKYFIGNUM 
BLOCKYFIGNUM = 0;
clc; help('BlockyIntro');
%choices('Blocky',header,labels,callbacks); 

NumberOfChoices=4;
choice=0;
while (choice<NumberOfChoices)
 choice=menu(header,'Change Settings', 'Run All Scripts','Run One Script', 'Clean Up and Exit');
 eval(callbacks(choice,:));
end

clear header labels callbacks


 
  BlockyFig(num)

% BlockyFig -- callback implementing functions of BlockyDemo
global BLOCKYFIGNUM
BlockyInit
if 1 <= num & num <=7 ,
 nihil = ~exist('BLOCKYFIGNUM');
 if ~nihil,
 empty = isempty(BLOCKYFIGNUM);
 else
 empty = 1;
 end
 if ~empty,
 nil = BLOCKYFIGNUM == 0;
 else
 nil = 1;
 end
 if nil,
 BlockyInit
 end
 figname = sprintf('aifig%g',num);
 clc; home; help(figname);
 eval(figname);
 RegisterPlot(figname,'{R}');
 drawnow;
 BLOCKYFIGNUM = num;
elseif num==19, % change default settings
 clc; home
 disp('There Are No Defaults to CHANGE');
elseif num==20, %run all figures
 BlockyInit
 for num=1:7,
 figname = sprintf('aifig%g',num);
 eval(figname);
 RegisterPlot(figname,'{R}');
 clc; home; help(figname);
 drawnow;
 BLOCKYFIGNUM = num;
 end
else
 disp('Bizarre argument in BlockyFig')
 disp(sprintf('num =%g',num))
end




  BlockyInit -- setup global data structures for BlockyDemo
global Fil2 Fil4 Fil6 Fil8
global EdgeFil2 EdgeFil4 EdgeFil6 EdgeFil8
%
 MakeAllFilters
 
 
 
  BlockyIntro -- Info about BlockyDemo
%
% The files in this directory may be used to reproduce the
% figures in the paper
%
% ``Smooth Wavelet Decompositions with Blocky Coefficient Kernels''
% by D.L. Donoho
%
% The paper may be obtained by WWW from 
% http://www-stat.stanford.edu/~donoho/Reports/1993/blocky.ps.Z.
%
% The paper shows how to make wavelets which satisfy prescribed
% boxcar averages and illustrates some applications of these.
%

 
 

 
  MakeAllFilters -- Create AI Filters and Edge Filters
global Fil2 Fil4 Fil6 Fil8
Fil2 = MakeAIFilter(2); Fil4 = MakeAIFilter(4);
Fil6 = MakeAIFilter(6); Fil8 = MakeAIFilter(8);
% MakeBdryFilters
global EdgeFil2 EdgeFil4 EdgeFil6 EdgeFil8
EdgeFil2 = MakeAIBdryFilter(2); 
EdgeFil4 = MakeAIBdryFilter(4);
EdgeFil6 = MakeAIBdryFilter(6); 
EdgeFil8 = MakeAIBdryFilter(8);
 
 
 
 
  ******* ./Wavelab850/Papers/Correl/ *******
  corfig05 -- Correl Figure 05: Denoised estimates of signal and 
  CorrelCleanup -- Clean up Correl global variables
  CorrelFig -- callback for display of individual figures in CorrelDemo
  CorrelInit -- Initialize all global variables for CorrelDemo
  CorrelIntro -- Info about CorrelDemo
  IMJPlotSpikes -- Plot 1-d signal as baseline with series of spikes
  LevelShrinkS -- Wavelet Shrinkage with universal threshold and 
  MultiSUREThresh -- 
  IMJPlotWaveCoeff -- Wavelet Coeffts Plot: IMJ version with 
  ScaleMaker -- Calculate level dependent estimates of scale
  StatShrink -- Wavelet shrinkage using Stationary wavelet transform
  StepPlot -- split signal into nlayers equal pieces for plotting
 
 
  ******* ./Wavelab850/Papers/Ideal/ *******
  claxis -- Clear plot axis
  IdealCleanup -- Clean up Ideal global variables
  IdealFig -- callback implementing functions of IdealDemo
  IdealInit -- setup all global variables
  IdealIntro -- Info about IdealDemo
  IdealKern -- Estimate using Ideal Kernel
  IdealWave -- Estimate using Ideal Wavelet Damping
  IdealWts -- Ideal weights for Diagonal Linear Estimates
  idfig01 -- Ideal Figure 01: Four Spatially Inhomogenoeus Signals
  idfig02 -- Ideal Figure 02: S8 Wavelet Coefficients of the Four Signals
  idfig03 -- Ideal Figure 03: Noisy Versions of Four Signals
  idfig04 -- Ideal Figure 04: S8 Wavelet Coefficients of Four Noisy datasets
  idfig05 -- Ideal Figure 05: Ideal Wavelet Damping of Noisy Data
  idfig06 -- Ideal Figure 06: Ideal Damping in Wavelet Domain
  idfig07 -- Ideal Figure 07: Wavelet Shrinkage Damping of Noisy Data
  idfig08 -- Ideal Figure 08: Wavelet Shrinkage in Wavelet Domain
  idfig09 -- Ideal Figure 09: Examples of VisuShrink Thresholding
  idfig10 -- Ideal Figure 10: Ideal Fourier Damping of Noisy Data
  IMJPlotWaveCoeff -- Wavelet Coeffts Plot: IMJ version with changed axis labelling
 
 
  ******* ./Wavelab850/Papers/MinEntSeg/ *******
  BloxPlot -- Make blox plot
  Compare2 -- Compare two compression curves, semilog plot
  FastAllSeg -- Fast All-Segmentations algorithm
  FastEntProfile -- Fast Entropy Profile
  MakeAIRightFilter -- Make Average Interpolating Right Filter
  MakeAISegFilter -- Extract coefficients for segmented B/W Transform
  MakeBumps -- Make Bumps Object for Segmentation Pursuit
  MakeMESTestSig -- Make MES Test Signal
  mefig201 -- Visualize Average-Interpolation
  mefig202 -- Portray Average-Interpolating (Boundary) Wavelets
  mefig203 -- Visualize Segmented Average Interpolation
  mefig204 -- Compare Ordinary AI & Segmented AI Refinements
  mefig205 -- Four Test Functions
  mefig206 -- Wavelet Transforms w/o segmentation
  mefig207 -- MRA plots without segmentation
  mefig208 -- Segmented wavelet transform for the Four Test Signals
  mefig209 -- MRA plots w/segmentation
  mefig210 -- Compression numbers
  mefig211 -- De-Noising "Ramp"; Segmented transform; ideal segmentation.
  mefig212 -- Segmented vs. Ordinary de-noising; object "Ramp"
  mefig213 -- De-Noising "Cusp"; Segmented Transform; ideal segmentation.
  Note -- this script makes calculations which are used later by mefig214.m
  mefig214 -- Segmented vs. Ordinary de-noising; object "Cusp" 
  Note -- this script uses results calculated in mefig213.m
  mefig302 -- Segmented Wavelet Transform at variable locations 
  mefig303 -- Entropy Profiles, Object "Ramp"
  mefig304 -- Segmented Wavelet Transform of "Cusp"; variable segmentations
  mefig305 -- Entropy Profiles, Object "Cusp"
  mefig306 -- Ideal and SURE Risk Profiles, objects "Ramp," "Cusp" 
  mefig401 -- Multi-Resolution Filter Bank output, object "Ramp"
  mefig402 -- Entropy Profiles for Object "Blocks"
  mefig601 -- Entropy Profiles for Object "Bumps"
  mefig602 -- MultiResolution Filter Bank Output, Object "Bumps"
  mefig603 -- Segmentation Components of object "Bumps"
  mefig604 -- Segmentation Pursuit, object "Bumps"
  mefig605 -- Residuals from Segmentation Pursuit, object "Bumps"
  MESCleanup -- Clean up global variables from MES Demo
  MESFig -- Callback implementing functions of MESDemo
  MESInit -- Setup global variables for MES Demo
  MESIntro -- Info for MESDemo
  off_filter -- Offset-index convolution
  PrintMEFigs -- Print All Figures of MinEntSegPaper
  Seg_Pursuit -- Seg Pursuit
  SharpPart -- Singularity Surfer
 
 
  ******* ./Wavelab850/Papers/MIPT/ *******
  mifig21 -- Visualize Median Interpolation, D=2.
  mifig22 -- Nonlinearity of \pi_{(2)}
  mifig23 -- Visualize Median Interpolation, D=4.
  mifig24 -- Nonlinearity of Median Refinement
  mifig31 -- MI vs LMI
  mifig41 -- Noiseless Data Transforms (MIPT coefficients)
  mifig42 -- Noiseless Data Transforms (MIPT coefficients)
  mifig51 -- De-Noising of Gaussian White Noise
  mifig52 -- De-Noising of Cauchy White Noise
  MedianCleanup -- remove all global variables
  MedianFig -- callback implementing functions of MedianDemo
  MedianInit -- setup all global variables
  MIPTIntro -- Info for MIPTDemo
  MIPTThreshold -- Determine MIPT thresholds for denoising
  mitable51 -- Theoretically Derived MIPT Thresholds
 
 
  ******* ./Wavelab850/Papers/RiskAnalysis/ *******
  Copyright -- Sudeshna Adak, Dept. of Statistics,
  Copyright -- Sudeshna Adak, Dept. of Statistics,
  Copyright -- Sudeshna Adak, Dept. of Statistics,
  Copyright -- Sudeshna Adak, Dept. of Statistics,
  Copyright -- Sudeshna Adak, Dept. of Statistics,
  Copyright -- Sudeshna Adak, Dept. of Statistics,
  Copyright -- Sudeshna Adak, Dept. of Statistics, 
  do_button -- Called by WLBrowser
  disp('Name of Data -- file name without extension')
  Copyright -- Sudeshna Adak, Dept. of Statistics,
  NewZoom -- Called by WLBrowser
  Copyright -- Sudeshna Adak, Dept. of Statistics,
  -- Data for Actions menu #3 --- %
  -- Data for Data menu # 1 --- %
  -- Data for Settings menu #3 --- %
  -- Data for Signals menu #2 --- %
  RiskFig -- This implements making of the figures by the RiskDemo
  -- Data for Actions menu #3 --- %
  -- Data for Data menu # 1 --- %
  -- Data for Settings menu #3 --- %
  -- Data for Signals menu #2 --- %
  sonify -- Called by WLBrowser
 
 
  ******* ./Wavelab850/Papers/ShortCourse/ *******
  FHT_Med -- Forward Med Transform
  FigNameString -- Append two-digit number to string
  FWT_SegAI -- Forward Segmented Transform
  HCoarsen -- Coarsen object
  IHT_Med -- Inverse Med Transform
  IWT_SegAI -- Inverse Segmented Transform
  MakeAISegmentFilter -- coefficients for segmented B/W Transform
  MedTriadDown -- Med Triad Down
  MedTriadUp -- Med Triad Up
  PlotSAIMultiRes -- Plot Segmented Multi-Resolution Analysis
  SAIDyadDown -- SAI Dyad Down
  SAIDyadUp -- SAI Dyad Up
  SCCleanup -- Clean up Short Course global variables
  scfig01 -- Short Course 01: De-Noising of NMR Signal
  scfig02 -- Short Course 02: Noisy Deconvolution
  scfig03 -- Short Course 03: Comparison of Wavelet and Packet DeNoising
  scfig04 -- Short Course 04: Four Spatially Inhomogenoeus Signals
  scfig05 -- Short Course 05: Noisy Versions of Four Signals
  scfig06 -- Short Course 06: Wavelet Shrinkage of the Four Noisy Signals
  scfig07 -- Short Course 07: AutoSpline Reconstructions from Noisy Data
  scfig08 -- Short Course 08: AutoTrunc Truncated Fourier Reconstructions
  scfig09 -- Short Course 09: WaveShrink of object yBlocks in Haar Basis
  scfig10 -- Short Course 10: Comparing Compression Abilities
  scfig11 -- Short Course 11: Compare Wavelet and DCT Partial Reconstructions
  scfig12 -- Short Course 12: DeNoising a 2-d object
  scfig13 -- Short Course 13: Smoothing Counts data by square roots; ESCA data
  scfig14 -- Short Course 14: Estimating Time Series Spectrum
  scfig15 -- Short Course 15: Noisy Differentiation by WVD
  scfig16 -- Short Course 16: Noisy Differentiation by WVD in Wavelet Domain
  scfig17 -- Short Course 17: Noisy Differentiation by Ideal Fourier Damping
  scfig18 -- Short Course 18: Depict Deconvolution in Wavelet Domain
  scfig19 -- Short Course 19: Display Vaguelette Kernels
  scfig20 -- Short Course 20: Four Time-Frequency test signals
  scfig21 -- Short Course 21: Four Noisy Time-Frequency test signals
  scfig22 -- Short Course 22: Wavelet Packet DeNoising
  ``Cosine Packet'' -- adaptively selected Cosine Packet basis
  ``Fourier'' -- discrete cosine transform
  scfig23 -- Short Course 23: Compare Four DeNoising methods
  ``Wavelet'' -- nearly symmetric wavelet with 8 vanishing moments
  ``Wavelet Packet'' -- adaptively selected Wavelet Packet basis
  scfig24 -- Short Course 24: Compare Segmented and Ordinary Refinement
  scfig25 -- Short Course 25: DeNoising a segmented transform
  scfig26 -- Short Course 26: Sine signal with Cauchy Contamination
  scfig27 -- Short Course 27: Linear Smoothing of Cauchy Noise
  scfig28 -- Short Course 28: Robust De-Noising
  SCFig -- Callback implementing functions of SCDemo
  SCInit -- Setup global datastructures for ShortCourse
  SCIntro -- Info for SCDemo
  SegRefine -- Perform refinement
  SRefine -- Perform refinment
  triadLength -- Return length of triad
 
 
  ******* ./Wavelab850/Papers/SpinCycle/ *******
  CSInit -- setup global datastructures for SpinCycle
  cspinf01 -- Four Spatially Inhomogeneous Signals
  cspinf02 -- Noisy Versions of Four Signals
  cspinf03 -- Ordinary S8 Wavelet Shrinkage
  cspinf04 -- Cycle Spinning S8 Wavelet Shrinkage
  cspinf05 -- Fully Translation-Invariant S8 Wavelet Shrinkage
  cspinf06 -- Ordinary Haar Wavelet Shrinkage
  cspinf07 -- Fully TI Haar Wavelet Shrinkage
  cspinf08 -- Closeups of Fully TI Haar Wavelet Shrinkage
  cspinf09 -- Fully TI, S8, Hard-Thresholding Shrinkage
  cspinf10 -- TI, Hard, Haar Shrinkage, Small Sample Size
  cspinf11 -- Cycle Spinning SURE Threshold Wavelet Shrinkage
  cspinf12 -- Two examples for time-frequency de-noising
  cspinf13 -- Ordinary Non-Spinning CP De-Noising
  cspinf14 -- CycleSpinning CP De-Noising
  cspinf15 -- Two Sinusoids for Fourier Denoising
  cspinf16 -- Non-Spinning Fourier De-Noising
  cspinf17 -- Cycle-Spinning Fourier De-Noising
  cspinf18 -- Translation-Invariant Table of Blocks
  cspinf19 -- Unscrambling TI Table to Stat Table
  cspinf20 -- Unscrambling TI Table 
  cspinf21 -- Thresholding Translation-Invariant Table of Blocks
  CSFig -- callback implementing functions of SpinDemo
  CSInit -- setup global datastructures for SpinCycle
  CSIntro -- Info for SpinDemo
  FourierDeNoise -- Hard thresholding in Fourier Domain
  SpinCycleInit -- setup global datastructures for SpinCycle
  CSFig -- callback implementing functions of SpinDemo
  WPDeNoise_TI -- Translation-invariant de-noising of object in specific WP basis
 
 
  ******* ./Wavelab850/Papers/Tour/ *******
  CompressoGram -- compare compress abilities of four transforms
  FastHist -- fast histogram tool
  LogoGram -- Time Series Spectrum Estimation
  NoiseMaker -- Add Noise to Signal
  Rootogram -- compute rootogram of counts data
  toufig01 -- Tour Figure 01: Four Spatially Inhomogenoeus Signals
  toufig02 -- Tour Figure 02: Noisy Versions of Four Signals
  toufig03 -- Tour Figure 03: Wavelet Shrinkage of the Four Noisy Signals
  toufig04 -- Tour Figure 04: AutoSpline Reconstructions from Noisy Data
  toufig05 -- Tour Figure 05: AutoTrunc Truncated Fourier Reconstructions
  toufig06 -- Tour Figure 06: WaveShrink of object yBlocks in Haar Basis
  toufig07 -- Tour Figure 07: Comparing Compression Abilities
  toufig08 -- Tour Figure 08: Root-Density example
  toufig09 -- Tour Figure 09: Estimating Time Series Spectrum
  toufig10 -- Tour Figure 10: Noisy Differentiation
  toufig11 -- Tour Figure 11: Noisy Differentiation in Wavelet Domain
  toufig12 -- Tour Figure 12: Noisy Differentiation by Ideal Fourier Damping
  toufig13 -- Tour Figure 13: Noisy Deconvolution
  toufig14 -- Tour Figure 14: Depict Deconvolution in Wavelet Domain
  disp('Sorry -- data not available for figure 15')
  toufig15 -- Tour Figure 15: Display Vaguelette Kernels
  TourCleanup -- Clean up Tour global variables
  TourFig -- callback implementing functions of TourDemo
  TourInit -- setup global datastructures for Tour
  TourIntro -- Info about TourDemo
  WaveCoeffLim -- Plot WT with non-standard limits of x-range
 
 
  ******* ./Wavelab850/Papers/VillardDeLans/ *******
  TreeDeNoise -- Selection of an Optimal Zero-Tree Threshold
  VdLCleanup -- Remove global datastructures for VdL
  vdlfig01 -- Plot of Doppler Signal
  vdlfig02 -- Plot of Meyer wavelet
  vdlfig03 -- 95% wavelet/dct compression of seismic signal
  vdlfig04 -- 95% wavelet/dct compression curves for Ingrid Daubechies
  vdlfig05 -- Use stationary Haar transform to denoise Blocks
  vdlfig06 -- Show WP and CP phase planes for Tweet signal
  vdlfig07 -- Show CW and l^1 entropy phase planes for Werner Sorrows signal
  vdlfig08 -- Show MP phase plane for Tweet signal
  vdlfig09 -- Minimum Entropy Focusing: Evaluate minimum
  vdlfig10 -- Minimum Entropy Focusing: Reconstruct blurry signal
  vdlfig11 -- Use tree-constrained thresholding to denoise Blocks
  VdLFig -- callback implementing functions of VdLDemo
  VdLInit -- setup global data structures for VdLDemo
  VdLIntro -- Info about VdLDemo
 
 
  ******* ./Wavelab850/Pursuit/ *******
  Pursuit:Contents v802 -- Matching Pursuit with Wavelet- and Cosine- Packet Dictionaries
  CPAtomicSynthesis -- Generate signal from list of CP-components
  CPImpulse -- CPAnalysis of an individual Cosine Packet
  cp packet table -- ignored [except for size]
  dcp packet table -- complete cosine packet analysis of
  CPMPursuitTour -- Cosine Packet atomic decomposition
  CPPursuitBF -- Matching Pursuit using Cosine Packet Dictionary with backfitting
  CPPursuit -- Matching Pursuit with Cosine Packet Dictionary
  CPPursuitTour -- Cosine Packet atomic decomposition
  FCPSynthesis -- Synthesize signal from full cosine packet table
  FWPAtomicSynthesis -- Fast Wavelet Packet Synthesis from Atom List
  FWPSynthesis -- Synthesize signal from full wavelet packet table
  ImageAtomicPhase -- Image time-frequency content of atom list
  ix2pkt -- Convert linear index to packet table index
  MakeSynthesisTable -- Convert atom list to packet table
  MonitorPursuit -- Display intermediate Pursuit progress
  pkt2ix -- Convert packet table index to linear index
  PlotAtomicPhase -- Plot time-frequency content of atom list
  PlotSynthesisTable -- Display entries in atomlist as WP or CP table
  WPAtomicSynthesis -- Wavelet Packet Synthesis from Atom List
  dwp packet table -- complete wavelet packet analysis of
  WPImpulse -- WPAnalysis of an individual Wavelet Packet
  WPMPursuitTour -- Wavelet Packet atomic decomposition
  WPPursuitBF -- Matching Pursuit using Wavelet Packet Dictionary with backfitting
  WPPursuit -- Atomic Decomposition into Wavelet Packets by Matching Pursuit
  WPPursuitTour -- Wavelet Packet atomic decomposition
 
 
  ******* ./Wavelab850/TimeFrequency/ *******
  CohenDist -- Choi William Distribution
  TimeFrequency:Contents v850 -- Time-Frequency Distribution Tools
  IWindowFT -- Inverse Window Fourier Transform
  MakeWindow -- Make artificial Window
  Ridge_WindowFT -- Find Ridges of Window Fourier Transform
  Spectrogram -- Time-Varying Spectrum
  TFDist_AF -- Alias-Free Generalized Discrete Time-Frequency Distribution
  WignerDist -- Alias-Free Wigner-Ville Distribution
  WindowFT -- Window Fourier Transform
  WVDist_AF -- Alias-Free Wigner-Ville Distribution
  WVDist -- Wigner-Ville Distribution
 
 
  ******* ./Wavelab850/Utilities/ *******
  AppendTitle -- Utility to Build Title String
  AutoImage -- Image display of object assuming arbitrary values
  Utilities:Contents v850 -- Utilties for Writing Scripts
  CutDyad -- Truncate signal to Dyadic length
  GrayImage -- Standard gray-scale image display
  HitAnyKey -- Tool for pausing in scripts
  ifprint -- Conditional printing to postscript file
  LockAxes -- Version-independent axis command
  MakeTiledFigures -- Tile the screen with figures
  PadDyad -- Zero-fill signal to Dyadic length
  RegisterPlot -- Add legend with file name, date, flag
  rnshift -- t circular right shift of 1-d signal
  ShapeAsRow -- Make signal a row vector
  ShapeLike -- Make 1-d signal with given shape
  UnlockAxes -- Version-independent axis command
  versaplot -- Version-independent plot routine
  WaitUntil -- Burn up CPU cycles until sec seconds elapse from oldclock
  WhiteNoise -- Version-independent white noise generator
 
 
  ******* ./Wavelab850/Workouts/ *******
  Workouts:Contents v850 -- WaveLab Workouts
 
 
  ******* ./Wavelab850/Workouts/BestOrthoBasis/ *******
  BBWork01 -- Illustrate Time-Frequency Tiling by Wavelet Packets
  Panel 1 -- Fourier Transform on side
  Panel 2 -- Time-Frequency Rectangle
  Panel 4 -- Signal
  BBWork02 -- Illustrate Time-Frequency Tiling by Cosine Packets
  Panel 1 -- Fourier Transform on side
  Panel 2 -- Time-Frequency Rectangle
  Panel 4 -- Signal
  BBWork03 -- Ilustrate Tiling <--> Tree Correspondence, WP Dict;
  BBWork04 -- Ilustrate Tiling <--> Tree Correspondence, CP Dict;
  BBWork11 -- Illustrate WP Methodology on LoSine, SlowChirp
  BBWork12 -- Illustrate CP Analysis Techniques on SlowChirp
  BBWork20 -- Illustrate Candidate signals
  BBWork21 -- Wavelet Packet Tour on Four Basic Signals
  BBWork22 -- Wavelet Packet Tour on Four Complex Signals
  BBWork31 -- Cosine Packet Tour on Four Basic Signals
  BBWork32 -- Cosine Packet Tour on Four Complex Signals
  BBWork41 -- Illustrate CP/WP Analysis Techniques on Seismic
  BBWork42 -- Illustrate CP/WP Analysis Techniques on Warble
  BBWork43 -- WP and CP Tours on Transients Signal
  BBWork61 -- Compression Numbers on Four Basic Signals
  Compare4Compress -- Compare WP and CP Compressions
  BBWork01 -- Illustrate Time-Frequency Tiling by Wavelet Packets
  BBWork02 -- Illustrate Time-Frequency Tiling by Cosine Packets
  BBWork03 -- Ilustrate Tiling <--> Tree Correspondence, WP Dict;
  BBWork04 -- Ilustrate Tiling <--> Tree Correspondence,CP Dict;
  BBWork11 -- Illustrate WP Methodology on LoSine, SlowChirp
  BBWork12 -- Illustrate CP Analysis Techniques on SlowChirp
  BBWork21 -- Wavelet Packet Tour on Four Basic Signals
  BBWork22 -- Wavelet Packet Tour on Four Complex Signals
  BBWork31 -- Cosine Packet Tour on Four Basic Signals
  BBWork32 -- Cosine Packet Tour on Four Complex Signals
  BBWork41 -- Illustrate CP/WP Analysis Techniques on Seismic
  BBWork42 -- Illustrate CP/WP Analysis Techniques on Warble
  BBWork61 -- Compression Numbers on Four Basic Signals
  BBWorkout:Contents -- BestBasis Workout
 
 
  ******* ./Wavelab850/Workouts/MatchingPursuit/ *******
  Compare4Decomp -- Compare four ways of decomposing object
  MPWork00 -- Basics of Matching Pursuit -- Illustrate Step-by-Step Behavior
  MPWork21 -- Matching-Pursuit, Wavelet Packet Tour on Four Basic Signals
  MPWork22 -- Matching-Pursuit, Wavelet Packet Tour on Four Complex Signals
  MPWork31 -- Matching-Pursuit, Cosine Packet Tour on Four Basic Signals
  MPWork32 -- Matching-Pursuit, Cosine Packet Tour on Four Complex Signals
  MPWork61 -- Matching-Pursuit, Best-Basis Compression # Comparison on Four Basic Signals
  MPWork62 -- Matching-Pursuit, Best-Basis Compression # Comparison on Four Complex Signals
  MPWorkout:Contents -- Matching Pursuit Workout
  CPPrintPursuit -- Print Intermediate Results, MP with Cosine Packet Dictionary
  MPWork00 -- Basics of Matching Pursuit, Illustrate Step-by-Step Behavior
  MPWork21 -- Matching-Pursuit, Wavelet Packet Tour on Four Basic Signals
  MPWork22 -- Matching-Pursuit, Wavelet Packet Tour on Four Complex Signals
  MPWork31 -- Matching-Pursuit, Cosine Packet Tour on Four Basic Signals
  MPWork32 -- Matching-Pursuit, Cosine Packet Tour on Four Complex Signals
  MPWork43 -- Matching-Pursuit, WP and CP Tours on Transients Signal
  MPWork61 -- Matching-Pursuit, Best-Basis Compression
  MPWork62 -- Matching-Pursuit, Best-Basis Compression
  PlotEquivDimension -- Plot Equivalent Dimensions
  PlotEquivRate -- Plot Equivalent Rates
  PlotSemiCompressNum -- Plot Compression Numbers on Semilog Scale
  WPPrintPursuit -- Print Intermediate Results, MP with Wavelet Packets
 
 
  ******* ./Wavelab850/Workouts/MultiFractal/ *******
  MFWork01 -- CWT Analysis of Cantor Signal
  MFWork02 -- CWT Analysis of Brownian Signal
  MultiFractal:Contents -- Continuous Wavelet Transform Workouts
 
 
  ******* ./Wavelab850/Workouts/Toons/ *******
  BB_CP -- Decomposition in Best Cosine Packet Basis
  BB_WP -- Decomposition in Best Wavelet Packet Basis
  4 subplots -- several approximations Best-m toon0524
  4 subplots -- several approximations First-m toon0523
  Toons:Contents v850 -- Cartoon Guide to Wavelets
  DoAllToons -- Run all toons scripts. 
  Usefulness -- run this, beginning to end, to see if all toons work
  KeepBiggest -- Kill (set to zero) all but N biggest elements
  KeepFirst -- Kill (set to zero) all but N first elements
  PlotLevelNorms -- Plot level-wise norms
  PlotOrderCoef -- Plot Ordered Coefficients
  PrintAllToons -- Print all toons, in sequence
  Coiflet C3 -- orthonormal wavelets system where both father and
  Daubechies D4 -- the first continuous, compactly supported
  Haar -- the first wavelet; a square-wave wavelet
  Symmlet S8 -- nearly-symmetric orthogonal wavelet of
  toon0111 -- Wavelet Families
  toon0112 -- Interpolating Wavelets
  toon0113 -- Average Interpolating Wavelets
  toon0114 -- Meyer Wavelets
  toon0121 -- Wavelets Come in Genders
  toon0131 -- Scale Families of Wavelets
  toon0132 -- Wavelets come at all different scales and positions
  toon0140 -- Illustrating Boundary Wavelets
  toon0141 -- Illustrating Boundary Wavelets
  toon0142 -- Illustrating Boundary Wavelets
  toon0151 -- Visualize wavelet decomposition of Ramp
  toon0152 -- Visualize wavelet decomposition of Doppler
  toon0161 -- Visualize multi-resolution decomposition
  toon0171 -- Illustrate smoothness of wavelets
  toon0172 -- Illustrate smoothness of wavelets
  toon0173 -- Illustrate smoothness of wavelets
  toon0174 -- Illustrate smoothness of wavelets
  toon0181 -- Scale Families of Wavelets
  toon0182 -- Scale Families of Wavelets
  toon0183 -- Scale Families of Wavelets
  toon0184 -- Scale Families of Wavelets
  toon0211 -- 2d Wavelet Families
  toon0212 -- 2d Wavelet Families
  toon0221 -- 2d Wavelet Families
  toon0222 -- 2d Wavelet Families
  toon0231 -- Visualize Wavelet Xform
  toon0232 -- Visualize Wavelet Xform
  toon0233 -- Visualize Wavelet Xform
  toon0311 -- Analysis of Sinusoid with Haar Wavelet
  toon0312 -- Analysis of Sinusoid with S8 Wavelet
  toon0331 -- Behavior of Wavelet Transform of Singularities
  toon0332 -- Behavior of Wavelet Transform of Singularities
  toon0333 -- Behavior of Wavelet Transform of Singularities
  toon0334 -- Behavior of Wavelet Transform of Singularities
  toon0335 -- Behavior of Wavelet Transform of Singularities
  toon0341 -- Compare Symmlet & Haar Coefficient Decay
  toon0342 -- Compare Wavelet & Haar Coefficient Decay
  toon0411 -- Compare First-n and Best-n Reconstructions
  toon0412 -- Compare First-n and Best-n Reconstructions
  toon0421 -- Compare Wavelet & Haar Partial Reconstructions
  toon0422 -- Compare Wavelet & Haar Partial Reconstructions
  toon0431 -- Compare First-n and Best-n Compression Numbers
  toon0432 -- First-m Compression Numbers of Four Objects
  toon0433 -- Best-m Compression Numbers of Four Objects
  toon0441 -- Best-m Wavelet & DCT Partial Reconstructions
  toon0442 -- First-m Wavelet & DCT Partial Reconstructions
  toon0443 -- Best-m Compression Numbers in DCT & DWT of Four Objects
  toon0444 -- First-m Compression Numbers in DCT & DWT of Four Objects
  toon0451 -- Unconditional Basis Property
  toon0511 -- Lincoln Refinement Scheme
  toon0512 -- Lincoln Refinement Scheme
  toon0513 -- Lincoln Refinement Scheme
  toon0514 -- Lincoln Refinement Scheme
  toon0521 -- Progressive Transmission
  toon0522 -- Progressive Transmission
  toon0523 -- Progressive Transmission
  toon0524 -- Progressive Transmission
  toon0531 -- Compress/Decompress Seismic Trace
  toon0532 -- Compress/Decompress Seismic Trace
  toon0533 -- Compare all four methods of analysis on seismic signal
  toon0541 -- Data Compression of Ingrid Daubechies
  toon0542 -- Data Compression of Ingrid Daubechies
  toon0543 -- Data Compression of Ingrid Daubechies
  toon0544 -- Data Compression of Ingrid Daubechies
  toon0545 -- Data Compression of Ingrid Daubechies
  toon0546 -- Data Compression of Ingrid Daubechies
  toon0547 -- Data Compression of Ingrid Daubechies
  toon0548 -- Data Compression of Ingrid Daubechies
  toon0551 -- Four Spatially Inhomogenoeus Signals
  toon0552 -- Noisy Versions of Four Signals
  toon0553 -- Wavelet Shrinkage of the Four Noisy Signals
  toon0554 -- AutoSpline Reconstructions from Noisy Data
  toon0561 -- Fast Algorithm for Matrix Multiplication
  toon0562 -- Fast Algorithm for Matrix Multiplication
  toon0563 -- Fast Algorithm for Matrix Multiplication
  toon0564 -- Fast Algorithm for Matrix Multiplication
  toon0565 -- Fast Algorithm for Matrix Multiplication
  toon0711 -- Show some Wavelet Packets
  toon0721 -- Show some Wavelet Packets in Frequency Domain
  Panel 1 -- Fourier Transform on side
  Panel 2 -- Time-Frequency Rectangle
  Panel 4 -- Signal
  toon0731 -- Illustrate Time-Frequency Tiling by Wavelet Packets
  Panel 1 -- Fourier Transform on side
  Panel 2 -- Time-Frequency Rectangle
  Panel 4 -- Signal
  toon0732 -- Illustrate Time-Frequency Tiling by Wavelet Packets
  toon0741 -- Illustrate WP Methodology on SlowChirp
  toon0742 -- Visualize wavelet packet decomposition of SlowChirp
  toon0811 -- Show some Cosine Packets
  toon0821 -- Some Cosine Packets in the Frequency Domain
  Panel 1 -- Fourier Transform on side
  Panel 2 -- Time-Frequency Rectangle
  Panel 4 -- Signal
  toon0831 -- Illustrate Time-Frequency Tiling by Cosine Packets
  Panel 1 -- Fourier Transform on side
  Panel 2 -- Time-Frequency Rectangle
  Panel 4 -- Signal
  toon0832 -- Illustrate Time-Frequency Tiling by Cosine Packets
  toon0841 -- Illustrate CP Methodology on SlowChirp
  toon0842 -- Visualize cosine packet decomposition of SlowChirp
  toon1001 -- Wavelet Packet Analysis on LoSine
  toon1002 -- Wavelet Packet Analysis on Ramp
  toon1003 -- Wavelet Packet Analysis on Doppler
  toon1004 -- Wavelet Packet Analysis on Dirac
  toon1005 -- Wavelet Packet Analysis on HiSine
  toon1006 -- Wavelet Packet Analysis on LinChirp
  toon1007 -- Wavelet Packet Analysis on QuadChirp
  toon1101 -- Cosine Packet Analysis on LoSine
  toon1102 -- Cosine Packet Analysis on Ramp
  toon1103 -- Cosine Packet Analysis on Doppler
  toon1104 -- Cosine Packet Analysis on Dirac
  toon1105 -- Cosine Packet Analysis on HiSine
  toon1106 -- Cosine Packet Analysis on LinChirp
  toon1107 -- Cosine Packet Analysis on QuadChirp
  toon1611 -- FBI Fingerprint Compression
  toon1612 -- FBI Fingerprint Compression
  toon1613 -- FBI Fingerprint Compression
