function [T,F,Z]=Wigner(D,fm);

% y: input signal
% D: the data matrix: D(:,1)=t; D(:,2)=y; 
% [0 fm] The Frequency range for the signal
t=D(:,1);
y=D(:,2);
N=length(y);
dt=t(2)-t(1);% time resolution
fs=1/dt; % the sampling frequency

%fm % the maximum freq

%t=0:dt:(N-1)*dt; % time axis

df=fs/N; % frequency resolution
M=fix(fm/df); % the maximum index for frequency axis
f=0:df:(M-1)*df;

ya=hilbert(y);% get the complex signal
yt=flipud(conj(ya)); % flip the conjugate of the complex signal
[mx,l]=max(yt.*ya);
%pause
k=1; % time variable index
n=0:N-1;n=n';
no=-N/2+1; % the time-shift index

while no<=N/2
   yno=shift2(ya,-no); % shift left by no
   ytno=shift2(yt,+no); % shift right by no
   PSD=(1/2)*(1/N)*fft(yno.*ytno).*exp(j*2*pi*n*(l-1)/N);   % X = (1/sqrt(N))*abs(fft(x));
   Y(:,k)=abs(PSD(1:M)).*sign(real(PSD(1:M)));
   no=no+1;
   k=k+1; % time index
end
%return 
T=t;
F=f/2; % scale down the frequency variable by 2
Z=Y;
%Z=abs(Y);
%Z=Z/max(max(Z));
%figure
%subplot(211)
%mesh(T,F,Y)
%subplot(212)
%contour(T,F,abs(Y))


function ys=shift2(y,no)
N=length(y);
ys=zeros(N,1);
if no<0 % shift left
   lo=abs(no);
   ys(1:N-lo)=y(lo+1:N);
else % shift right
   lo=abs(no);
   ys(lo+1:N)=y(1:N-lo);
end