Practical cross-spectrum estimation using Blackman-Tukey approach
Let $$\hat{\phi}_{yu}(\omega) = \sum_{k = -M}^{M}  w(k) \ \hat{r}_{yu}(k) e^{-i\omega k}\tag{1}$$
apply a change of variable $$k'=k+M+1\Rightarrow k=k'-M-1$$

and $(1)$ becomes

$$\begin{align}\hat{\phi}_{yu}(\omega) &amp;= \sum_{k' = 1}^{N}  w(k'-M-1) \ \hat{r}_{yu}(k'-M-1) e^{-i\omega (k'-M-1)}\\[10pt]
&amp;=e^{i\omega (M+1)}\sum_{k' = 1}^{N}  w(k'-M-1) \ \hat{r}_{yu}(k'-M-1) e^{-i\omega k'}\\
&amp;=e^{i\omega (M+1)}\sum_{k' = 1}^{N}  w'(k') \ \hat{r}'_{yu}(k') e^{-i\omega k'}
\end{align}$$

where $N=2M+1$ and $w'$ and $\hat{r}'_{yu}(k')$ are the shifted vectors as you did in your own try.



The above answer is when you want to use linear shift. If you wan to use circular shift, then do it like the following:     

 x=w.*rhat;
 phi=fft(circshift(x',M+1)'); 


where M is even and w and rhat are just the symmetric vectors around M+1 (for instance w=[2 4 6 4 2] is symmetric around M=3).