Auto Correlation for Time Frequency Analysis Given a signal $x(t)$, how do I implement a form of autocorrelation function defined as $a(t,T) = x(t-T)x(t+T)$, where $T$ is an arbitrary constant? 

(a fast implementation would be ideal)

Edit: 
This kind of signal I came across from seeing a "Parametric Symmetric Autocorrelation function", defined as above.

It is used in time-frequency analysis methods like WVD,...etc. 
$ R(t,\tau) = x(t + \frac{\tau}{2})x^{*}(t-\frac{\tau}{2})$

Thus far, I have implemented the steps as below for an example chirp:
but the output of the fft2 at the end is wrong. (not a correct frequency)

At the output of the autocorrelation function (the PSIAF variant):


The final output of the LVD is wrong (should be point like):

Solved: Will look at using some of the already published C codes as per answer below to compute the $ R(t,\tau) $
This line seems wrong:

X = signal(X1_signal_indices).*conj(X1conj_signal_indices); 


shouldn't it be

X = signal(X1_signal_indices).*conj(signal(X1conj_signal_indices)); 


??

Note that there is some C code for implementing the WVD and other distributions here.  That code calculates your $R(t,\tau)$ first before convolving various 2D functions with it and then taking the FFT in order to generate the different distributions.