function [output] = audiomod(input,rate,fs)
    %audio: audio signal to be modified (e.g., speech or music).
    %output: output audio
    %rate: Playback speed, rate > 0
    %fs: sample rate

    %Audio signal prep
    [audio,fs] = audioread(input);

    %Experiments with Block Length led us to a BlockLength of 1000
    %and a hop length of 200 for the best outputs. This results in
    %an overlap of about 80%.
    iBlockLength = 1000;
    iHopLength=iBlockLength/5;
    numBlocks = floor(length(audio)/iBlockLength);
    fftlen= 2*iBlockLength;

    %Keep track of the Phase of the audio signal
    phase=zeros(1,iBlockLength);

    %Set Output to appropriate length with some additional padding
    output=zeros(1,floor(length(audio)/rate+iBlockLength));

    %Iterate through the signal at appropriate rate of hopLength * rate
    for i = 1:round(iHopLength*rate):length(audio)-(iBlockLength+iHopLength)+1
        %Two windows of samples
        win1=audio(i:i+iBlockLength-1);
        win2=audio(i+iHopLength:i+iBlockLength-1+iHopLength);

        %FFT of each window, element wise multiplication to hamming windows
        s1=fft(win1.*hamming(iBlockLength))';
        s2=fft(win2.*hamming(iBlockLength))';

        %Modify the phase in accordance to the iteration
        phase=mod((phase+angle(s1./s2)),2*pi);

        %Synthesis
        ablock=real(ifft(abs(s1).*exp(1j*phase)))';
        n=round(i/rate);
        output(n:n+iBlockLength-1)= output(n:n+iBlockLength-1) + (ablock.*hamming(iBlockLength))';
    end
    soundsc(output,fs);
end