clear
"[filename, pathname] = uigetfile('*.*', 'EMG file');"
cd(pathname)
a=load(filename);
pause

"%assume we get 1 sec of data, if sampling frequency =1000 then window size=1000"
windowsize=1024;
ii = 1:windowsize ;

"[r,c]=size(a);"

"x=a(:,1);"

% use the hanning window
windowfr=hanning(windowsize);

% perform the short time Fourier Transform. Also caluclate the median frequency
k=0;
iii=k*windowsize+ii;
while (k*windowsize+11)<r
    k=k+1;
    iii=k*windowsize+ii;
if iii<r
    signal=x(iii);
"    Pxx=psd(signal, windowsize, 1024, windowfr);"
    plength=length(Pxx);
    for p=1:plength
"        Pxx_all(k,p)=Pxx(p);"
    end
    median_freq(k)=median(Pxx);
end
end


%Show graphs for the TF distribution
figure(1)
%contour(Pxx_all)
mesh(Pxx_all)
"%view(-22 , 90);"
axis tight
xlabel('frequency (Hz)');
ylabel('Time (sec)');
%zlabel('Amplitude (V)');

% Show graph - the development of median frequency
figure(2)
plot(median_freq)
xlabel('Time (sec)');
ylabel('Normalised Frequency');