function signal = chorus(x, amp, Fs, winLen)
% x - input signal audio
% Fs - Sample Rate of the incoming signal
% amp - amplitude of the chorus components
% 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 a simulated 'robot' voice that attempts to create
% synthesized speech effect.  So that the input sounds as if it was 
% generated from a computer.
% written by John Arroyo, ja2124@columbia.edu

if nargin < 4
    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 = freq - 1; % find the index of the first chorus component
    shiftIndex2 = freq + 1; % find the index of the second chorus component

    scaleFactor = freq - 1; % 2 hz %(4 steps) 1.18 (3 steps, ie A4 to C5)
    shiftIndex1 = floor(N*scaleFactor/Sr); % find the index of the first harmonic
  
    scaleFactor2 = freq + 1; %+ 2 Hz 
    shiftIndex2 = floor(N*scaleFactor2/Sr); % find the index of the second harmonic
    
    % get the frequency components around the fundamental frequency 
    % (or strongest frequency)
    Xmaga = zeros(N/2+1,1);
    Xmagb = zeros(N/2+1,1);
    Xmaga(shiftIndex1-sampleIndex:shiftIndex1+sampleIndex) = amp*X(myIndex-sampleIndex:myIndex+sampleIndex);
    Xmagb(shiftIndex2-sampleIndex:shiftIndex2+sampleIndex) = amp*X(myIndex-sampleIndex:myIndex+sampleIndex); %compensate for overlap
    Xmag2 = Xmaga + Xmagb;
    
    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)';