Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Commit [f86e6c] Maximize Restore History

nonsepidgt and idgt functionality integrated into standard dgt and idgt

Peter L. Soendergaard Peter L. Soendergaard 2012-09-18

changed comp/comp_window.m
changed comp/gabpars_from_window.m
changed comp/gabpars_from_windowsignal.m
changed gabor/Contents.m
changed gabor/dgt.m
changed gabor/dgtlength.m
changed gabor/gabwin.m
changed gabor/idgt.m
changed gabor/idgtreal.m
changed gabor/wilwin.m
changed nonsepgab/Contents.m
changed testing/test_nonsepdgt.m
copied nonsepgab/inonsepdgt.m -> comp/comp_inonsepdgt.m
copied nonsepgab/latticetype2matrix.m -> gabor/latticetype2matrix.m
copied nonsepgab/matrix2latticetype.m -> gabor/matrix2latticetype.m
copied nonsepgab/nonsepdgt.m -> comp/comp_nonsepdgt.m
comp/comp_window.m Diff Switch to side-by-side view
Loading...
comp/gabpars_from_window.m Diff Switch to side-by-side view
Loading...
comp/gabpars_from_windowsignal.m Diff Switch to side-by-side view
Loading...
gabor/Contents.m Diff Switch to side-by-side view
Loading...
gabor/dgt.m Diff Switch to side-by-side view
Loading...
gabor/dgtlength.m Diff Switch to side-by-side view
Loading...
gabor/gabwin.m Diff Switch to side-by-side view
Loading...
gabor/idgt.m Diff Switch to side-by-side view
Loading...
gabor/idgtreal.m Diff Switch to side-by-side view
Loading...
gabor/wilwin.m Diff Switch to side-by-side view
Loading...
nonsepgab/Contents.m Diff Switch to side-by-side view
Loading...
testing/test_nonsepdgt.m Diff Switch to side-by-side view
Loading...
nonsepgab/inonsepdgt.m to comp/comp_inonsepdgt.m
--- a/nonsepgab/inonsepdgt.m
+++ b/comp/comp_inonsepdgt.m
@@ -1,5 +1,5 @@
-function [f,g]=inonsepdgt(coef,g,a,lt,varargin)
-%INONSEPDGT  Inverse discrete Gabor transform
+function f=comp_inonsepdgt(coef,g,a,lt,do_timeinv,alg)
+%COMP_INONSEPDGT  Compute Inverse discrete Gabor transform
 %   Usage:  f=inonsepdgt(c,g,a,lt);
 %           f=inonsepdgt(c,g,a,lt,Ls);
 %
@@ -8,7 +8,8 @@
 %         g     : Window function.
 %         a     : Length of time shift.
 %         lt    : Lattice type
-%         Ls    : length of signal.
+%         do_timeinv : Do a time invariant phase ?
+%         alg   : Choose algorithm
 %   Output parameters:
 %         f     : Signal.
 %
@@ -17,42 +18,13 @@
 %   lattice type *lt*. The number of channels is deduced from the size of
 %   the coefficients *c*.
 %
-%   `inonsepdgt(c,g,a,lt,Ls)` does as above but cuts or extends *f* to length
-%   *Ls*.
+%     * $alg=0$ : Choose the fastest algorithm
 %
-%   `[f,g]=inonsepdgt(...)` additionally outputs the window used in the
-%   transform. This is useful if the window was generated from a description
-%   in a string or cell array.
+%     * $alg=0$ : Always choose multi-win
 %
-%   For perfect reconstruction, the window used must be a dual window of the
-%   one used to generate the coefficients.
+%     * $alg=1$ : Always choose shear
 %
-%   The window *g* may be a vector of numerical values, a text string or a
-%   cell array. See the help of |gabwin|_ for more details.
-%
-%   If *g* is a row vector, then the output will also be a row vector. If *c* is
-%   3-dimensional, then `inonsepdgt` will return a matrix consisting of one column
-%   vector for each of the TF-planes in *c*.
-%
-%   Assume that `f=inonsepdgt(c,g,a,lt,L)` for an array *c* of size $M\times N$.
-%   Then the following holds for $k=0,\ldots,L-1$:
-% 
-%   ..          N-1 M-1          
-%     f(l+1)  = sum sum c(m+1,n+1)*exp(2*pi*i*m*l/M)*g(l-a*n+1)
-%               n=0 m=0          
-%
-%   .. math:: f(l+1) = \sum_{n=0}^{N-1}\sum_{m=0}^{M-1}c(m+1,n+1)e^{2\pi iml/M}g(l-an+1)
-%
-%   `inonsepdgt` takes the following flags at the end of the line of input
-%   arguments:
-%
-%     'freqinv'  Compute an `inonsepdgt` using a frequency-invariant phase. This
-%                is the default convention described above.
-%
-%     'timeinv'  Compute an `inonsepdgt` using a time-invariant phase. This
-%                convention is typically used in filter bank algorithms.
-%
-%   See also:  dgt, gabwin, dwilt, gabtight
+%   This is a computational subroutine, do not call it directly.
 
 %   AUTHOR : Nicki Holighaus and Peter Soendergaard
 %   TESTING: TEST_NONSEPDGT
@@ -60,34 +32,14 @@
 
 % Check input paramameters.
 
-if nargin<3
-  error('%s: Too few input parameters.',upper(mfilename));
-end;
-
-if numel(g)==1
-  error('g must be a vector (you probably forgot to supply the window function as input parameter.)');
-end;
-
-definput.keyvals.Ls=[];
-definput.flags.nsalg={'multiwin','shear'};
-[flags,kv,Ls]=ltfatarghelper({'Ls'},definput,varargin);
-
-wasrow=0;
 
 M=size(coef,1);
 N=size(coef,2);
 W=size(coef,3);
-
-% use assert_squarelat to check a and the window size.
-assert_squarelat(a,M,1,'INONSEPDGT');
-
 L=N*a;
 
-[g,info] = comp_window(g,a,M,L,lt,'INONSEPDGT');
-
-assert_L(L,size(g,1),L,a,M,'INONSEPDGT');
-
-if flags.do_multiwin
+if (alg==1) || (alg==0 && lt(2)<=2) 
+    
     % ----- algorithm starts here, split into sub-lattices ---------------
     
     mwin=comp_nonsepwin2multi(g,a,M,lt);
@@ -110,11 +62,8 @@
         f=f+comp_idgt(sub,mwin(:,ii+1),lt(2)*a,M,L,0);  
     end;
 
-end;
-
-
-if flags.do_shear
-
+else
+    
     b = L/M;
     s = b*lt(1)/lt(2);
     [s0,s1,br] = shearfind(a,b,s,L);
@@ -166,12 +115,3 @@
 
 end;
     
-% Cut or extend f to the correct length, if desired.
-if ~isempty(Ls)
-  f=postpad(f,Ls);
-else
-  Ls=L;
-end;
-
-f=comp_sigreshape_post(f,Ls,wasrow,[0; W]);
-
nonsepgab/matrix2latticetype.m to gabor/matrix2latticetype.m
--- a/nonsepgab/matrix2latticetype.m
+++ b/gabor/matrix2latticetype.m
@@ -6,6 +6,14 @@
 %   description into the standard description of a lattice using the *a*,
 %   *M* and *lt*. The conversion is *only* valid for the specified transform
 %   length *L*.
+%
+%   The lattice type *lt* is a $1 \times 2$ vector $[lt_1,lt_2]$ denoting an
+%   irreducible fraction $lt_1/lt_2$. This fraction describes the distance
+%   in frequency (counted in frequency channels) that each coefficient is
+%   offset when moving in time by the time-shift of *a*. Some examples:
+%   `lt=[0 1]` defines a square lattice, `lt=[1 2]` defines the quinqux
+%   (almost hexagonal) lattice, `lt=[1 3]` describes a lattice with a
+%   $1/3$ frequency offset for each time shift and so forth.
 %
 %   An example:::
 %
@@ -17,9 +25,9 @@
 %   The following code generates plots which show the coefficient layout
 %   and enumeration of the first 4 lattices in the time-frequecy plane:::
 %
-%     a=4;
+%     a=6;
 %     M=6;
-%     L=72;
+%     L=36;
 %     b=L/M;
 %     N=L/a;
 %     cw=3;
@@ -44,7 +52,6 @@
 %       axis('square');
 %       title(sprintf('lt=[%i %i]',lt1(fignum),lt2(fignum)),'Fontsize',ftz);
 %     end;
-
 %
 %   See also: latticetype2matrix
 
nonsepgab/nonsepdgt.m to comp/comp_nonsepdgt.m
--- a/nonsepgab/nonsepdgt.m
+++ b/comp/comp_nonsepdgt.m
@@ -1,110 +1,29 @@
-function [c,Ls,g]=nonsepdgt(f,g,a,M,lt,varargin)
-%NONSEPDGT  Non-separable Discrete Gabor transform
+function c=comp_nonsepdgt(f,g,a,M,lt,do_timeinv,alg)
+%NONSEPDGT  Compute Non-separable Discrete Gabor transform
 %   Usage:  c=nonsepdgt(f,g,a,M,lt);
-%           c=nonsepdgt(f,g,a,M,lt,L);
-%           [c,Ls]=nonsepdgt(f,g,a,M,lt);
-%           [c,Ls]=nonsepdgt(f,g,a,M,lt,L);
 %
 %   Input parameters:
 %         f     : Input data.
 %         g     : Window function.
 %         a     : Length of time shift.
+%         M     : Number of channels.
 %         lt    : Lattice type
-%         M     : Number of channels.
-%         L     : Length of transform to do.
+%         do_timeinv : Do a time invariant phase ?
+%         alg   : Choose algorithm
 %   Output parameters:
 %         c     : $M \times N$ array of coefficients.
-%         Ls    : Length of input signal.
 %
 %   `nonsepdgt(f,g,a,M,lt)` computes the non-separable discrete Gabor
 %   transform of the input signal *f* with respect to the window *g*,
-%   time-shift *a*, number of channels *M* and lattice type *lt*. The output is a
-%   vector/matrix in a rectangular layout.
+%   time-shift *a*, number of channels *M* and lattice type *lt*.
 %
-%   The lattice type *lt* is a $1 \times 2$ vector $[lt_1,lt_2]$ denoting an
-%   irreducible fraction $lt_1/lt_2$. This fraction describes the distance
-%   in frequency (counted in frequency channels) that each coefficient is
-%   offset when moving in time by the time-shift of *a*. Some examples:
-%   `lt=[0 1]` defines a square lattice, `lt=[1 2]` defines the quinqux
-%   (almost hexagonal) lattice, `lt=[1 3]` describes a lattice with a
-%   $1/3$ frequency offset for each time shift and so forth.
+%     * $alg=0$ : Choose the fastest algorithm
 %
-%   The length of the transform will be the smallest multiple of *a* and *M*
-%   that is larger than the signal. *f* will be zero-extended to the length of
-%   the transform. If *f* is a matrix, the transformation is applied to each
-%   column. The length of the transform done can be obtained by
-%   `L=size(c,2)*a;`
+%     * $alg=0$ : Always choose multi-win
 %
-%   The window *g* may be a vector of numerical values, a text string or a
-%   cell array. See the help of |gabwin|_ for more details.
+%     * $alg=1$ : Always choose shear
 %
-%   `nonsepdgt(f,g,a,M,lt,L)` computes the Gabor coefficients as above, but
-%   does a transform of length *L*. f will be cut or zero-extended to length
-%   *L* before the transform is done.
-%
-%   `[c,Ls]=nonsepdgt(...)` additionally returns the length of the input
-%   signal *f*. This is handy for reconstruction::
-%
-%               [c,Ls]=nonsepdgt(f,g,a,M,lt);
-%               fr=inonsepdgt(c,gd,a,Ls);
-%
-%   will reconstruct the signal *f* no matter what the length of *f* is,
-%   provided that *gd* is a dual window of *g*.
-%
-%   `[c,Ls,g]=nonsepdgt(...)` additionally outputs the window used in the
-%   transform. This is useful if the window was generated from a description
-%   in a string or cell array.
-%
-%   The non-separable discrete Gabor transform is defined as follows:
-%   Consider a window *g* and a one-dimensional signal *f* of length *L* and
-%   define $N=L/a$.  The output from `c=nonsepdgt(f,g,a,M,lt)` is then given
-%   by:
-%
-%   ..              L-1 
-%      c(m+1,n+1) = sum f(l+1)*conj(g(l-a*n+1))*exp(-2*pi*i*(m+w(n))*l/M), 
-%                   l=0  
-%
-%   .. math:: c\left(m+1,n+1\right)=\sum_{l=0}^{L-1}f(l+1)\overline{g(l-an+1)}e^{-2\pi il(m+w(n))/M}
-%
-%   where $m=0,\ldots,M-1$ and $n=0,\ldots,N-1$ and $l-an$ is computed
-%   modulo *L*.  The additional offset $w$ is given by
-%   $w(n)=\mod(n\cdot lt_1,lt_2)/lt_2$ in the formula above.
-%
-%   Coefficient layout:
-%   -------------------
-%
-%   The following code generates plots which show the coefficient layout
-%   and enumeration of the first 4 lattices in the time-frequecy plane:::
-%
-%     a=6;
-%     M=6;
-%     L=36;
-%     b=L/M;
-%     N=L/a;
-%     cw=3;
-%     ftz=12;
-%     
-%     [x,y]=meshgrid(a*(0:N-1),b*(0:M-1));
-%
-%     lt1=[0 1 1 2];
-%     lt2=[1 2 3 3];
-%
-%     for fignum=1:4
-%       figure;
-%       z=y;
-%       if lt2(fignum)>0
-%         z=z+mod(lt1(fignum)*x/lt2(fignum),b);
-%       end;
-%       for ii=1:M*N
-%         text(x(ii)-cw/4,z(ii),sprintf('%2.0i',ii),'Fontsize',ftz);
-%         rectangle('Curvature',[1 1], 'Position',[x(ii)-cw/2,z(ii)-cw/2,cw,cw]);
-%       end;
-%       axis([-cw L -cw L]);
-%       axis('square');
-%       title(sprintf('s=[%i %i]',lt1(fignum),lt2(fignum)),'Fontsize',ftz);
-%     end;
-%
-%   See also:  nonsepgabdual
+%   This is a computational subroutine, do not call it directly.
 
 %   AUTHOR : Nicki Holighaus and Peter L. Soendergaard
 %   TESTING: TEST_NONSEPDGT
@@ -112,49 +31,15 @@
 
 % Assert correct input.
 
-if nargin<5
-  error('%s: Too few input parameters.',upper(mfilename));
-end;
-
-definput.keyvals.L=[];
-definput.flags.nsalg={'multiwin','shear'};
-[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);
-
-% Change f to correct shape.
-[f,Ls,W,wasrow,remembershape]=comp_sigreshape_pre(f,upper(mfilename),0);
-
-if isempty(L)
-  L=dgtlength(Ls,a,M,lt);
-else
-  Lcheck=dgtlength(L,a,M,lt);
-  if Lcheck~=L
-    error('%s: Invalid transform size L',upper(mfilename));
-  end;
-end;
-
+L=size(f,1);
+W=size(f,2);
 b=L/M;
 N=L/a;
 
-[g,info] = comp_window(g,a,M,L,lt,'NONSEPDGT');
-
-if (info.isfir)  
-  if info.istight
-    g=g/sqrt(2);
-  end;  
-end;
-
-f=postpad(f,L);
-
-% If the signal is single precision, make the window single precision as
-% well to avoid mismatches.
-if isa(f,'single')
-  g=single(g);
-end;
-
-
 % ----- algorithm starts here, split into sub-lattices ---------------
 
-if flags.do_multiwin
+if (alg==1) || (alg==0 && lt(2)<=2) 
+
     mwin=comp_nonsepwin2multi(g,a,M,lt);
     
     % simple algorithm: split into sublattices
@@ -179,9 +64,7 @@
         c(:,:,w) = c(:,:,w).*repmat(E,M,1);
     end;
     
-end;
-
-if flags.do_shear
+else
 
     b=L/M;
     N=L/a;