CODEWITHSUNDEEP
pythonscipy
Scoodood
Scoodood

Reputation: 659

Python Implementation of Bartlett Periodogram

I am trying to implement Periodogram in Python based on the description from Bartlett's method, and compared the result with those from Scipy, by setting overlap=0, use window='boxcar' (rectangle window). However, my result is off by some scale factor. Can someone points out what was wrong with my code? Thanks

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


def my_bartlett_periodogram(x, fs, nperseg, nfft):    
    nsegments = len(x) // nperseg
    psd = np.zeros(nfft)
    for segment in x.reshape(nsegments, nperseg):
        psd += np.abs(np.fft.fft(segment))**2 / nfft
    psd[0] = 0   # important!!
    psd /= nsegments
    psd = psd[0 : nfft//2]
    freq = np.linspace(0, fs/2, nfft//2)
    return freq, psd

def plot_output(t, x, f1, psd1, f2, psd2):
    fig, axs = plt.subplots(3,1, figsize=(12,15))
    axs[0].plot(t[:300], x[:300])
    axs[1].plot(freq1, psd1)
    axs[2].plot(freq2, psd2)
    axs[0].set_title('Input (len=8192, fs=512)')
    axs[1].set_title('Bartlett Periodogram (nfft=512, zero-overlap, no-window)')
    axs[2].set_title('Scipy Periodogram (nfft=512, zero-overlap, no-window)')
    axs[0].set_xticks([])
    axs[2].set_xlabel('Freq (Hz)')
    plt.show()

# Run
fs = nfft = nperseg = 512
t = np.arange(8192) / fs
x = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*100*t) + np.sin(2*np.pi*150*t)
freq1, psd1 = my_bartlett_periodogram(x, fs, nperseg, nfft)
freq2, psd2 = signal.welch(x, fs, nperseg=nperseg, nfft=nfft, window='boxcar', noverlap=0)
plot_output(t, x, freq1, psd1, freq2, psd2)

Output

Upvotes: 7

Views: 3312

Answers (1)

MB-F
MB-F

Reputation: 23637

TL;DR:

Nothing wrong with the code. But welch returns the power spectral density, which is the power spectrum times fs and it compensates for cutting away half the spectrum by multiplying with 2.

To compensate, psd2 * fs / 2 should be very similar to psd.

According to Wikipedia the calculation of psd seems correct:

  1. The original N point data segment is split up into K (non-overlapping) data segments, each of length M
  2. For each segment, compute the periodogram by computing the discrete Fourier transform (DFT version which does not divide by M), then computing the squared magnitude of the result and dividing this by M.
  3. Average the result of the periodograms above for the K data segments.

So whom shall we trust more, Wikipedia or scipy? I would tend towards the latter, but we can find out for ourselves. According to Parseval's theorem the integral over the squared signal should be the same as the integral over the sqared FFT magnitude. Since the Periodogram is obtained from the squared FFT the theorem should hold approximately.

print(np.mean(y**2))  # 1.499727698431174
print(np.mean(psd))  # (1.4999999999999991+0j)
print(np.mean(psd2))  # 0.0058365758754863788

That's close enough for psd, so let's assume it's correct. But I refuse to believe that scipy should be so blatantly wrong! Let's take a closer look at the documentation and see what they have to say about the scaling argument (emphasis mine):

Selects between computing the power spectral density (‘density’) where Pxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Pxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’

Uh-huh! welch's result is the power spectral density, which means it has units of Power per Hz. However, we compared it against the signal power. If we multiply psd2 with the sampling rate to get rid of the 1/Hz units it's the same as psd. Well, except for a factor 2. This factor is meant to compensate for cutting away half the spectrum. If we set return_onesided=False to get the full spectrum that factor is gone.

Upvotes: 5

Related Questions