function fft_ifft_peak( filePath, fileName )
%--------------------------------------------------------------------------
% FFT to IFFT Peak Tracking Phase Locked Vocoder
%
% adapted from:
% VX_tstretch_real_pv_phaselocked.m   [DAFXbook, 2nd ed., chapter 7]
%
% Cooper Baker - 2014
%--------------------------------------------------------------------------

close all;

% Settings
%--------------------------------------------------------------------------
windowSize    = 1024;
overlap       = 4;
stretchFactor = 2;
window        = hann( windowSize, 'periodic' );
tag           = 'peak';

% Initializations
%--------------------------------------------------------------------------
if any( exist( 'fileName' ) ~= 1 )
    [ fileName, filePath ] = uigetfile( '*.wav', 'Audio File' );
end

analysisHop   = windowSize / overlap;
synthHop      = analysisHop * stretchFactor;
[ input, sr ] = audioread( [ filePath, fileName ] );
halfWinSize   = windowSize / 2;
inputSize     = length( input );
input         = [ zeros( windowSize, 1 ) ; input ; zeros( windowSize - mod( inputSize, analysisHop ), 1 ) ] / max( abs( input ) );
output        = zeros(windowSize + ceil( length( input ) * stretchFactor ), 1 );
omega         = 2 * pi * analysisHop * [ 0 : halfWinSize ]' / windowSize;
phi0          = zeros( halfWinSize + 1, 1 );
psi           = zeros( halfWinSize + 1, 1 );
psi2          = psi;
numPrevPeaks  = 0;
analysisIndex = 0;
analysisMax   = length( input ) - windowSize;
synthIndex    = 0;

% create progress bar dialog box
bar = waitbar( 0, '0%', 'Name', sprintf( '%s: processing %s...', mfilename, fileName ) );

% Processing Loop
%--------------------------------------------------------------------------
while analysisIndex < analysisMax

    % copy and window the input frame
    frame = input( analysisIndex + 1 : analysisIndex + windowSize ) .* window;

    % perform fft and save relevant half of spectrum
    spec = fft( fftshift( frame ) );
    halfSpec = spec( 1 : halfWinSize + 1 );

    % cartesian to polar conversion
    mag   = abs  ( halfSpec );
    phase = angle( halfSpec );

    % find spectral peaks ( local maxima )
    peakLoc  = zeros( halfWinSize + 1, 1 );
    numPeaks = 0;

    for b = 3 : halfWinSize - 1
        if( mag( b ) > mag( b - 1 ) && mag( b ) > mag( b - 2 ) && mag( b ) > mag( b + 1 ) && mag( b ) > mag( b + 2 ) )
            numPeaks = numPeaks + 1;
            peakLoc( numPeaks ) = b;
            b = b + 3;
        end
    end

    % propagate peak phases and compute spectral bin phases
    if( analysisIndex == 0 )
        psi = phase;
    elseif( numPeaks > 0 && numPrevPeaks > 0 )
        prevPeak = 1;

        for p = 1 : numPeaks

            % connect current peak to the previous closest peak
            peak2 = peakLoc( p );
            while( prevPeak < numPrevPeaks && abs( peak2 - prevPeakLoc( prevPeak + 1 ) ) < abs( peak2 - prevPeakLoc( prevPeak ) ) )
                prevPeak = prevPeak + 1;
            end
            peak1 = prevPeakLoc( prevPeak );

            % propagate peak's phase assuming linear frequency
            % variation between connected peak1 and peak2
            % ( avgPeak is a 1-based indexing spectral bin )
            avgPeak           = ( peak1 + peak2 ) * 0.5;
            peakOmega         = 2 * pi * analysisHop * ( avgPeak - 1.0 ) / windowSize;
            peakDeltaPhase    = peakOmega + mod( ( phase( peak2 ) - phi0( peak1 ) - peakOmega ) + pi, -2 * pi ) + pi;
            peakTargetPhase   = mod( ( psi( peak1 ) + peakDeltaPhase * stretchFactor ) + pi, -2 * pi ) + pi;
            peakPhaseRotation = mod( ( peakTargetPhase - phase( peak2 ) ) + pi, -2 * pi ) + pi;

            % rotate phases of all bins around the current peak
            if( numPeaks == 1 )
                bin1 = 1;
                bin2 = halfWinSize + 1;
            elseif( p == 1 )
                bin1 = 1;
                bin2 = halfWinSize + 1;
            elseif( p == numPeaks )
                bin1 = round( ( peakLoc( p - 1 ) + peak2 ) * 0.5 );
                bin2 = halfWinSize + 1;
            else
                bin1 = round( ( peakLoc( p - 1 ) + peak2 ) * 0.5 ) + 1;
                bin2 = round( ( peakLoc( p + 1 ) + peak2 ) * 0.5 );
            end
                psi2( bin1 : bin2 ) = mod( ( phase( bin1 : bin2 ) + peakPhaseRotation ) + pi, -2 * pi ) + pi;
        end
            psi = psi2;
    else
        deltaPhase = omega + mod( ( phase - phi0 - omega ) + pi, -2 * pi ) + pi;
        psi = mod( ( psi + deltaPhase * stretchFactor ) + pi, -2 * pi ) + pi;
    end

    % polar to cartesian conversion
    spec = ( mag .* exp( 1i * psi ) );

    % reconstruct entire spectrum
    spec = [ spec( 1 : halfWinSize + 1 ) ; conj( spec( halfWinSize : -1 : 2 ) ) ];

    % perform an ifft on the spectrum
    frame = ifft( spec );

    % discard imaginary data
    frame = real( frame );

    % shift zero frequency component to center of spectrum
    frame = fftshift( frame );

    % apply the window
    frame = frame .* window;

    % overlap-add the synthesized frame into the output buffer
    output( synthIndex + 1 : synthIndex + windowSize ) = output( synthIndex + 1 : synthIndex + windowSize ) + frame;

    % store values for next frame
    phi0          = phase;
    prevPeakLoc   = peakLoc;
    numPrevPeaks  = numPeaks;

    % increment analysis and synthesis indices
    analysisIndex = analysisIndex + analysisHop;
    synthIndex    = synthIndex + synthHop;

    % update the progress bar
    progress = ( analysisIndex / analysisMax );
    waitbar( progress, bar, sprintf( '%2.3f%%', progress * 100 ) );
end

% close progress bar dialog box
close( bar );

% Output
%--------------------------------------------------------------------------

% crop output buffer
output = output( windowSize + 1 : length( output ) );

% normalize output buffer
normCoeff = 1 / max( abs( output ) );
output    = output * normCoeff;

% plot output
timevector = linspace( 0, length( output )/sr, length( output ) );
plot( timevector, output );

% create comment string
commentString = sprintf( 'WinSize: %.0f,\nOverlap: %.0f,\nStretch: %.2f,\nWindow: %s', windowSize, overlap, stretchFactor, 'Hann' );

% write audio file to disk
[ a, fileName, b ] = fileparts( fileName );
outFile = sprintf( '%s.%s.wav', fileName, tag );
audiowrite( outFile, output, sr, 'BitsPerSample', 32, 'Artist', 'NormCoeff', 'Title', num2str( normCoeff ), 'Comment', commentString );

% EOF