function [f,relres,iter]=isgram(s,g,a,varargin)
%ISGRAM Spectrogram inversion
% Usage: f=isgram(c,g,a);
% f=isgram(c,g,a,Ls);
% [f,relres,iter]=isgram(...);
%
% Input parameters:
% c : Array of coefficients.
% g : Window function.
% a : Length of time shift.
% Ls : length of signal.
% Output parameters:
% f : Signal.
% relres : Vector of residuals.
% iter : Number of iterations done.
%
% ISGRAM(s,g,a) attempts to invert a spectrogram computed by :
%
% s = abs(dgt(f,g,a,M)).^2;
%
% using an iterative method.
%
% ISGRAM(c,g,a,Ls) does as above but cuts or extends f to length Ls.
%
% If the phase of the spectrogram is known, it is much better to use
% idgt.
%
% [f,relres,iter]=ISGRAM(...) additionally return the residuals in a
% vector relres and the number of iteration steps iter.
%
% Generally, if the spectrogram has not been modified, the iterative
% algorithm will converge slowly to the correct result. If the
% spectrogram has been modified, the algorithm is not guaranteed to
% converge at all.
%
% ISGRAM takes the following parameters at the end of the line of input
% arguments:
%
% 'lt',lt Specify the lattice type. See the help on MATRIX2LATTICETYPE.
%
% 'zero' Choose a starting phase of zero. This is the default
%
% 'rand' Choose a random starting phase.
%
% 'int' Construct a starting phase by integration. Only works
% for Gaussian windows.
%
% 'griflim' Use the Griffin-Lim iterative method, this is the
% default.
%
% 'bfgs' Use the limited-memory Broyden Fletcher Goldfarb
% Shanno (BFGS) method.
%
% 'tol',t Stop if relative residual error is less than the specified tolerance.
%
% 'maxit',n Do at most n iterations.
%
% 'print' Display the progress.
%
% 'quiet' Don't print anything, this is the default.
%
% 'printstep',p If 'print' is specified, then print every p'th
% iteration. Default value is p=10;
%
% The BFGS method makes use of the minFunc software. To use the BFGS method,
% please install the minFunc software from:
% http://www.cs.ubc.ca/~schmidtm/Software/minFunc.html.
%
% Examples:
% ---------
%
% To reconstruct the phase of 'greasy', use the following:
%
% % Setup the problem and the coefficients
% f=greasy;
% g='gauss';
% a=20; M=200;
% c=dgt(f,g,a,M);
% s=abs(c).^2;
% theta=angle(c);
%
% % Reconstruct and get spectrogram and angle
% r=isgram(s,g,a);
% c_r=dgt(r,g,a,M);
% s_r=abs(c_r).^2;
% theta_r=angle(c_r);
%
% % Compute the angular difference
% d1=abs(theta-theta_r);
% d2=2*pi-d1;
% anglediff=min(d1,d2);
%
% % Plot the difference in spectrogram and phase
% figure(1);
% plotdgt(s./s_r,a,16000,'clim',[-10,10]);
% colormap([bone;flipud(bone)])
% title('Relative difference in spectrogram');
%
% figure(2);
% plotdgt(anglediff,a,16000,'lin');
% colormap(bone);
% title('Difference in angle');
%
% See also: isgramreal, frsynabs, dgt
%
% References:
% R. Decorsière and P. L. Søndergaard. Modulation filtering using an
% optimization approach to spectrogram reconstruction. In Proceedings of
% the Forum Acousticum, 2011.
%
% D. Griffin and J. Lim. Signal estimation from modified short-time
% Fourier transform. IEEE Trans. Acoust. Speech Signal Process.,
% 32(2):236--243, 1984.
%
% D. Liu and J. Nocedal. On the limited memory BFGS method for large
% scale optimization. Mathematical programming, 45(1):503--528, 1989.
%
%
% Url: http://ltfat.github.io/doc/gabor/isgram.html
% Copyright (C) 2005-2023 Peter L. Soendergaard <peter@sonderport.dk> and others.
% This file is part of LTFAT version 2.6.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% AUTHOR : Remi Decorsiere and Peter L. Søndergaard.
% REFERENCE: OK
% 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.keyvals.lt=[0 1];
definput.keyvals.tol=1e-6;
definput.keyvals.maxit=100;
definput.keyvals.printstep=10;
definput.flags.method={'griflim','bfgs'};
definput.flags.print={'quiet','print'};
definput.flags.startphase={'zero','rand','int'};
[flags,kv,Ls]=ltfatarghelper({'Ls','tol','maxit'},definput,varargin);
M=size(s,1);
N=size(s,2);
W=size(s,3);
if ~isnumeric(a) || ~isscalar(a)
error('%s: "a" must be a scalar',upper(mfilename));
end;
if rem(a,1)~=0
error('%s: "a" must be an integer',upper(mfilename));
end;
L=N*a;
Ltest=dgtlength(L,a,M,kv.lt);
if Ltest~=L
error(['%s: Incorrect size of coefficient array or "a" parameter. See ' ...
'the help of DGTLENGTH for the requirements.'], ...
upper(mfilename))
end;
sqrt_s=sqrt(s);
if flags.do_zero
% Start with a phase of zero.
c=sqrt_s;
end;
if flags.do_rand
c=sqrt_s.*exp(2*pi*i*rand(M,N));
end;
if flags.do_int
if kv.lt(2)>1
error(['%s: The integration initilization is not implemented for ' ...
'non-sep lattices.'],upper(mfilename));
end;
c=constructphase(s,g,a);
end;
% gabwin is called after constructphase above, because constructphase
% needs the window as a string/cell
g=gabwin(g,a,M,L,kv.lt,'callfun',upper(mfilename));
gd = gabdual(g,a,M,L);
% For normalization purposes
norm_s=norm(s,'fro');
relres=zeros(kv.maxit,1);
if flags.do_griflim
for iter=1:kv.maxit
%c = comp_dgt_proj(c,g,gd,a,M,L);
if kv.lt(2)==1
f=comp_idgt(c,gd,a,kv.lt,0,0);
c=comp_dgt(f,g,a,M,kv.lt,0,0,0);
else
f=comp_idgt(c,gd,a,kv.lt,0,0);
c=comp_dgt(f,g,a,M,kv.lt,0,0,0);
end;
relres(iter)=norm(abs(c).^2-s,'fro')/norm_s;
c=sqrt_s.*exp(i*angle(c));
if flags.do_print
if mod(iter,kv.printstep)==0
fprintf('ISGRAM: Iteration %i, residual = %f\n',iter,relres(iter));
end;
end;
if relres(iter)<kv.tol
relres=relres(1:iter);
break;
end;
end;
f=comp_idgt(c,gd,a,kv.lt,0,0);
end;
if flags.do_bfgs
if exist('minFunc')~=2
error(['To use the BFGS method in ISGRAMREAL, please install the minFunc ' ...
'software from http://www.cs.ubc.ca/~schmidtm/Software/minFunc.html.']);
end;
% Setting up the options for minFunc
opts = struct;
opts.display = kv.printstep;
opts.maxiter = kv.maxit;
% Don't limit the number of function evaluations, just the number of
% time-steps.
opts.MaxFunEvals = 1e9;
opts.usemex = 0;
f0 = comp_idgt(c,gd,a,kv.lt,0,0);
[f,fval,exitflag,output]=minFunc(@objfun,f0,opts,g,a,M,s,kv.lt);
% First entry of output.trace.fval is the objective function
% evaluated on the initial input. Skip it to be consistent.
relres = output.trace.fval(2:end)/norm_s;
iter = output.iterations;
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,0,[0; W]);
% Subfunction to compute the objective function for the BFGS method.
function [f,df]=objfun(x,g,a,M,s,lt);
L=size(s,2)*a;
c=comp_dgt(x,g,a,M,lt,0,0,0);
inner=abs(c).^2-s;
f=norm(inner,'fro')^2;
df=4*real(conj(comp_idgt(inner.*c,g,a,lt,0,0)));