Practical cross-spectrum estimation using Blackman-Tukey approach I would like to estimate the cross-spectrum of two signals using the (lag-windowed) Blackman-Tukey approach but I'm having difficulties with proper practical implementation. As defined in equation 2.8.32 of the the book Spectral Analysis of Signals (2005, Prentice Hall) by Stoica&amp;Moses, the Blackman-Tukey periodogram estimate of the cross-spectrum between two stationary signals $ y(t)$ and $ u(t)$ (for the purpose of my question we may also assume that both signals are real-valued) is given by

\begin{align*}
\hat{\phi}_{yu}(\omega) = \sum_{k = -M}^{M}  w(k) \ \hat{r}_{yu}(k) e^{-i\omega k}
\end{align*}

where $ w$ denotes some appropriate lag-window (symmetric with $ 2M'+1$ non-zero observations, $ M' &lt; M$ ) and $ \hat{r}_{yu}$ estimated cross-covariance function with M lags (also in total $ 2M+1$ observations). I try to implement the estimation via Discrete Fourier Transform (DFT) using a typical Fast Fourier Transform (FFT) algorithm available in most statistical packages, say MATLAB's fft function (documentation here). 

What I am not sure about is how to perform the DFT on windowed cross-covariance values since index $k$ is taking both negative and positive values. When using MATLAB's fft function we feed in a vector representing the signal (with some length $N$) we want to transform and the algorithm performs the transform with index $k$ running from $1$ to $N$. As I understand it, if we were dealing with a DFT of an autocovariance function, say $ \hat{r}_{yy}$, then the negative indices would be no problem as $\hat{r}_{yy} $ is symmetric so the "two sides" around $ k = 0$ are essentially the same and one could just perform an "one-sided" DFT by feeding values $ w(m) \ \hat{r}_{yy}(m) \ , \ m = 0,1,...,M$ to the fft algorithm. As symmetry need not be the case with cross-covariances, my understanding is that this approach does not work. My initial solution was just to cram the windowed cross-covariance values 

$$  w(k) \ \hat{r}(k) \ , \ k = -M,...,-1,0,1,...M $$ 

into a vector of length $ 2M + 1$ and feed it to the fft function as it is. However, I fear there are problems with this approach due to exponents of the DFT taking up only negative values. My question is therefore how should one perform the DFT with a vector of (windowed) cross-covariance values in order to obtain the desired estimate of the cross-spectrum? 

Just to be clear I explicitly want to use the Blackman-Tukey / cross-covariance approach and not e.g. Welch's approach in estimation. Help is greatly appreciated!
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).