How do I get the phase angle from Cross Correlation? Im trying to cross correlate two signals in matlab and get the phase difference between the signals.

For cross correlation (the idea is to do it without xcorr) I used:

Cxx=fftshift(ifft(fft(x,N).*conj(fft(y,N))))/(norm(x) * norm(y));


I get the result and not sure about the reference point I have to take for phase calculation Now which point should I take as zero on time scale and how do I get the phase difference from this result. 

Thanks.
There is no such thing as the phase angle between two signals unless they both consist of a single sinusoid at the same frequency, that is, $x(t) = A\cos(\omega t+\psi)$
and $y(t) = B\cos(\omega t + \phi)$.  If you have $N$ samples of these signals $x(t)$ and
$y(t)$, taken at times $0$, $T$, $2T, \ldots$, $(N-1)T$, so that
$$x[n] = x(nT), ~~ y[n] = y(nT), 0 \leq n &lt; N,$$ and either $N\omega T$ is an
integer multiple of $2\pi$ or $N\omega T \gg 1$
then the phase angle between the two sinusoids is
$$\theta 
= \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} x[n]y[n]}{\sqrt{\displaystyle 
\sum_{n=0}^{N-1}(x[n])^2\displaystyle \sum_{n=0}^{N-1}(y[n])^2}}\right).
\tag{1}$$
For complex-valued sinusoids $x(t) = Ae^{j(\omega t+\psi)}$ and
$y(t) = BAe^{j(\omega t+\phi)}$, $(1)$ should be replaced by
$$\theta 
= \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} x[n](y[n])^*}{\sqrt{\displaystyle 
\sum_{n=0}^{N-1}|x[n]|^2\displaystyle \sum_{n=0}^{N-1}|y[n]|^2}}\right).
\tag{2}$$

Of course, these formula can be used for arbitrary signals, not just for pure sinusoids,
but then, what you get is the _angle between the two vectors 
$\mathbf x = (x[0], x[1], \ldots, x[N-1])$ and 
$\mathbf y =(y[0], y[1], \ldots, y[N-1])$ in the $N$-dimensional spaces
$\mathbb R^N$ or $\mathbb C^N$, and not a phase angle between two sinusoids 
at the same frequency.
Note that the three points $\mathbf 0 = (0,0,\ldots,0)$, $\mathbf x$ and
$\mathbf y$ lie in a (two-dimensional) plane in $N$-dimensional space and
the $\theta$ that you get is the angle between the line segments with endpoints
$\mathbf 0$ and $\mathbf x$ and endpoints $\mathbf 0$ and $\mathbf y$ which
lie in this plane. Another way to think about this is that
$$\langle\mathbf x, \mathbf y\rangle = ||\mathbf x||\cdot||\mathbf y||\cdot\cos(\theta)$$
and thus $(1)$ is obtained from
$$\cos(\theta) 
= \frac{\langle\mathbf x, \mathbf y\rangle}{||\mathbf x||\cdot||\mathbf y||}
= \frac{\displaystyle \sum_{n=0}^{N-1} x[n](y[n])^*}{\sqrt{\displaystyle 
\sum_{n=0}^{N-1}|x[n]|^2}\sqrt{\displaystyle \sum_{n=0}^{N-1}|y[n]|^2}}.$$

If you must use FFTs because that's the way you have been told to do it, then you have
$$\theta = \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} X[n](Y[n])^*}{\sqrt{\displaystyle 
\sum_{n=0}^{N-1}|X[n]|^2\displaystyle \sum_{n=0}^{N-1}|Y[n]|^2}}\right)\tag{3}$$
so that you have the ineffable pleasure of not only needing to
compute two FFTs first, but also of using complex multiplications in $(3)$
instead of the real multiplications in $(1)$ (for real-valued
signals).  This is overkill in my estimation,
but as usual, YMMV, and what your boss insists on is always right, regardless
of what people write on Internet forums.