RTPGHI - Real-Time Phase Gradient Integration

function [c,newphase,tgrad,fgrad]=rtpghi(s,gamma,a,M,varargin)
%RTPGHI Real-Time Phase Gradient Integration
%   Usage:  c=rtpghi(s,gamma,a,M);
%           [c,newphase,tgrad,fgrad] = rtpghi(...);
%
%   Input parameters:
%         s        : Initial coefficients.
%         gamma    : Window width factor.
%         a        : Hop factor.
%         M        : Number of channels.
%   Output parameters:
%         c        : Coefficients with the constructed phase.
%         newphase : Just the (unwrapped) phase.
%         tgrad    : Relative time phase derivative.
%         fgrad    : Relative frequency phase derivative.
% 
%   RTPGHI(s,gamma,a,M) creates complex DGTREAL coefficients from their
%   absolute values s using the Real-Time Phase Gradient Heap Integration
%   algorithm. s must have been obtained as:
%
%       c = dgtreal(f,g,a,M,'timeinv');
%       s = abs(c);
%
%   and the algorithm attempts to recover c. Parameter gamma is window 
%   g specific and it can be computed using PGHI_FINDGAMMA.
%
%   This function works entirely simiral to PGHI except it is using
%   the real-time version of the algorithm. Please see help of PGHI 
%   (resp. CONSTRUCTPHASEREAL from LTFAT) for more details.
%
%   Algorithm version:
%
%       'normal'    1 frame delay version of the algorithm RTPGHI(1)
%
%       'causal'    No delay version of the algorithm RTPGHI(0)
%
%   See also: dgtreal, idgtreal, pghi
%
%   References:
%     Z. Průša, P. Balazs, and P. L. Søndergaard. A Noniterative Method for
%     Reconstruction of Phase from STFT Magnitude. IEEE/ACM Trans. on Audio,
%     Speech, and Lang. Process., 25(5), May 2017.
%     
%     Z. Průša and P. L. Søndergaard. Real-Time Spectrogram Inversion Using
%     Phase Gradient Heap Integration. In Proc. Int. Conf. Digital Audio
%     Effects (DAFx-16), Sep 2016.
%     
%
%
%   Url: http://ltfat.github.io/phaseret/doc/gabor/rtpghi.html

% Copyright (C) 2016 Zdenek Prusa <zdenek.prusa@gmail.com>.
% This file is part of PHASERET version 0.2.1
%
% 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/>.

%   AUTHORS: Zdenek Prusa
%
thismfilename = upper(mfilename);
complainif_notposint(a,'a',thismfilename);
complainif_notposint(M,'M',thismfilename);

if ~isscalar(gamma)
    error('%s: gamma must be a scalar.',upper(mfilename));
end

definput.keyvals.tol=1e-6;
definput.flags.phase={'timeinv','freqinv'};
definput.flags.variant={'normal','causal'};
[flags,kv]=ltfatarghelper({},definput,varargin);
tol = kv.tol;
[M2,N,W] = size(s);

if W>1
    error('%s: *s* must not be 3 dimensional.',thismfilename);
end

M2true = floor(M/2) + 1;

if M2true ~= M2
    error('%s: Mismatch between *M* and the size of *s*.',thismfilename);
end

abss = abs(s);

[tgrad, fgrad, logs] = comp_pghiphasegrad( abss, gamma, a, M, flags.do_timeinv, flags.do_causal);

newphase = zeros(M2,N);
c = zeros(M2,N);

for n=1:N
    idx = mod( n-1-1:n-1, N ) + 1;
    nprev = mod( n-2, N ) + 1;
    newphase(:,n) = comp_rtpghiupdate(logs(:,idx),tgrad(:,idx),fgrad(:,n),newphase(:,nprev),tol,M);
    c(:,n)=abss(:,n).*exp(1i*newphase(:,n));
end