function [tot, waxis] = wigner (x0,title)

% based on the wig2 rotuine from the HOSA (higher Order Spectrum Analysis) toolbox
% Adapted by AJL 5/2000

showcontour = 1;

%WIG2	 Computes the Wigner Distribution (second-order)
%	[wx,waxis] = wig2 (x, nfft,flag)
%	x     - time series, must be a vector
%	nfft  - FFT length to use; default is the power of 2 just larger
%	        than twice the length of x.
%	flag  - By default, if signal 'x' is real valued, its analytic form
%	        is used to compute the WD; this helps supress cross terms
%	        around D.C.;  if flag is 0, the analytic form is not used.
%	wx    - the Wigner distribution
%	        rows correspond to time, columns to frequencies
%	        time increases with row number, frequencies with col number
%	waxis - the frequency axis associated with the WD
%       It is recommended that the analytic form of the signal be used.

%  Copyright (c) 1991-1999 by United Signals & Systems, Inc. and The Mathworks, Inc. All Rights Reserved.
%       $Revision: 1.7 $
%  A. Swami   January 20, 1995

%     RESTRICTED RIGHTS LEGEND
% Use, duplication, or disclosure by the Government is subject to
% restrictions as set forth in subparagraph (c) (1) (ii) of the
% Rights in Technical Data and Computer Software clause of DFARS
% 252.227-7013.
% Manufacturer: United Signals & Systems, Inc., P.O. Box 2374,
% Culver City, California 90231.
%
%  This material may be reproduced by or for the U.S. Government pursuant
%  to the copyright license under the clause at DFARS 252.227-7013.

% --------------------- parameter checks --------------------------
[p, m] = size(x0); % If x0 is not a 1xm vector we assume we want to average the spectrum results 
if (min(m,p) ~= 1)
   disp 'This might take awhile'
end

if (m == 1) x0 = x0'; [p,m]=size(x0); end; % assume they made a mistake passing the vector transposed

lfft = 2^nextpow2(m);                      % minimum FFT length

if (m < lfft)
   disp 'WARNING:FFT will be slower since length not power of two'
end

tot = zeros(m,2*m); % deliberately swapped for later
waxis = [-m:m-1] / (2*2*m);
taxis = 1:m;

% Now loop over the number of lines
for (numlines = 1:p)
 x = zeros(m*2,1);   x(1:m) = x0(numlines,:);          cx = conj(x);
 wx = zeros(m*2,m);
 L1 = m-1;

 % --------- compute r(tau,t) = cx(t-tau/2) * x(t+tau/2) -------------

 for n=0:L1
   indm = max(-n,-L1+n) : min(n,L1-n);
   indy = indm + (indm < 0) * 2 * m ;   % output indices y(m;n)
   y = zeros(2*m,1);
   y(indy + 1) = x(n+indm + 1) .* cx(n-indm + 1);
   wx(:,n+1)   = y;
 end

 % ----------- WD(f,t) = FT (tau-->f) r(tau,t) ---------------------
 wx = fft(wx);
 wx = real(wx.');                  % force it to be real
 wx = abs(wx); %AJL
 
 % ----------- display the WD ---------------------------------------
 %  note the frequency scaling by 2M
 %  WD(f,t) =  X(2f,t),  where X(f,t) = FT ( r(tau,t) )


 wx = wx(:,[m+1:2*m,1:m]) ;

%remove these for batch (non-display work)
%   if (showcontour == 1)
%     clf;
%     contour(waxis,taxis, abs(wx),8), grid on 
%     ylabel('time in samples')
%     xlabel('frequency')
%     set(gcf,'Name','Individual Wigner Ville Distribution')
%     drawnow;
%    else
%     myimagedisplay(wx, 'Wigner Ville Distribution of slice');
%     xlabel('time/space');
%     ylabel('frequency/spatial frequency')
%    end;

  % total the results
  tot = tot + wx;
end;


     myimagedisplay(tot,title);
     save( title, 'waxis', 'taxis', 'tot'); % save away for later work to show contour
%     xlabel('time/space')
%     ylabel('frequency/spatial frequency')
%     figure;
%     contour(waxis,taxis, abs(tot),8), grid on 
%     ylabel('time in samples')
%     xlabel('frequency')
%     set(gcf,'Name','Contour Aggregate Wigner Ville Distribution')
%     drawnow;