function rs = ispect_cos(instruct, outfile)
% iSPECT_COS	Interactive version of the namesake script 
%     which demonstrates the effects of window length on 
%     spectral resolution.
%
%     iSPECT_COS(WINLEN) plots the DFTs of a series of signals 
%     which consist of 2 frequency components with varying 
%     amounts of separation.  WINLEN specifies the length of
%     a rectangular window with which the signal is multiplied.
%     The length of the window determines the frequency resolution 
%     in the resulting DFT.
%
% Author:      J. Greenberg
% Modifications: D. Tran
% Class : HST.582J
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set values of constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global siglen fixedfftlen; % Specify variables as global
siglen = 1000;		   % length of signal
fixedfftlen = 1024;	   % fixed length used to calculate DTFT

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STAGE 1	Initialize the return string
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rs = char('');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STAGE 2	Set directory for storage of graphic files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change working dir to temp dir on ACS machine
% NOTE: We won't make use of the Matweb variable MLDIR
tempdir = '/var/local/httpd/htdocs/matweb/tmp';

% Uncomment below when running script from a stand-alone test file
%tempdir = '/afs/athena/course/6/6.555/matweb/tmp';
cd(tempdir);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STAGE 3	Get HTML form input variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Get unique identifier (to form file name)
mlid = getfield(instruct, 'mlid');

% Show time signal? Field non-existent if box unchecked
try
   showtime = getfield(instruct, 'showtime');
catch
   showtime = ' no';
end

% Show previous plot?
try
   compareplot = getfield(instruct, 'compareplot');
catch
   compareplot = ' no';
end

% Get current plot's color
ccolor = lower(getfield(instruct, 'ccolor'));
ccolorshape = [ccolor 'o'];
% Get saved plot's color
scolor = lower(getfield(instruct, 'scolor'));
scolorshape = [scolor 'd'];
% Get standard parameters
[A0,A1,f0,f1,winlen,wintype,fftlen] = getParams(instruct);

% Get saved parameters
if compareplot == 'yes'
   [sA0,sA1,sf0,sf1,swinlen,swintype,sfftlen] = getSaved(instruct);

   % Create plot title for debugging purpose
%   stitle = sprintf('Saved values: A0=%.2f, A1=%.2f, f0=%.2f, f1=%.2f, winlen=%d, wintype=%s, fftlen=%d', sA0, sA1, sf0, sf1, swinlen, swintype, sfftlen); 

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STAGE 4	Perform MATLAB computations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find FFT for standard parameters
[w,v,W,Vf,Vk] = compFFT(A0,A1,f0,f1,winlen,wintype,fftlen);

% Find FFT for saved parameters
if compareplot == 'yes'
   [sw,sv,sW,sVf,sVk] = compFFT(sA0,sA1,sf0,sf1,swinlen,swintype,sfftlen); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate graphics
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Specify time axes
upperxlim = floor(winlen/2);    % upper limit of window
if (mod(winlen,2) == 1)		% if window length is odd
   lowerxlim = -upperxlim;	% lower limit of window
else
   lowerxlim = -(upperxlim-1);
end
% axis for window
axt1 = [lowerxlim-1 upperxlim+1 min(w) max(w)+0.2*abs(max(w))];
% axis for windowed signal 
axt2 = [axt1(1) axt1(2) min(v)-0.2*abs(min(v)) max(v)+0.2*abs(max(v))];
% NOTE: We don't plot the saved time plots for now 
% as it would make the plot unnecessarily cluttered

% Specify frequency axes
fixedincr = 1./fixedfftlen;               % fixed frequency increment
fixedf = [-0.5:fixedincr:0.5-fixedincr]'; % fixed frequency range
incr = 1./fftlen;		% frequency increment
f = [-0.5:incr:0.5-incr]';	% frequency scale
if (compareplot == 'yes')
   sincr = 1./sfftlen;	                % saved frequency increment
   sf = [-0.5:sincr:0.5-sincr]';	% saved frequency scale
   axf1 = [-0.5 0.5 0 1.1*max([ max(W) max(sW) ])];   % axis for FFT of window
   axf2 = [-0.5 0.5 0 1.1*max([ max(Vf) max(sVf) ])]; % axis for FFT of signal
else
   axf1 = [-0.5 0.5 0 1.1*max(W)];
   axf2 = [-0.5 0.5 0 1.1*max(Vf)];
end

% NOTE: do not use 'hold' alone as it prints output which upsets MATWEB
% hold on;    % Turn hold on for all subplots

% Perform housekeeping
wscleanup('ml*spect_cos.jpeg', 1); % Delete jpegs older than 1 hour

% Create "invisible" figure 
fhandle = figure('visible','off');
clf;

% Plot
if showtime == 'yes'

   % Adjust the size and position of the JPEG image
   set(fhandle, 'PaperPosition', [0 0 8 6])

   % Plot w[n]
   subplot(3,2,1)
   coolplot([lowerxlim:upperxlim]',w(1:winlen),ccolorshape)
   axis(axt1)
   ylabel('w[n]')

   % Plot |W(f)|
   subplot(3,2,2)
   if (compareplot == 'yes')
      slegend = plot(fixedf,sW,scolor);
      hold on; % placement may cause problem with display of plot border
   end
   curlegend = plot(fixedf,W,ccolor);
   axis(axf1)
   set(gca,'YAxisLocation','right')
   ylabel('|W(f)|')

   % Plot v[n]
   subplot(3,2,3)
   coolplot([lowerxlim:upperxlim]',v(1:winlen),ccolorshape)
   axis(axt2) 
   ylabel('v[n]')

   % Plot |V(f)|
   subplot(3,2,4)
   if (compareplot == 'yes')
      plot(fixedf,sVf,scolor)
      hold on;
   end
   plot(fixedf,Vf,ccolor)
   axis(axf2)
   set(gca,'YAxisLocation','right')
   ylabel('|V(f)|')

   % Empty plot
   subplot(3,2,5)  % empty plot
   ah = gca;       % get handle so we can put a legend on this axis
   axis(axt2);	   % but show labels
   set(ah, 'TickDir', 'out')       % tick direction is outward
   %set(ah, 'XTickLabel', '')      % kill x-axis
   set(ah, 'YTickLabel', '', 'YTick', [])      % kill y-axis
   set(ah, 'Box', 'on')			       % turn on box

   % Plot |V[k]|
   subplot(3,2,6)
   if (compareplot == 'yes')
      coolplot(sf,sVk,scolorshape);
      hold on;
   end
   coolplot(f,Vk,ccolorshape);
   axis(axf2);
   set(gca,'YAxisLocation','right');
   xlabel('frequency'), ylabel('|V[k]|')

   subplotmx(3,2);   % compress the subplots: 3 rows, 2 columns

else

   % Adjust the size of the JPEG image
   set(fhandle, 'PaperPosition', [0 0 8 6])

   % Plot |W(f)|
   subplot(3,1,1)
   ah = gca;   % Get handle so we can put a legend on this axis
   if (compareplot == 'yes')
      slegend = plot(fixedf,sW,scolor);
      hold on; % placement may cause problem with display of plot border
   end
   curlegend = plot(fixedf,W,ccolor);
   axis(axf1);
   ylabel('|W(f)|')

   % Plot |V(f)|
   subplot(3,1,2)
   if (compareplot == 'yes')
     plot(fixedf,sVf,scolor);
     hold on;
   end
   plot(fixedf,Vf,ccolor);
   axis(axf2);
   ylabel('|V(f)|')

   % Plot |V[k]|
   subplot(3,1,3)
   if (compareplot == 'yes')
      coolplot(sf,sVk,scolorshape);
      hold on;
   end
   coolplot(f,Vk,ccolorshape);
   axis(axf2);
   xlabel('frequency'), ylabel('|V[k]|')

   subplotmx(3,1);   % compress subplots: 3 rows, 1 column

end

% Create legend.  Must be _after_ compressing subplots as 
% plots with legend are not affected by subplotmx   
if (compareplot == 'yes')
  lh = legend(ah,[curlegend, slegend], 'Current', 'Saved',0);
  % The following if statement should not be needed, it was
  % necessary to add it after the version of Matlab running on the
  % server changed.  I believe this is a Matlab bug, and will
  % hopefully resolve in the future.  Right now, this prevents
  % Matlab from automatically placing the legend box to avoid
  % conflicts.
  if (strcmp(showtime, 'yes'))
    set(lh,'Position',[0.25 0.22 0.13 0.06])  
  else
    set(lh,'Position',[0.75 0.8 0.13 0.06])
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STAGE 5	Put output variables in an output structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Render jpeg and write to file.
drawnow;
s.GraphFileName = sprintf('%sspect_cos.jpeg', mlid);
wsprintjpeg(fhandle, s.GraphFileName);
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STAGE 6	Call HTMLREP with output struct and output HTML
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Put name of graphic file into HTML template file.
templatefile = '/afs/athena/course/6/6.555/matweb/samp/ispect2_cos.html';
if (nargin == 1)
   rs = htmlrep(s, templatefile);
elseif (nargin == 2)
   rs = htmlrep(s, templatefile, outfile);
end



%-----------------------------------------------------------------
% Subfunctions
%-----------------------------------------------------------------

% getParams(instruct) gets the values of standard variables from 
% the input structure
function [A0,A1,f0,f1,winlen,wintype,fftlen] ...
= getParams(instruct)

   % Get amplitudes of cosines
   A0 = str2double(getfield(instruct,'A0'));
   A1 = str2double(getfield(instruct,'A1'));

   % Get frequencies of cosines
   f0 = str2double(getfield(instruct,'f0'));
   f1 = str2double(getfield(instruct,'f1'));

   % Get window length
   winlen = str2double(getfield(instruct, 'winlen'));
   % Get window type
   wintype = lower(getfield(instruct, 'wintype'));

   % Get FFT length
   fftlen = str2double(getfield(instruct, 'fftlen'));

% getSaved(instruct) gets the values of saved variables from 
% the input structure.  Saved variable names begin with 's'.
function [sA0,sA1,sf0,sf1,swinlen,swintype,sfftlen] ...
= getSaved(instruct)

   % Get amplitudes of cosines
   sA0 = str2double(getfield(instruct,'sA0'));
   sA1 = str2double(getfield(instruct,'sA1'));

   % Get frequencies of cosines
   sf0 = str2double(getfield(instruct,'sf0'));
   sf1 = str2double(getfield(instruct,'sf1'));

   % Get window length
   swinlen = str2double(getfield(instruct, 'swinlen'));
   % Get window type
   swintype = lower(getfield(instruct, 'swintype'));

   % Get FFT length
   sfftlen = str2double(getfield(instruct, 'sfftlen'));

% compFFT(A0,A1,f0,f1,winlen,wintype,fftlen) find the FFT of the
% window and the windowed signal
function [w,v,W,Vf,Vk] = compFFT(A0,A1,f0,f1,winlen,wintype,fftlen) 

   % Use global parameters
   global siglen fixedfftlen;

   % Determine window 
   switch wintype
      case 'bartlett'
         w = bartlett(winlen);
      case 'blackman'
         w = blackman(winlen);
      case 'hamming'
         w = hamming(winlen);
      case 'hanning'
         w = hanning(winlen);
      case 'boxcar'
         w = boxcar(winlen);
   end
   w = [w; zeros(siglen-winlen,1)];  % zero-padding

   % Create signal
   n = [1:siglen]';
   cos0 = cos(2*pi*f0*n);
   cos1 = cos(2*pi*f1*n);
   x = A0*cos0 + A1*cos1;
   % Window signal
   v = w .* x;				  
   % Calculate FFT
   W = fftshift(abs(fft(w,fixedfftlen)));	  % DTFT of window
   Vf = fftshift(abs(fft(v,fixedfftlen))); % DTFT of windowed signal
   if fftlen == fixedfftlen
      Vk = Vf;				  
   else
      Vk = fftshift(abs(fft(v,fftlen)));     % DFT of windowed signal
   end
