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

nameoffile = fullfile(path,file);

[signal,fs] = audioread(nameoffile);

signal = signal(:,1);

velocity = 1; %Average Peak Systolic Velocity, '1' : yes | '0' : no

type = 'velocity';      %'freq' or 'velocity' on y-axis

time = 5;               %plot for specific time period (in seconds)

t=1:fix(time*fs);

%signal=signal(t);

frame = 21.3;             %size for spectrogram (ms)

s_length = length(signal);

sampleTime = ( 1:s_length )/fs;

frameSize = fix(frame*0.001*fs);

figure()

if ( strcmp(type, 'freq') )

        [B,f,T] = spectrogram(signal,hanning(frameSize),round(frameSize/2),frameSize*2,fs);

        B = 20*log10(abs(B));

        imagesc(T,f,B);axis xy;colorbar;ylim([0 5000]);caxis([-25 25]);%datacursormode on;

        t=colormap(gray);

        colormap(1-t);

        h = zoom;

        h.motion = 'horizontal';

        h.enable = 'on';

        xlabel('Time (s)');ylabel('Frequency (Hz)'); 

end

if ( strcmp(type, 'velocity') )

        [B,f,T] = spectrogram(signal,hann(frameSize),round(frameSize/2),frameSize*2,fs);

        B = 20*log10(abs(B));

        v = f*0.040526;

        imagesc(T,v,B);axis xy;ylim([0 100]);colorbar;caxis([-25 25]);

        t=colormap(gray);

        colormap(1-t);

        g = zoom;

        g.motion = 'horizontal';

        g.enable = 'on';

        xlabel('Time (s)');ylabel('Velocity (cm/s)');

end

set(gcf, 'Position', get(0,'Screensize'));

if velocity

[a,b] = find(B>-5);

vel = v(a);

Time = T(b);

[p,loc] = findpeaks(vel,'MinPeakHeight',0.7*max(vel),'MinPeakDistance',500);

figure()

plot(Time,vel,Time(loc),p,'o');xlabel('Time (s)');ylabel('velocity (cm/s)');ylim([0 0.040526*4000]);

set(gcf, 'Position', get(0,'Screensize'));

max_velocity = max(p);

Velocity = mean(p);

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

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

end