Some doubts in the matlab code of computing Wigner distribution
The FFT comes from the summing operation and multiplication by the
complex pointer term in the formula. Compare it to the formula for
the DFT (or FFT):

$$X_k = \sum_{n=0}^{N-1} x_n\cdot e^{-\frac {2\pi j}{N}kn}$$
This is not similar, real just returns the real part of the array.
The Wigner distribution is real, but because of rounding errors, there is always a very small complex part in the result. This operation removes it and just returns the real part.