function [stftfcn,t,t2,f1]=stftfcn(data,fsamp,tsub_goal);
tsamp=1/fsamp;
points_per_subwindow=floor(tsub_goal/tsamp);
tsubwindow=points_per_subwindow*tsamp;
subwindows=floor(length(data)/points_per_subwindow);
n1=[0:points_per_subwindow-1];
t=[0:length(data)-1]*tsamp;  % Whole time axis
t1=n1*tsamp; % Time axis within sample
n2=[0:subwindows];
t2=n2*points_per_subwindow*tsamp; % Time axis of samples
f1=n1*(1/tsubwindow);
fmin=1/tsubwindow
sdata=reshape(data(1:points_per_subwindow*subwindows),...
              points_per_subwindow,subwindows);
hdata=diag(hamming(length(f1)))*sdata;
stftfcn=fft(hdata);
