%
% QFT computes the amplitude spectrum of a data 
% sequence. The user is prompted for an input signal (a vector)
% and whether drift removal, tapering, and mean removal should be done.
% A vector QFT_out of the amplitude spectrum is created along with
% a vector QFT_f of frequencies in units of cy/sample. 
%
% Enter the name of the (workspace) signal vector in the 'SIGNAL NAME? box, click on
% the check boxes as appropriate and then click OK. The cycle can be repeated any 
% number of times with different options. Click END to finish.
%
% Plots will be drawn of the raw data, the modified data, and the amplitude
% spectrum. The amplitude apectrum can be redrawn with plot(QFT_f,QFT_out).

%***********************************************************************************
%
%                              Jim Merriam
%                       University of Saskatchewan
%                          jim.merriam@usask.ca
%                              (Feb, 1993)  
%                         (last Oct 1998)
%***********************************************************************************
clf

disp(' ')
disp(' QFT is a FT m THAT ALLOWS INTERACTIVE CONTROL OF ')
disp(' TAPERING, DRIFT AND DC REMOVAL ')
disp(' OPERATE IN THE CURRENT FIGURE WINDOW')
disp(' ')

clear QFT_out QFT_f x_mod
clear x_mod MeanOut N_samples TaperEnds taper DriftOut
clear OK c f h1 h2 h3 h4 he s M_taper qft_in

% uicontrol for the name of the signal vector
he=uicontrol('style','edit','units','normalized','position',[.7 .325 .2 .05], ... 
   'string','F?','Callback','qft_in=eval(get(he,''string''));');

% uicontrol to keep cycling through options list
h5=uicontrol('style','pushbutton','units','normalized','position',[.7 .01 .1 .05], ...
'string','END','Callback','keepon=''N'';');
keepon='Y';
while keepon=='Y' 
 
% uicontrol to remove mean 
MeanOut='N';
h1=uicontrol('style','checkbox','units','normalized','position',[.7 .1 .2 .05], ... 
'string','REMOVE MEAN?','Callback','MeanOut=''Y'';');

% uicontrol for drift linear removal
DriftOut='N';
h2=uicontrol('style','checkbox','units','normalized','position',[.7 .175 .2 .05],...
'string','REMOVE DRIFT?','Callback','DriftOut=''Y'';');

% uicontrol for ends tapering 
TaperEnds='N';
h3=uicontrol('style','checkbox','units','normalized','position',[.7 .25 .2 .05], ...
'string','TAPER ENDS?','Callback','TaperEnds=''Y'';');

% uicontrol to pause while options are selected
OK='PAUSEaBIT';
h4=uicontrol('style','pushbutton','units','normalized','position',[.8 .01 .1 .05], ...
'string','OK','Callback','OK=''continue '';');
while OK=='PAUSEaBIT' & keepon=='Y'
   pause(0.1)
end

if(keepon=='N')
  break		% changed from 'return', in order to close the calculations properly - I.M. Oct 06 2014
end % return to keyboard 
  
N_samples=length(qft_in);
 if(N_samples~=floor(N_samples/2)*2); % if N is odd
 N_samples=N_samples-1;               % truncate
 qft_in(1)=[];                        % by deleting first sample
 disp(' TRUNCATING BY DELETING FIRST SAMPLE')  
 end
 
 [n,m]=size(qft_in);
 if(m>n)
 qft_in=qft_in';% make sure qft_in is a column vector
 end
 
 s=[1:N_samples]'; % set up sample vector

subplot(221)
plot(s,qft_in,'.','MarkerSize',6)
title('RAW DATA'),xlabel('SAMPLE NUMBER')

x_mod=qft_in; % copy input vector to working vector

if DriftOut == 'Y'
  c=polyfit(s,x_mod,1); % fit and remove a linear drift
  f=polyval(c,s);
  x_mod=x_mod-f;
  end

  if TaperEnds == 'Y'
  M_taper=floor(.1*N_samples); % taper first and last ten percent to their mean 

 % taper=hann(2*M_taper);    % with cosine taper
 taper=0.5*cos(-pi*[-M_taper:1:M_taper]'/M_taper) + 0.5; % Hanning taper formula corrected - I.M. Oct 6 2014
  mid_point=(x_mod(1)+x_mod(N_samples))/2; % avg of first and last 
  x_mod=x_mod-mid_point; % subtract avg of first and last% close the figure window - I.M. Oct 06 2014
  x_mod(1:M_taper)=x_mod(1:M_taper).*taper(1:M_taper);
  x_mod(N_samples-M_taper+1:N_samples)= ... 
     x_mod(N_samples-M_taper+1:N_samples).*taper(M_taper+1:2*M_taper);
  x_mod=x_mod+mid_point; % add avg of first and last back in
  else;
  end;
  
  if MeanOut == 'Y'
  x_mod=x_mod-mean(x_mod); % remove the mean
  end;

subplot(222)% close the figure window - I.M. Oct 06 2014
plot(s,x_mod,'.','MarkerSize',6) % plot modified data vector
title('MODIFIED')
xlabel('SAMPLE NUMBER')

P=fft(x_mod);
P=2*sqrt(P.*conj(P))/N_samples;
QFT_out(1:N_samples/2)=P(1:N_samples/2); % amplitude spectrum
QFT_out=QFT_out';
QFT_f=(0:N_samples/2-1)'/(N_samples); % frequencies cy/sample
clear P;

subplot(223)
plot(QFT_f,QFT_out,'.','MarkerSize',6) % plot amplitude spectrum
title('AMPLITUDE SPECTRUM')
xlabel('FREQUENCY (cy/sample)'),ylabel('AMPLITUDE');


end % end while keepon =='Y'

close;		% close the figure window - I.M. Oct 06 2014

clear x_mod MeanOut N_samples TaperEnds taper DriftOut
clear OK c f h1 h2 h3 h4 he s M_taper qft_in