boblsturm@0: function y = istft(spec,parameter, varargin)
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: % Name: istft
boblsturm@0: % Date: 03-2014
boblsturm@0: % Programmer: Jonathan Driedger
boblsturm@0: % http://www.audiolabs-erlangen.de/resources/MIR/TSMtoolbox/
boblsturm@0: %
boblsturm@0: % Computing the 'inverse' of the stft according to the paper "Signal 
boblsturm@0: % Estimation from Modified Short-Time Fourier Transform" by Griffin and
boblsturm@0: % Lim.
boblsturm@0: %
boblsturm@0: % Input:   spec             a complex spectrogram generated by stft.
boblsturm@0: %          parameter.
boblsturm@0: %             synHop        hop size of the synthesis window.
boblsturm@0: %             win           the synthesis window.
boblsturm@0: %             zeroPad       number of zeros that were padded to the
boblsturm@0: %                           window to increase the fft size and therefore
boblsturm@0: %                           the frequency resolution.
boblsturm@0: %             numOfIter     number of iterations the algorithm should
boblsturm@0: %                           perform to adapt the phase.
boblsturm@0: %             origSigLen    original length of the audio signal such that
boblsturm@0: %                           the output can be trimmed accordingly.
boblsturm@0: %             restoreEnergy when windowing the synthesis frames, there is a
boblsturm@0: %                           potential for some energy loss. This option
boblsturm@0: %                           will rescale every windowed synthesis frame to
boblsturm@0: %                           compensate for this energy leakage.
boblsturm@0: %             fftShift      in case the stft was computed with an fftShift,
boblsturm@0: %                           setting this parameter to 1 will compensate for
boblsturm@0: %                           that by applying the same shifting operation 
boblsturm@0: %                           again to each frame after the ifft.
boblsturm@0: %
boblsturm@0: % Output:  y                the time-domain signal.
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: % Reference: 
boblsturm@0: % If you use the 'TSM toolbox' please refer to:
boblsturm@0: % [DM14] Jonathan Driedger, Meinard Mueller
boblsturm@0: %        TSM Toolbox: MATLAB Implementations of Time-Scale Modification 
boblsturm@0: %        Algorithms
boblsturm@0: %        Proceedings of the 17th International Conference on Digital Audio  
boblsturm@0: %        Effects, Erlangen, Germany, 2014.
boblsturm@0: %
boblsturm@0: % License:
boblsturm@0: % This file is part of 'TSM toolbox'.
boblsturm@0: % 
boblsturm@0: % 'TSM toolbox' is free software: you can redistribute it and/or modify it
boblsturm@0: % under the terms of the GNU General Public License as published by the
boblsturm@0: % the Free Software Foundation, either version 3 of the License, or (at
boblsturm@0: % your option) any later version.
boblsturm@0: % 
boblsturm@0: % 'TSM toolbox' is distributed in the hope that it will be useful, but
boblsturm@0: % WITHOUT ANY WARRANTY; without even the implied warranty of
boblsturm@0: % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
boblsturm@0: % Public License for more details.
boblsturm@0: % 
boblsturm@0: % You should have received a copy of the GNU General Public License along
boblsturm@0: % with 'TSM toolbox'. If not, see http://www.gnu.org/licenses/.
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: 
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: % check parameters
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: if nargin<2
boblsturm@0:     parameter=[];
boblsturm@0: end
boblsturm@0: 
boblsturm@0: if ~isfield(parameter,'synHop')
boblsturm@0:     parameter.synHop = 2048;
boblsturm@0: end
boblsturm@0: if ~isfield(parameter,'win')
boblsturm@0:     parameter.win = win(4096,2); % hann window
boblsturm@0: end
boblsturm@0: if ~isfield(parameter,'zeroPad')
boblsturm@0:     parameter.zeroPad = 0;
boblsturm@0: end
boblsturm@0: if ~isfield(parameter,'numOfIter')
boblsturm@0:     parameter.numOfIter = 1;
boblsturm@0: end
boblsturm@0: if ~isfield(parameter,'origSigLen')
boblsturm@0:     parameter.origSigLen = -1; % no trimming
boblsturm@0: end
boblsturm@0: if ~isfield(parameter,'restoreEnergy')
boblsturm@0:     parameter.restoreEnergy = 0;
boblsturm@0: end
boblsturm@0: if ~isfield(parameter,'fftShift')
boblsturm@0:     parameter.fftShift = 0;
boblsturm@0: end
boblsturm@0: 
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: % some pre calculations
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: numOfFrames = size(spec,2);
boblsturm@0: numOfIter = parameter.numOfIter;
boblsturm@0: 
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: % audio calculation
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: % first iteration
boblsturm@0: Yi = spec;
boblsturm@0: yi = LSEE_MSTFT(Yi,parameter);
boblsturm@0: 
boblsturm@0: % remaining iterations
boblsturm@0: parStft.win = parameter.win;
boblsturm@0: parStft.numOfFrames = numOfFrames;
boblsturm@0: parStft.zeroPad = parameter.zeroPad;
boblsturm@0: parStft.anaHop = parameter.synHop;
boblsturm@0: for j = 2 : numOfIter
boblsturm@0:     Yi = abs(spec) .* exp(1j*angle(stft(yi,parStft)));
boblsturm@0:     yi = LSEE_MSTFT(Yi,parameter);
boblsturm@0: end
boblsturm@0: y = yi';
boblsturm@0: y = y./max(y);
boblsturm@0: 
boblsturm@0: % if the original Length of the signal is known, also remove the zero
boblsturm@0: % padding at the end
boblsturm@0: if parameter.origSigLen > 0
boblsturm@0:     y = y(1:parameter.origSigLen);
boblsturm@0: end
boblsturm@0: 
boblsturm@0: end
boblsturm@0: 
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: % the Griffin Lim procedure
boblsturm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
boblsturm@0: function x = LSEE_MSTFT(X,parameter)
boblsturm@0: % some pre calculations
boblsturm@0: w = parameter.win;
boblsturm@0: w = w(:);
boblsturm@0: zp = parameter.zeroPad;
boblsturm@0: w = [zeros(floor(zp/2),1);w;zeros(floor(zp/2),1)];
boblsturm@0: winLen = length(w);
boblsturm@0: winLenHalf = round(winLen/2);
boblsturm@0: synHop = parameter.synHop;
boblsturm@0: numOfFrames = size(X,2);
boblsturm@0: winPos = (0:numOfFrames-1) * synHop + 1;
boblsturm@0: signalLength = winPos(end) + winLen - 1;
boblsturm@0: 
boblsturm@0: x = zeros(signalLength,1); % resynthesized signal
boblsturm@0: ow = zeros(signalLength,1); % sum of the overlapping windows
boblsturm@0: for i = 1 : numOfFrames
boblsturm@0:     currSpec = X(:,i);
boblsturm@0:     
boblsturm@0:     % add the conjugate complex symmetric upper half of the spectrum
boblsturm@0:     Xi = [currSpec;conj(currSpec(end-1:-1:2))];
boblsturm@0:     xi = real(ifft(Xi));
boblsturm@0:     if parameter.fftShift == 1
boblsturm@0:         xi = fftshift(xi);
boblsturm@0:     end
boblsturm@0:     xiw = xi .* w;
boblsturm@0:     
boblsturm@0:     if parameter.restoreEnergy == 1
boblsturm@0:         xiEnergy = sum(abs(xi));
boblsturm@0:         xiwEnergy = sum(abs(xiw));
boblsturm@0:         xiw = xiw * (xiEnergy/(xiwEnergy+eps));
boblsturm@0:     end
boblsturm@0:     
boblsturm@0:     x(winPos(i):winPos(i)+winLen-1) = ...
boblsturm@0:         x(winPos(i):winPos(i)+winLen-1) + xiw;
boblsturm@0:     
boblsturm@0:     ow(winPos(i):winPos(i)+winLen-1) = ...
boblsturm@0:         ow(winPos(i):winPos(i)+winLen-1) + w.^2;
boblsturm@0: end
boblsturm@0: ow(ow<10^-3) = 1; % avoid potential division by zero
boblsturm@0: x = x ./ ow;
boblsturm@0: 
boblsturm@0: % knowing the zeropads that were added in the stft computation, we can
boblsturm@0: % remove them again now. But since we do not know exactly how many
boblsturm@0: % zeros were padded at the end of the signal, it is only safe to remove
boblsturm@0: % winLenHalf zeros.
boblsturm@0: x = x(winLenHalf+1:end-winLenHalf);
boblsturm@0: end