help cleanup before workshop

Kevin2599 · Oct 2, 2012 · 4f22266 · 4f22266
1 parent 914550e
commit 4f22266
Show file tree

Hide file tree

Showing 16 changed files with 263 additions and 143 deletions.
diff --git a/comp/comp_window.m b/comp/comp_window.m
@@ -50,7 +50,7 @@
     if isrect
       g=gabdual(g,a,M);
     else
-      g=gabdual(g,a,M,[],lt);
+      g=gabdual(g,a,M,'lt',lt);
     end;
     info.isdual=1;
     info.tfr=a*M/L;
@@ -59,7 +59,7 @@
     if isrect
       g=gabtight(a,M,L);
     else
-      g=gabtight(a,M,L,lt);
+      g=gabtight(a,M,L,'lt',lt);
     end;
     info.tfr=a*M/L;
     info.istight=1;
@@ -94,15 +94,15 @@
     if isrect
       g = gabdual(g,a,M,L);
     else
-      g = gabdual(g,a,M,L,lt);
+      g = gabdual(g,a,M,L,'lt',lt);
     end;
     info.isdual=1;
    case {'tight'}
     [g,info.auxinfo] = comp_window(g{2},a,M,L,lt,callfun);    
     if isrect
       g = gabtight(g,a,M,L);
     else
-      g = gabtight(g,a,M,L,lt);
+      g = gabtight(g,a,M,L,'lt',lt);
     end;
     info.istight=1;
    case firwinnames

diff --git a/demos/demo_audiocompression.m b/demos/demo_audiocompression.m
@@ -1,17 +1,19 @@
 %DEMO_AUDIOCOMPRESSION  Audio compression using N-term approx
 %
 %   This demos shows how to do audio compression using best N-term
-%   approximation of WMDCT transform.
+%   approximation of an |wmdct|_ transform.
 %
-%   The signal is transformed using an orthonormal WMDCT transform.
-%   Then approximations with fixed number N of coefficients are obtained
+%   The signal is transformed using an orthonormal |wmdct|_ transform.
+%   Then approximations with a fixed number *N* of coefficients are obtained
 %   by:
 %
-%     * Linear approximation     - N coefficients with lowest frequency index
+%     * Linear approximation: The *N* coefficients with lowest frequency
+%       index are kept.
 %
-%     * Non-linear approximation - N largest coefficients (in magnitude)
+%     * Non-linear approximation: The *N* largest coefficients (in
+%     magnitude) are kept.
 %
-%   The corresponding approximated signal is computed with inverse WMDCT.
+%   The corresponding approximated signal can be computed using |iwmdct|_.
 %
 %   .. figure::
 %
@@ -20,39 +22,41 @@
 %      The figure shows the output Signal to Noise Ratio (SNR) as a function
 %      of the number of retained coefficients.
 %
-%      Note: actually, inverse WMDCT is not used, as Parseval theorem states
-%      that the norm of a signal equals the norm of the sequence of its
-%      wmdct coefficients. The latter is used for computing SNRs.
+%   Note: The inverse WMDCT is not needed for computing computing
+%   SNRs. Instead Parseval theorem states that the norm of a signal equals
+%   the norm of the sequence of its |wmdct|_ coefficients.
 
 % Load audio signal
 % Use the 'glockenspiel' signal.
 sig=gspi;
 
 % Shorten signal
-SigLength = 2^16;
-sig = sig(1:SigLength);
+L = 2^16;
+sig = sig(1:L);
 
-% Initializations
-NbFreqBands = 1024;
-NbTimeSteps = SigLength/NbFreqBands;
+% Number of frequency channels
+M = 1024;
+
+% Number of time steps
+N = L/M;
 
 % Generate window
-gamma = wilorth(NbFreqBands,SigLength);
+gamma = wilorth(M,L);
 
 % Compute wmdct coefficients
-c = wmdct(sig,gamma,NbFreqBands);
+c = wmdct(sig,gamma,M);
 
 
 % L2 norm of signal
 InputL2Norm = norm(c,'fro');
 
 % Approximate, and compute SNR values
-kmax = NbFreqBands;
+kmax = M;
 kmin = kmax/32;             % 32 is an arbitrary choice
 krange = kmin:32:(kmax-1);  % same remark
 
 for k = krange,
-    ResL2Norm_NL = norm(c-largestn(c,k*NbTimeSteps),'fro');
+    ResL2Norm_NL = norm(c-largestn(c,k*N),'fro');
     SNR_NL(k) = 20*log10(InputL2Norm/ResL2Norm_NL);
     ResL2Norm_L = norm(c(k:kmax,:),'fro');
     SNR_L(k) = 20*log10(InputL2Norm/ResL2Norm_L);
@@ -63,8 +67,8 @@
 figure(1);
 
 set(gca,'fontsize',14);
-plot(krange*NbTimeSteps,SNR_NL(krange),'x-b',...
-     krange*NbTimeSteps,SNR_L(krange),'o-r');
+plot(krange*N,SNR_NL(krange),'x-b',...
+     krange*N,SNR_L(krange),'o-r');
 axis tight; grid;
 legend('Best N-term','Linear');
 xlabel('Number of Samples', 'fontsize',14);

diff --git a/demos/demo_auditoryfilterbank.m b/demos/demo_auditoryfilterbank.m
@@ -8,7 +8,7 @@
 %   Each channel is subsampled by a factor of 8 (a=8), and to generate a
 %   nice plot, 4 channels per Erb has been used.
 %
-%   The filterbank cover only the positive frequencies, so we must use
+%   The filterbank covers only the positive frequencies, so we must use
 %   |filterbankrealdual|_ and |filterbankrealbounds|_.
 %
 %   .. figure:: 
@@ -26,22 +26,14 @@
 %      gammatone filters on an Erb scale.  The dynamic range has been set to
 %      50 dB, to highlight the most important features.
 %
-%   .. figure:: 
-%
-%      Auditory filterbank representation using Gaussian filters
-%
-%      Auditory filterbank representation of the spoken sentense using
-%      Gaussian filters on an Erb scale.  The dynamic range has been set to
-%      50 dB, to highlight the most important features.
-%
 %   See also: freqtoaud, audfiltbw, gammatonefir, ufilterbank, filterbankrealdual
 %
 %   References: glasberg1990daf
 
 
 % Use part of the 'cocktailparty' spoken sentense.
 f=cocktailparty;
-f=f(1:64000,:);
+f=f(20000:80000,:);
 fs=44100;
 a=8;
 channels_per_erb=4;
@@ -86,34 +78,3 @@
 
 figure(2);
 plotfilterbank(coef_gam,a,fc,fs,dynrange_for_plotting,'audtick');
-
-
-%% ---------------  Gauss filters ----------------------------
-
-
-g_gauss=cell(M,1);
-
-% Account for the bandwidth specification to PGAUSS is for 79% of the ERB(integral).  
-bw_gauss=audfiltbw(fc)/0.79*1.5;
-
-for m=1:M
-  g_gauss{m}=pgauss(filterlength,'fs',fs,'bw',bw_gauss(m),'cf',fc(m));
-end;
-
-disp('Frame bound ratio for Gaussian filterbank, should be close to 1 if the filters are choosen correctly.');
-filterbankrealbounds(g_gauss,a,filterlength)
-
-% Synthesis filters
-gd_gauss=filterbankrealdual(g_gauss,a,filterlength);
-
-% Analysis transform
-coef_gauss=ufilterbank(f,g_gauss,a);
-
-% Synthesis transform
-r_gauss=2*real(ifilterbank(coef_gauss,gd_gauss,a));
-
-disp('Error in reconstruction, should be close to zero.');
-norm(f-r_gauss)/norm(f)
-
-figure(3);
-plotfilterbank(coef_gauss,a,fc,fs,dynrange_for_plotting,'audtick');
diff --git a/filterbank/filterbank.m b/filterbank/filterbank.m
@@ -14,9 +14,22 @@
 %   The output coefficients are stored a cell array. More precisely, the
 %   n'th cell of *c*, `c{m}`, is a 2D matrix of size $M(n) \times W$ and
 %   containing the output from the m'th channel subsampled at a rate of
-%   $a(m)$.  `c{m}(n,l)` is thus the value of the coefficient for time index
+%   $a(m)$.  *c\{m\}(n,l)* is thus the value of the coefficient for time index
 %   *n*, frequency index *m* and signal channel *l*.
 %
+%   The coefficients *c* computed from the signal *f* and the filterbank
+%   with windows *g_m* are defined by
+%
+%   ..            L-1
+%      c_m(n+1) = sum f(l+1) * g_m (a(m)n-l+1)
+%                 l=0
+%
+%   .. math:: c_m\left(n+1\right)=\sum_{l=0}^{L-1}f\left(l+1\right)g\left(a_mn-l+1\right)
+%
+%   where $an-l$ is computed modulo $L$.
+%
+%   See also: ufilterbank, ifilterbank, pfilt
+%
 %   References: bohlfe02
 
 if nargin<3

diff --git a/filterbank/ufilterbank.m b/filterbank/ufilterbank.m
@@ -14,6 +14,16 @@
 %   column in *f* is filtered by the corresponding filter in *g*. If *f* is
 %   a matrix, the output will be 3-dimensional, and the third dimension will
 %   correspond to the columns of the input signal.
+%
+%   The coefficients *c* computed from the signal *f* and the filterbank
+%   with windows *g_m* are defined by
+%
+%   ..              L-1
+%      c(n+1,m+1) = sum f(l+1) * g_m (an-l+1)
+%                   l=0
+%
+%   .. math:: c\left(n+1,m+1\right)=\sum_{l=0}^{L-1}f\left(l+1\right)g\left(an-l+1\right)
+
 %
 %   See also: ifilterbank, filterbankdual
 %

diff --git a/fourier/pfilt.m b/fourier/pfilt.m
@@ -12,6 +12,17 @@
 %   `pfilt(f,g,a,dim)` filters along dimension dim. The default value of
 %   [] means to filter along the first non-singleton dimension.
 %
+%   The coefficients obtained from filtering a signal *f* by a filter *g* are
+%   defined by
+%
+%   ..          L-1
+%      c(n+1) = sum f(l+1) * g(an-l+1)
+%               l=0
+%
+%   .. math:: c\left(n+1\right)=\sum_{l=0}^{L-1}f\left(l+1\right)g\left(an-l+1\right)
+%
+%   where $an-l$ is computed modulo $L$.
+%
 %   See also: pconv
 
 

diff --git a/frames/frsyniter.m b/frames/frsyniter.m
@@ -1,5 +1,5 @@
 function [f,relres,iter]=frsyniter(F,c,varargin)
-%FRSYNITER  Iterative filter bank inversion
+%FRSYNITER  Iterative analysis frame inversion
 %   Usage:  f=frsyniter(F,c);
 %
 %   Input parameters:
@@ -29,45 +29,94 @@
 %
 %     'maxit',n    Do at most n iterations.
 %
+%     'unlocbox'   Use unlocbox. This is the default.
+%
+%     'lsqr'       Use the Matlab `lsqr`function.
+%
+%     'print'      Display the progress.
+%
+%     'quiet'      Don't print anything, this is the default.
+%
 %   Examples
 %   --------
 %
 %   The following example shows how to rectruct a signal without ever
 %   using the dual frame:::
 %
-%      F=newframe('dgtreal','gauss','dual',10,20);
+%      F=newframe('dgtreal','gauss','none',10,20);
 %      c=frana(F,bat);
 %      r=frsyniter(F,c);
 %      norm(bat-r)/norm(bat)
 %
 %   See also: newframe, frana, frsyn
 
+% AUTHORS: Nathanael Perraudin and Peter L. Soendergaard
+
   if nargin<2
     error('%s: Too few input parameters.',upper(mfilename));
   end;
 
   definput.keyvals.Ls=[];
   definput.keyvals.tol=1e-6;
   definput.keyvals.maxit=100;
-
+  definput.flags.alg={'unlocbox','lsqr'};
+  definput.keyvals.printstep=10;
+  definput.flags.print={'quiet','print'};
+
   [flags,kv,Ls]=ltfatarghelper({'Ls'},definput,varargin);
 
   % Determine L from the first vector, it must match for all of them.
   L=framelengthcoef(F,size(c,1));
 
-  % Set up the persisten variable
-  afun(1, 'dummy', F);
 
-  % It is possible to specify the initial guess
-  [f,flag,~,iter,relres]=lsqr(@afun,c,kv.tol,kv.maxit);
-
+  if flags.do_lsqr
+      % Set up the persisten variable
+      afun(1, 'dummy', F);
+
+      % It is possible to specify the initial guess
+      [f,flag,~,iter,relres]=lsqr(@afun,c,kv.tol,kv.maxit);
+
+  end;
+
+  if flags.do_unlocbox
+
+      % Get the upper frame bound (Or an estimation bigger than the bound)
+      [~,B]=framebounds(F,L,'a'); 
+
+      % Set the parameter for the fast projection on a B2 ball
+      param.At=@(x) franaadj(F,x);  % adjoint operator
+      param.A=@(x)  frana(F,x);     % direct operator
+      param.y=c;                    % coefficient
+      param.tight=0;                % It's not a tight frame
+      param.max_iter=kv.maxit;
+      param.tol=kv.tol; 
+      param.nu=B;
+
+      % Display parameter 0 nothing, 1 summary at convergence, 2 all
+      % steps
+      if flags.do_print
+          param.verbose=1;
+      else
+          param.verbose=0;
+      end;
+
+      % Make the projection. Requires UNLocBOX
+      [f, ~] = fast_proj_B2(zeros(L,1), 0, param);
+
+      % compute the residue
+      res = param.A(f) - param.y; norm_res = norm(res(:), 2);
+      relres=norm_res/norm(c(:), 2);
+
+      iter=0; % The code of the fast_proj_B2 is not yet compatible with this
+  end;
+
   % Cut or extend f to the correct length, if desired.
   if ~isempty(Ls)
-    f=postpad(f,Ls);
+      f=postpad(f,Ls);
   else
-    Ls=L;
+      Ls=L;
   end;
-
+      
 end
 
 % This is a nested function, as it must use variables from the scope