Practical cross-spectrum estimation using Blackman-Tukey approach

Question

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&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' < 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!

Answer 1

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) &= \sum_{k' = 1}^{N} w(k'-M-1) \ \hat{r}_{yu}(k'-M-1) e^{-i\omega (k'-M-1)}\\[10pt] &=e^{i\omega (M+1)}\sum_{k' = 1}^{N} w(k'-M-1) \ \hat{r}_{yu}(k'-M-1) e^{-i\omega k'}\\ &=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).