# Blackman-Tukey PSD in Python

I am trying to calculate the Blackman-Tuckey (BT) PSD in Python to check my understanding (getting started with signal processing). I have tried making the calculation myself and compare it with Scipy's periodogram and pyspectral CORRELOGRAMPSD. I expected some differences between my result and the periodogram (it is a different estimator after all), but there still is a large difference with CORRELOGRAMPSD (which is supposed to be the same). I will try to provide a minimal working example:

Imports and signal generation:

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from scipy.signal import windows as win

from spectrum.correlation import CORRELATION, xcorr
from spectrum import CORRELOGRAMPSD

# Composite signal
sample_intv = 1
sample_freq = 1/sample_intv
t = np.linspace(1, 1000, 1000)
count = t.size

def make_signal(t, amplitude, period, phase=0):
return amplitude * np.sin(2 * np.pi * t/period - phase)

x_time = (make_signal(t, amplitude=2.0, period=50) +
make_signal(t, amplitude=1.0, period=15) +
make_signal(t, amplitude=0.5, period=5))


Calculate the auto-correlation (so far, everything seems ok, I checked with statsmodels and pyspectral and the result seems to match)

lags = signal.correlation_lags(count, count)
valid = (np.abs(lags) < count -1) & (lags >= 0)
lags = lags[valid]

unscaled_xcorr = signal.correlate(x_time, x_time, mode='full')[valid]
xcorr = unscaled_xcorr/(count - np.abs(lags))


Calculate the BT PSD from the autocorrelation (Stoica, 2005, p.37-38):

def blakman_tukey(xcorr, sample_freq, lags=None, window='boxcar'):

# Calculate peridogram resolution (Stoica, 2005, p.37-38)
f_res = 1/xcorr.size
f_max = np.floor(1/f_res) * f_res
f_full = np.arange(0, f_max, f_res)
over_half = f_full > .5
f_half = f_full[~over_half]

# Frequency-domain representation of the data
weighted_xcorr = xcorr * win.get_window(window, xcorr.size)

if lags is None:
x_freq = np.fft.fft(weighted_xcorr)
else:
xcorr_2d, f_2d = np.meshgrid(weighted_xcorr, f_full)
lags_2d, f_2d = np.meshgrid(lags, f_full)
x_freq = np.sum(xcorr_2d * np.exp(-1j * 2 * np.pi * f_2d * lags_2d/sample_freq), axis=1)

# Calculate the PSD from the frequency-domain representation
psd = (np.abs(x_freq)**2)/(xcorr.size * sample_freq) #
psd = psd[~over_half]
psd[1:] *= 2

return f_half, psd

f_per, x_per = signal.periodogram(x_time)
f_a, psd_a = blakman_tukey(xcorr, sample_freq)
f_b, psd_b = blakman_tukey(xcorr, sample_freq, lags)

psd_spec = CORRELOGRAMPSD(x_time, norm='unbiased', lag=999, NFFT=1000, window='rectangle')
f_spec = np.linspace(0, 1, psd_spec.size)
psd_spec = psd_spec[f_spec <= .5]
f_spec = f_spec[f_spec <= .5]
psd_spec[1:] *=2


Plot

plt.plot(f_per, x_per, label='periodogram')
plt.plot(f_spec, psd_spec, label='BT (pyspectral)')
plt.plot(f_a, psd_a, label='BT (FFT)');
plt.plot(f_b, psd_b, label='BT (calc)');
plt.legend();


I do not really understand why the two smaller peaks are so much smaller compared to pyspectral implementation. Is it an issue with the taper? It is supposed to be like that? I've been scratching my head for days and I cannot quite figure out what may be wrong.

Ignacio

Edit 1: plot when negative lags are allowed

Edit 2: plot with log10 yaxis scaling

• They BT (pyspectral) line appears to be negative, which is strange. There might be some issues with the high lag you've chosen, according to the docs.
– Ash
Commented Mar 26 at 16:54
• I also thought having a negative part was weird, but I do not really know why it happens. I reduced the lag to a tenth of the length of the input vector (100). The peaks are lower, but the negative part is still present. Commented Mar 26 at 17:04
• Why do you take only lags >= 0? You're taking only the positive lags of the autocorrelation sequence, but you'll notice in the Stoica and Moses book it uses both negative and positive lags. Commented Mar 26 at 18:19
• When I use both positive and negative lags (removed the >=0 part), the first peak (from the left) is even larger, while the other two remain smaller than expected. I've edited the original plot to make the plot available Commented Mar 26 at 22:30
• Can you also plot on log10 scale? Commented Mar 27 at 3:09

The reason the scaling is so far off is because you are squaring the correlogram. The definition of the periodogram is (per Stoica and Moses)

$$$$\hat{\phi}_{p} = \frac{1}{N}\lvert\sum_{t=0}^{N-1}y(t)e^{-j\omega t}\rvert^{2}$$$$

The definition of the correlogram in Blackman-Tukey form is

$$$$\hat{\phi}_{k} = \sum_{k=-N+1}^{N-1}w(k)\hat{r}(k)e^{-j\omega k}$$$$

By squaring and scaling the values of the correlogram you are messing up the relative powers.

Additionally, if you want a comparable output to the periodogram from your Blackman-Tukey estimate, you need to use a biased estimate of the autocorrelation sequence. There's a bigger discussion to be had if you are interested, but the basic idea can be seen in the discussion preceding equations 2.2.3 and 2.2.4 in the Stoica and Moses book. Essentially, spectral estimators are more concerned with the structure of the data as opposed to the scaling. The unbiased estimate of the autocorrelation sequence has higher statistical variance in the higher lag estimates, and therefore will have poorer structure, causing issues with spectral estimators based on an unbiased ACS estimate.

Here is the code I used that produced comparable estimates. Also, something is not right within your else conditional inside the BT function. I may have messed something up, but I'm not getting coherent estimates out of that method. It's also not really worth it to debug, just use the FFT and zero pad your window on either side to the appropriate length.

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from scipy.signal import windows as win

from spectrum.correlation import CORRELATION, xcorr
from spectrum import CORRELOGRAMPSD

# Composite signal
sample_intv = 1
sample_freq = 1/sample_intv
t = np.linspace(1, 1000, 1000)
count = t.size

def make_signal(t, amplitude, period, phase=0):
return amplitude * np.sin(2 * np.pi * t/period - phase)

x_time = (make_signal(t, amplitude=2.0, period=50) +
make_signal(t, amplitude=1.0, period=15) +
make_signal(t, amplitude=0.5, period=5))

lags = signal.correlation_lags(count, count)
valid = (np.abs(lags) < count)
lags = lags[valid]

unscaled_xcorr = signal.correlate(x_time, x_time, mode='full')[valid]
corr = unscaled_xcorr/(count - np.abs(lags))
corr_test,lags_test = xcorr(x_time,x_time,999,norm='biased')

def blakman_tukey(xcorr, sample_freq, lags=None, window='boxcar'):

# Calculate peridogram resolution (Stoica, 2005, p.37-38)
#f_res = 1/xcorr.size
f_res = len(xcorr)
#f_max = (f_res-1) / f_res
f_full = np.linspace(0, (f_res-1)/(f_res), f_res)
over_half = f_full > .5
f_half = f_full[~over_half]

# Frequency-domain representation of the data
weighted_xcorr = xcorr * win.get_window(window, len(xcorr))

if lags is None:
x_freq = np.fft.fft(weighted_xcorr)
#x_freq = np.fft.fftshift(x_freq)
else:
xcorr_2d, f_2d = np.meshgrid(weighted_xcorr, f_full)
lags_2d, f_2d = np.meshgrid(lags, f_full)
x_freq = np.sum(xcorr_2d * np.exp(-1j * 2 * np.pi * f_2d * lags_2d/sample_freq), axis=1)

# Calculate the PSD from the frequency-domain representation
psd = (np.abs(x_freq)) #
psd = psd[~over_half]
psd[1:] *= 2

return f_half, psd

f_per, x_per = signal.periodogram(x_time)
f_a, psd_a = blakman_tukey(corr_test, sample_freq)
f_b, psd_b = blakman_tukey(corr_test, sample_freq, 999)

psd_spec = CORRELOGRAMPSD(x_time, norm='unbiased', lag=999, NFFT=1000, window='rectangle')
f_spec = np.linspace(0, 1, psd_spec.size)
psd_spec = psd_spec[f_spec <= .5]
f_spec = f_spec[f_spec <= .5]
psd_spec[1:] *=2

plt.figure(1) plt.plot(f_per, x_per)

plt.figure(2) plt.plot(f_spec, psd_spec)

plt.figure(3) plt.plot(f_a,psd_a)

plt.figure(4) plt.plot(f_b,psd_b)

• Thanks a lot for your answer! I am making some tests with your code (It is taking me a while). I think the problem with the else is that 'lags' can accept both arrays and scalars, but the output is wrong when given the latter Commented Mar 27 at 19:48
• @opengisapprendice okay, let me know how it goes! Commented Mar 27 at 22:46
• I have spent the day reviewing your code, pyspectrum documentation and Stoica's book. I re-did most of the code, making sepate functions for the BT correlogram and the BT periodogram. I would like to ask you what is the appropiate course of action when there are such extensive changes: should I edit the original question? post an answer? I want to do the right thing because you helped me a lot Commented Mar 28 at 16:20
• @opengisapprendice if my answer gave you the guidance you were looking for, accepting and upvoting it is fine! If you feel it deserves another answer, then you can also post that and upvote mine! Commented Mar 28 at 17:04
• Done! I will post the new code soon Commented Mar 28 at 17:07