function [grain_out, grain_phi, grain_psi] = autotune(grain, ...
                                                grain_old_phi,  ...
                                                grain_old_psi, ...
                                                comparisonPitches, FS)

% PHASE VOCODER ===========================================================
s_win    = length(grain);
hs_win   = s_win/2;
Ra       = s_win/4;

w1       = hanning(s_win, 'periodic');   % Analysis Window
w2       = w1;                           % Synthesis Window

omega    = 2*pi*Ra*[0:hs_win-1]'/s_win;   % 2piR k/N

psi      = grain_old_psi;
phi0     = grain_old_phi;

% GRAIN PROCESSING ========================================================
grain = grain .* w1;
f     = fft(grain);
f     = f(1:hs_win);
r     = abs(f);
phi   = angle(f);

% Calculation of delta_phi as phase difference between current segment
% and the previous segment during analysis stage.
% This phase difference takes part in calculating the result phase in the
% synthesis stage, since the synthesized phase of the current segment is
% the previous synthesized phase + delta_phi*realPitchShift.
% The obtained effect is a shift in the frequency content of the signal.
delta_phi = omega + princarg(phi-phi0-omega);

% TONE FREQ. EXTRACTION ================================================
% First individuation of fund. frequency based on magnitude
[r_max, rMaxIndex] = max(r);

estimFrequency = 1/(2*pi)*delta_phi(rMaxIndex)*FS/(Ra);

% Calculate difference between 'correct' tones and measured frequency
estimDifference = abs(comparisonPitches - estimFrequency);

% Get argmin to obtain the closest frequency index value
[min_val, estimNearIndex] = min(estimDifference);

% FREQ. SHIFTING ========================================================
% Calculate the pitch shift ratio for the current tone adjustement
shiftRatio = comparisonPitches(estimNearIndex) / estimFrequency;

if isinf(shiftRatio) || isnan(shiftRatio)
    shiftRatio = 1;
end

psi   = princarg(psi+delta_phi*shiftRatio);

% RESYNTHESIS ===========================================================
ft    = (r.* exp(1i*psi));                     % Rebuilding FFT spectrum
ft = [ ft(1:hs_win) ; conj(ft(hs_win:-1:1))];

% Resynthesizing frequency shifted grain
grain_out = real(ifft(ft)).*w2;
grain_phi = phi;
grain_psi = psi;
