Some doubts in the matlab code of computing Wigner distribution I found that a discrete form of the Wigner distribution was
$$WD(n,k) = \sum_{m = -N/2}^{N/2}f(n+m)f^*(n-m)e^{-j\frac{2\pi}{N+1}2mk}$$
where $f(n)$ is the signal and time limited within $|n|\le N/2$. I manage to find the matlab code of the Wigner distribution (Sorry that I am not very familiar with the format to put the codes):

 1. function [WD] = EVD_tfrwv(x0)
 2. X=fft(x);
 3. X=[X(1:N/2+1);zeros(N,1);X(N/2+2:N+1)];
 4. x=2*ifft(X);
 5. x=[zeros(N,1);x;zeros(N,1)];
 6. X=zeros(2*N+1);
 7. for k=1:2*N+1   X(:,k)=x(k+(0:2*N));  end 
 8. ww=X.*conj(flipud(X));
 9. WW=fft(ww([N+1:2*N+1,1:N],:)); 
 10. WW=real(WW([N+2:2*N+1,1:N+1],:));


My question is


I understand line 8 tries to get the correlation $f(n+m)f^*(n-m)$ (probably) for different time $n$, but I get confused when line 9 performs the FFT after rearranging the order of these correlation samples. Could anyone tell me why we have to do this?
Line 10 has similar operation.
Why there is a real operator in line 10? I do not see any real operator in the Wigner distribution formula although I understand Wigner distribution is a real distribution. 


Thanks for anybody who tries to help!
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.