% Download zip in current directory?
if (~exist('Granular','dir'))
    downloadZipInCurrentDirectory=1;
    if downloadZipInCurrentDirectory
        unzip('http://users.cs.cf.ac.uk/Dave.Marshall/Multimedia/Lecture_Examples/Synthesis/Granular.zip');
        addpath('Granular/');
    end
    
    cd('Granular/');
end

[filename,path] = uigetfile({'*.wav;*.waV;','Wav Files'; ...
            '*.*', 'All files (*.*)'}, ...
            'Select a sound file');
		if isequal(filename,0) || isequal(path,0)
			disp('No file selected!')
			return;
		end
filenamepath = [path filename];

[x, fs] = audioread(filenamepath);

pitchmarks = PitchMarker(x);

%
% Show pitchmarks detected in the audio (small section shown)

% Create an array to mark corrosponding pitch mark position in audio

pitch_bins = zeros(1,length(x));
max_x = max(x);
for i = 1:length(pitchmarks)
pitch_bins(pitchmarks(i)) = max_x;
end


plot(x(1000:6000))
hold on
stem(pitch_bins(1000:6000),'g')

%
% Set up Parameters and do PSOLA:

alpha =   0.9; % Time stretching factor alpha 
beta =   3; % Pitch stretching factor beta 
gamma=1; % Formant shift factor -- more suitable for human voice manipulation

x_psola = psolaF(x,pitchmarks,alpha,beta,gamma);

%
% Plot and play the resultant pitch shifted stretched audio:

plot(x_psola)

Play_x_psola = 1;

if Play_x_psola
  figure
  
  sound(x_psola,fs);
end
 
% audiowrite('x1_PSOLA', x_out,fs);
	

%
%
% Supporting Functions
% PitchMarker.m:
%

function [ pitch ] = PitchMarker(inaudio)
% PitchMarker.m (From Digital Audio Effects)
% Finds all the pitch marks in the input file and returns the
% markings in a matrix


% initial settings

mark=[1:length(inaudio)]*0;
last_pos=1;
place=1;

blocksize=300;

i=1;

while last_pos+floor(blocksize*1.7) < length(inaudio)
    % grabs the next block to examine
    temp=inaudio(last_pos+50:last_pos+floor(blocksize*1.7));
    % finds the high point in the block
    [mag,place]=max(temp);
    % checks to see if there is really a signal in this block
    if mag < 0.01
        place=length(temp);
        mode = 0;
        mark(place+last_pos+50)=1;
        pitch(i)=place+last_pos+50;
    else
        mode = 1;
    end
    % checks to see if there is a pitch mark before the current pitch mark
    while mode == 1
        % finds largest point in block from beginning to current pitch mark
        [mag2,place2]=max(temp(1:place-50));
        % checks to see if high mark is suffincent size to be a pitch mark
        if mag2 > 0.90*mag
            mag=mag2;
            place=place2;
        else
            mode = 0;
            mark(place+last_pos+50)=1;
            pitch(i)=place+last_pos+50;
        end
    end
    % starts the next block to be examind 50 samples after this block
    blocksize=place+50;
    % makes sure next blocksize is of sufficent size
    if blocksize < 150
        blocksize=150;
    end
    last_pos=place+last_pos+50;
    i=i+1;
end

end

% PsolaF.m

function out=psolaF(in,m,alpha,beta,gamma)
%     in    input signal
%     m     pitch marks
%     alpha time stretching factor
%     beta  pitch shifting factor
%     gamma newFormantFreq/oldFormantFreq

P = diff(m);   %compute pitch periods

if m(1)<=P(1), %remove first pitch mark
    m=m(2:length(m));
    P=P(2:length(P));
end

if m(length(m))+P(length(P))>length(in) %remove last pitch mark
    m=m(1:length(m)-1);
else
    P=[P P(length(P))];
end

Lout=ceil(length(in)*alpha);
out=zeros(1,Lout); %output signal


tk = P(1)+1;                  %output pitch mark
while round(tk)<Lout
    [minimum i]=min(abs(alpha*m-tk) );    % find analysis segment
    pit=P(i);pitStr=floor(pit/gamma);
    gr=in(m(i)-pit:m(i)+pit).*hanning(2*pit+1);
    gr=interp1(-pit:1:pit,gr,-pitStr*gamma:gamma:pit);% stretch segm.
    iniGr=round(tk)-pitStr;endGr=round(tk)+pitStr;
    if endGr>Lout, break; end
    out(iniGr:endGr)=out(iniGr:endGr)+gr; % overlap new segment
    tk=tk+pit/beta;
end % end of while
end