File:Bartl.png

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

Nsub = 100 # number of subsequences

w_1 = 30 # frequency of the 1st component of the signal (Hz)
w_2 = 40 # frequency of the 2nd component of the signal (Hz)

a = 0.7 # magnitude of the 1st component of the signal
b = 0.4 # magnitude of the 2nd component of the signal

t = np.array([i for i in range(1,301)])/1000 # time samples (s)
fs = 1 / (t[1]-t[0]) # sampling frequency (Hz)

x = a*np.cos(2*np.pi*w_1*t) + b*np.sin(2*np.pi*w_2*t) # considered signal

y_mat = np.dot(np.ones((Nsub, 1)), x.reshape((1, len(x)))) # assume that subsequences are identical

Pxx = np.empty((Nsub, int((len(x)/2)+1)))
for i in range(np.shape(y_mat)[0]):
    y_mat[i,:] = y_mat[i,:] # + 2*np.random.randn(len(t))
    f, Pxx[i,:] = signal.periodogram(y_mat[i,:],\
                                     fs=fs, scaling='spectrum')
Pxx_bart = np.mean(Pxx, axis=0)

plt.subplots(1, 1, figsize=(6, 4), dpi=250)
plt.stem(f[1:20], Pxx_bart[1:20], '-')#, linewidth=2, color='b')
plt.ylabel('Spectrum')
plt.xlabel('Frequencies (Hz)')
plt.title('Bartlett\'s method')
plt.grid(True)
plt.show()