function signal = harmonize(x, scale1, scale2, amp, Fs, winLen)
% x - input signal audio
% Fs - Sample Rate of the incoming signal
% winLen - the length of the Hanning Window and the FFT size 
% hopFactor - factor by which the Hanning Window moves. hopFactor +
% [Window Overlap] = 1
% This function creates simulated harmonics around the fundamental
% frequency of the signal.  A better implementation of this would be to
% find the strongest 3 peaks and shift those frequencies.  This could also
% lead to some formant effects
% written by John Arroyo, ja2124@columbia.edu

if nargin < 6
    winLen = 1024; % a window length of 256 gives a nice synthesized speech sound
end

hopFactor = 1/4; %1/4 for best reconstruction

% a 75% window overlap needed for a smooth reconstruction.
% the integer is equal to 1-x%=hopFactor, 50% - 75% overlap is optimal which
% gives the vocoder the ability to handle fractional time scaling shifts 
% (fractions of the hamming window)
hopSize = winLen*hopFactor;

% encode the signal using stft (using a hann window)
X = encode(x, winLen, hopSize);

% Interpolate as is done with time shifting and create the robot effect
% Xint = interpolate(X, time, hopSize);
Xrows = size(X,1);
N = 2*(Xrows-1);

% Calculate the expected phase advance
expPhaseAdv = zeros(N/2+1,1);
expPhaseAdv(2:(1 + N/2),1) = (2*pi*hopSize)./(N./(1:(N/2))');

% Phase accumulator
% Preset to phase of first frame
% in case of 1:1 time scaling
phaseX = angle(X(:,1));

winIndex = 1; %begin the window index with 1
sizeX = size(X,2)-2;

%number of samples to grab around the harmonic note
windowScale = .033; %percentage of the window (bin) that will be used to create the width of the harmonic
numSamples = ceil(size(X,1)*windowScale); % should be odd since size of X is odd, is 17 in the case of a 1024 hann window

for i = 1:sizeX-1
    % magnitude of the input signal
    Xmag = abs(X(:,i));
    
    % get the largest frequency and find out the shifted frequency value
    sampleIndex = floor(numSamples/2);
    peak = max(Xmag);
    myIndex = find(Xmag == peak);
    freq = (myIndex/N)*(Fs); %frequency of the peak

    shiftIndex1 = floor(N*freq*scale1/Fs); % find the index of the first harmonic
    shiftIndex2 = floor(N*freq*scale2/Fs); % find the index of the second harmonic
    
    % get the frequency components around the fundamental frequency 
    % (or strongest frequency)
    Xmag2 = zeros(N/2+1,1);
    Xmag2(shiftIndex1-sampleIndex:shiftIndex1+sampleIndex) = amp*X(myIndex-sampleIndex:myIndex+sampleIndex);
    Xmag2(shiftIndex2-sampleIndex:shiftIndex2+sampleIndex) = amp*X(myIndex-sampleIndex:myIndex+sampleIndex); %compensate for overlap
  
    phaseX2 = zeros(N/2+1,1);
    phaseX2(shiftIndex1-sampleIndex:shiftIndex1+sampleIndex) = phaseX(myIndex-sampleIndex:myIndex+sampleIndex);
    phaseX2(shiftIndex2-sampleIndex:shiftIndex2+sampleIndex) = phaseX(myIndex-sampleIndex:myIndex+sampleIndex);
    
    % Generate the the current column of Xint 
    % add the scaled harmonics to the original signal
    Xint(:,winIndex) = Xmag .* exp(j*phaseX) + Xmag2(1:N/2+1) .* exp(j*phaseX2(1:N/2+1)); %x = mag*e^jw
    winIndex = winIndex+1;
    
    % Generate the phase advance
    % the phase advance is the difference of the original phases minus the
    % expected phase advance
    phaseAdv = angle(X(:,i+2)) - angle(X(:,i+1)) - expPhaseAdv;
    
    % Generate Phase Accumulator
    phaseX = phaseX + expPhaseAdv + phaseAdv;
end

% turn back into a waveform, inverse window and fft
signal = decode(Xint, winLen, hopSize)';