function time_freq_plot(x,fs,plot_option,wind_option)
% TIME_FREQ_PLOT(x,fs,plot_option,wind_option) plots both the time signal,x, 
%   and the FFT of x (X) for a given sampling frequency of fs.  Remaining 
%   input arguments are options and must all be specified or not.  If x is 
%   complex, X is plotted from 0 to fs else X is plotted from 0 to fs/2.
%
%   For normalized frequency set fs = 1.  
%   For normalized radian frequency set fs = 2Pi.
%   Any other fs is assumed to be in units of Hertz.
%
%   If plot_option = 1, the magnitude response in dB is given. (default)
%   If plot_option = 2, the magnitude response is given.
%   If plot_option = 3, the magnitude squared response is given.
%
%   If wind_option = 0 a rectangular window is used.
%   If wind_option = 1 a Hanning window applied.  (default)

%------------------
% Check valid input
%------------------
if ((nargin ~= 2) & (nargin ~= 4))
   error('Must have 2 or 4 input arguments.');
end;

if (nargin == 2)
   plot_option = 1;
   wind_option = 1;
else  % check valid input options
   if ((plot_option ~= 1) & (plot_option ~= 2) & (plot_option ~= 3))
      error('plot_option must be 1, 2, or 3.');
   end;
   if ((wind_option ~= 0) & (wind_option ~= 1))
      error('wind_option must be 0 or 1.');
   end;
end;

%-----------
% Initialize
%-----------
if (any(abs(imag(x)) > 10*eps))
   complex = 1;
else
   x = real(x);
   complex = 0;
end;

%----------
% Time plot
%----------
subplot(211),

plotline = 0; % 0 solid line, 1 -- , 2 -. , 3 :

if (complex)
   disp('Warning: x is complex.  Abs(x) will be displayed.');
   if (plotline == 0)
      plot(abs(x)),
   elseif (plotline == 1)
      plot(abs(x),'--'),
   elseif (plotline == 2)
      plot(abs(x),'-.'),
   else
      plot(abs(x),':'),
   end;
   ylabel('|X(w)| (dB)'),
else
   if (plotline == 0)
      plot(x),
   elseif (plotline == 1)
      plot(x,'--'),
   elseif (plotline == 2)
      plot(x,'-.'),
   else
      plot(x,':'),
   end;
   ylabel('|X(w)| (dB)'),
   axis([1,length(x),min(x),max(x)]);
end;

title('x (Time Domain)'),
xlabel('k'),
ylabel('x(k)'),
grid;

%-----------
% Freq. Plot
%-----------
axis;  % set back to auto-scaling
N = length(x);

% Optional Hanning Window
if (wind_option)
   hann_wind = hanning(N);
   hann_gain = 3/8;
   x = x .* hann_wind / sqrt(hann_gain);  % scale by window gain
end;

% Transform to frequency domain
X = fft(x);
if (complex)
   f = fs*(0:N-1)/N;  % frequency scale
   X = X(1:N)/sqrt(N/2);  % plot 0 to fs
else
   f = fs*(0:N/2-1)/N;  % frequency scale
   X = X(1:N/2)/sqrt(N/2);  % plot 0 to fs/2
end;

subplot(212),

if (plot_option == 1)                            % dB magnitude plot
   if (plotline == 0)
      plot(f,10*log10(X .* conj(X))),
   elseif (plotline == 1)
      plot(f,10*log10(X .* conj(X)),'--'),
   elseif (plotline == 2)
      plot(f,10*log10(X .* conj(X)),'-.'),
   else
      plot(f,10*log10(X .* conj(X)),':'),
   end;
   ylabel('|X(w)| (dB)'),
elseif (plot_option == 2)                        % magnitude plot
   if (plotline == 0)
      plot(f,abs(X)),
   elseif (plotline == 1)
      plot(f,abs(X),'--'),
   elseif (plotline == 2)
      plot(f,abs(X),'-.'),
   else
      plot(f,abs(X),':'),
   end;
   ylabel('|X(w)|'),
else                                             % magnitude squared plot
   if (plotline == 0)
      plot(f,X .* conj(X)),
   elseif (plotline == 1)
      plot(f,X .* conj(X),'--'),
   elseif (plotline == 2)
      plot(f,X .* conj(X),'-.'),
   else
      plot(f,X .* conj(X),':'),
   end;
   ylabel('|X(w)|^2'),
end;

if (fs == 1)
   xlabel('Normalized Frequency (f/fs)'),
elseif (fs == 2*pi)
   xlabel('Normalized Radian Frequency (theta)'),
else
   xlabel('Frequency (Hz)'),
end;

if (wind_option)
   title('X (Frequency Domain) (Hanning Windowed)'),
else
   title('X (Frequency Domain)'),
end;

grid;

