%% Signal Definition

[file, path] = uigetfile('*.wav', 'Select a wave file');
nameoffile = fullfile(path,file);

[y,fs] = audioread(nameoffile);
y = y(:,1);
y = y(fix(10*fs):fix(20*fs));

%% Frames and Windows

frame_size=20;
frame_shift=10;
window_length=(frame_size/1000)*fs;
sample_shift=(frame_shift/1000)*fs;
no_of_windows = (floor((length(y))/sample_shift)-ceil(window_length/sample_shift));
f0=zeros(1,no_of_windows);
maxamp = 0;
flag_y= envelop_hilbert(y);

%% Autocorrelation

for i=1: no_of_windows
    
    sig = y(i*sample_shift:(i*sample_shift)+window_length);
    thres = flag_y(i*sample_shift:(i*sample_shift)+window_length);
    activity_flag = max(abs(thres));    
    
    if activity_flag == 1                     % discarding low amplitude freq
        cor = xcorr(sig);
        [pks,locs] = findpeaks(cor);
        [midval, midloc] = max(pks);
        [maxval,maxloc] = max(pks(midloc+1:end));
        
        if isequal(maxval,maxloc)
            f0(i) = 0;
        else
            period = (abs(locs(midloc)-locs(maxloc+midloc)));
            f0(i) = fs / period;
        end
    else
        f0(i) = 0;
    end
    
end

v = f0 .* 4.0526e-2;

%% plots
t=1/fs:1/fs:(length(y)/fs);

[peaks,locs] = findpeaks(f0,'MinpeakHeight',0.5*max(f0),'MinPeakDistance',50);ylim([0 1200]);xlim([0 2000]);

plot(1:no_of_windows,f0*4.0526e-2,locs,peaks*4.0526e-2,'o');title('Pitch of each window');ylabel('Velocity (cm/s)');xlabel('Windows');
ylim([0 100]);
g = zoom;
g.enable = 'on';
set(gcf, 'Position', get(0,'Screensize'));

max_velocity = max(peaks)*4.0526e-2;
Velocity = mean(peaks)*4.0526e-2;

fprintf('\nMax Systolic Peak Velocity     = %.2f cm/s\n',max_velocity);
fprintf('\nAverage Systolic Peak Velocity = %.2f cm/s\n\n',Velocity);