From 4f222663ab1ed007a6eb9776a949406e74576e29 Mon Sep 17 00:00:00 2001 From: "Peter L. Soendergaard" Date: Tue, 2 Oct 2012 18:54:31 +0200 Subject: [PATCH] help cleanup before workshop --- comp/comp_window.m | 8 ++-- demos/demo_audiocompression.m | 44 ++++++++++---------- demos/demo_auditoryfilterbank.m | 43 +------------------- filterbank/filterbank.m | 15 ++++++- filterbank/ufilterbank.m | 10 +++++ fourier/pfilt.m | 11 +++++ frames/frsyniter.m | 71 ++++++++++++++++++++++++++++----- gabor/dgt.m | 27 +++++++++++-- gabor/dwilt.m | 3 +- gabor/gabdual.m | 48 +++++++++++++++------- gabor/gabdualnorm.m | 4 +- gabor/gabframebounds.m | 25 ++++++------ gabor/gabtight.m | 48 +++++++++++++++------- gabor/plotwmdct.m | 2 +- gabor/wmdct.m | 21 ++++++++-- testing/test_nonsepdgt.m | 26 ++++++------ 16 files changed, 263 insertions(+), 143 deletions(-) diff --git a/comp/comp_window.m b/comp/comp_window.m index 31aa1327..dfa4be0f 100644 --- 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,7 +94,7 @@ 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'} @@ -102,7 +102,7 @@ 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 index 05fab9a1..38500281 100644 --- 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 index 76b159d1..297c4ae9 100644 --- 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,14 +26,6 @@ % 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 @@ -41,7 +33,7 @@ % 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 index 070c6390..9e59eecc 100644 --- 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 index daad03c0..67917960 100644 --- 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 index d4ab8b8e..89e11f5c 100644 --- 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 index f1465b02..0fa64363 100644 --- 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,19 +29,29 @@ % % '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; @@ -49,25 +59,64 @@ 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 diff --git a/gabor/dgt.m b/gabor/dgt.m index 091c0222..a1df82bf 100644 --- a/gabor/dgt.m +++ b/gabor/dgt.m @@ -2,9 +2,8 @@ %DGT Discrete Gabor transform % Usage: c=dgt(f,g,a,M); % c=dgt(f,g,a,M,L); -% c=dgt(f,g,a,M,L,lt); -% [c,Ls]=dgt(f,g,a,M); -% [c,Ls]=dgt(f,g,a,M,L); +% c=dgt(f,g,a,M,'lt',lt); +% [c,Ls]=dgt(...); % % Input parameters: % f : Input data. @@ -96,6 +95,28 @@ % 'timeinv' Compute a DGT using a time-invariant phase. This % convention is typically used in filter bank algorithms. % +% Examples: +% --------- +% +% In the following example we create a Hermite function, which is a +% complex-valued function with a circular spectrogram, and visualize +% the coefficients using both `imagesc` and |plotdgt|_::: +% +% a=10; +% M=40; +% L=a*M; +% h=pherm(L,4); % 4th order hermite function. +% c=dgt(h,'gauss',a,M); +% +% % Simple plot: The squared modulus of the coefficients on linear scale +% figure(1); +% imagesc(abs(c).^2); +% +% % Better plot: zero-frequency is displayed in the middle, and the +% % coefficients are show on a logarithmic scale. +% figure(2); +% plotdgt(c,a,'dynrange',50); +% % See also: idgt, gabwin, dwilt, gabdual, phaselock % % Demos: demo_dgt diff --git a/gabor/dwilt.m b/gabor/dwilt.m index 4ffaac88..9121fb2d 100644 --- a/gabor/dwilt.m +++ b/gabor/dwilt.m @@ -2,8 +2,7 @@ %DWILT Discrete Wilson transform % Usage: c=dwilt(f,g,M); % c=dwilt(f,g,M,L); -% [c,Ls]=dwilt(f,g,M); -% [c,Ls]=dwilt(f,g,M,L); +% [c,Ls]=dwilt(...); % % Input parameters: % f : Input data diff --git a/gabor/gabdual.m b/gabor/gabdual.m index dc93ea31..a783e153 100644 --- a/gabor/gabdual.m +++ b/gabor/gabdual.m @@ -2,7 +2,7 @@ %GABDUAL Canonical dual window of Gabor frame % Usage: gd=gabdual(g,a,M); % gd=gabdual(g,a,M,L); -% gd=gabdual(g,a,M,L,lt); +% gd=gabdual(g,a,M,'lt',lt); % % Input parameters: % g : Gabor window. @@ -28,13 +28,33 @@ % of length *L*. Unless the dual window is a FIR window, the dual window % will have length *L*. % -% `gabdual(g,a,M,[],lt)` or `gabdual(g,a,M,L,lt)` does the same for a -% non-separable lattice specified by *lt*. Please see the help of -% |matrix2latticetype|_ for a precise description of the parameter *lt*. +% `gabdual(g,a,M,'lt',lt)` does the same for a non-separable lattice +% specified by *lt*. Please see the help of |matrix2latticetype|_ for a +% precise description of the parameter *lt*. % % If $a>M$ then the dual window of the Gabor Riesz sequence with window % *g* and parameters *a* and *M* will be calculated. % +% Examples: +% --------- +% +% The following example shows the canonical dual window of the Gaussian +% window::: +% +% a=20; +% M=30; +% L=300; +% g=pgauss(L,a*M/L); +% gd=gabdual(g,a,M); +% +% % Simple plot in the time-domain +% figure(1); +% plot(gd); +% +% % Frequency domain +% figure(2); +% magresp(gd,'dynrange',100); +% % See also: gabtight, gabwin, fir2long, dgt % AUTHOR : Peter Soendergaard. @@ -50,7 +70,7 @@ definput.keyvals.L=[]; definput.keyvals.lt=[0 1]; definput.keyvals.nsalg=0; -[flags,kv,L,lt]=ltfatarghelper({'L','lt'},definput,varargin); +[flags,kv,L]=ltfatarghelper({'L'},definput,varargin); %% ------ step 2: Verify a, M and L if isempty(L) @@ -63,13 +83,13 @@ end; % ----- step 2b : Verify a, M and get L from the window length ---------- - L=dgtlength(Ls,a,M,lt); + L=dgtlength(Ls,a,M,kv.lt); else % ----- step 2a : Verify a, M and get L - Luser=dgtlength(L,a,M,lt); + Luser=dgtlength(L,a,M,kv.lt); if Luser~=L error(['%s: Incorrect transform length L=%i specified. Next valid length ' ... 'is L=%i. See the help of DGTLENGTH for the requirements.'],... @@ -80,7 +100,7 @@ %% ----- step 3 : Determine the window -[g,info]=gabwin(g,a,M,L,lt,'callfun',upper(mfilename)); +[g,info]=gabwin(g,a,M,L,kv.lt,'callfun',upper(mfilename)); if LM$ then an orthonormal window of the Gabor Riesz sequence with % window *g* and parameters *a* and *M* will be calculated. % +% Examples: +% --------- +% +% The following example shows the canonical tight window of the Gaussian +% window. This is calculated by default by |gabtight|_ if no window is +% specified::: +% +% a=20; +% M=30; +% L=300; +% gt=gabtight(a,M,L); +% +% % Simple plot in the time-domain +% figure(1); +% plot(gt); +% +% % Frequency domain +% figure(2); +% magresp(gt,'dynrange',100); +% % See also: gabdual, gabwin, fir2long, dgt % AUTHOR : Peter Soendergaard. @@ -75,7 +95,7 @@ definput.keyvals.L=[]; definput.keyvals.lt=[0 1]; definput.keyvals.nsalg=0; -[flags,kv,L,lt]=ltfatarghelper({'L','lt'},definput,varargin); +[flags,kv,L]=ltfatarghelper({'L'},definput,varargin); %% ------ step 2: Verify a, M and L @@ -89,13 +109,13 @@ end; % ----- step 2b : Verify a, M and get L from the window length ---------- - L=dgtlength(Ls,a,M,lt); + L=dgtlength(Ls,a,M,kv.lt); else % ----- step 2a : Verify a, M and get L - Luser=dgtlength(L,a,M,lt); + Luser=dgtlength(L,a,M,kv.lt); if Luser~=L error(['%s: Incorrect transform length L=%i specified. Next valid length ' ... 'is L=%i. See the help of DGTLENGTH for the requirements.'],... @@ -107,7 +127,7 @@ %% ----- step 3 : Determine the window -[g,info]=gabwin(g,a,M,L,lt,'callfun',upper(mfilename)); +[g,info]=gabwin(g,a,M,L,kv.lt,'callfun',upper(mfilename)); if L