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

%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');