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

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).

• Thank's alot for your answer! Just to be clear, with "shifted vectors" do you mean a circular shift of the original vector so that the first M values are taken out and shifted to the end, so that the shifted vector starts with the value corresponding to M+1th value in the original vector? Or maybe something else? Now, practically speaking, if I were to calculate the part inside the sum notation with fft should I scale the result with a complex exponent $e^{i \omega (M+1)}$? – vvv Oct 7 '16 at 10:37
• You're welcome @vvv. We need $w'(k)=w(k-M-1)$ and $\hat{r}'_{yu}(k)=\hat{r}_{yu}(k-M-1)$. It means, both $w(k)$ and $\hat{r}_{yu}(k)$ should be shifted from $[-M, M]$ to $[1, 2M+1]$ such that the start of both are from $1$ rather than $-M$. When you did that, the inside of sum is consistent with the definition of FFT and you can calculate it. At the end, scale it with $e^{i\omega(M+1)}$ – msm Oct 7 '16 at 10:47
• Thanks again. I'm almost there; I circularly shifted both vectors and ran the FFT. What I am wondering though is the scaling term. Consider, for example, that we are dealing with an autocorrelation function, i.e. $\hat{r}_{yy}$. My understanding is that in this case the periodogram values should be real-valued. This indeed is the case after running the FFT (up to very small imaginary terms due to numerical imprecisions, I suppose). However, if I scale these results with $e^{i \omega (M+1)}$ the results are obviously complex valued figures. Any idea what's up with this? – vvv Oct 7 '16 at 18:18
• @vvv The answer I provided will require you to use linear shift, not circular 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). – msm Oct 8 '16 at 3:03
• Okay got it now, many thanks for your great help! Just a thought, it might be helpful to add your previous comment to your original answer as well, at least for me it did the trick (was confused about how to shift vectors). – vvv Oct 8 '16 at 9:21