Python: Spectrum's BURG Algorithm and Plotting

Question

I am trying to visualize a frequency spectrum using the BURG algroithm. The data that I am trying to visualize is the distance between heartbeats in milliseconds (e.g: [700, 650, 689, ..., 702]). Time distance is measured from R peak to R peak of next heartbeat.

Now I would like to visualize the frequency band with python's spectrum library (I'm a total noob). The minimum frequency that I am trying to display is 0.0033Hz, so all time differences in my dataset summarized are 5 Minutes long.

My approach was to first take the reciprocal of each value, then multiply by 1000, and then multiply by 60. This should get me the Bpm for each heartbeat.

This is what it looks like: [67.11409396 64.72491909 ... 64.58557589]

Afterwards I use spectrum's burg algorithm to create the PSD. The "data" list contains my BpM for each heartbeat.

AR, rho, ref = arburg(data.tolist(), 7)
PSD = arma2psd(AR, rho=rho, NFFT=1024)
PSD = PSD[len(PSD):len(PSD)//2:-1]
plot(linspace(0, 0.5, len(PSD)), 10*log10(abs(PSD)*2./(1.*pi)))
pylab.legend(['PSD estimate of x using Burg AR(7)'])

The graph that I get looks like this: 5 Minutes Spectrogram

This specific data already exists as a 3D-Spectrogram (Graph above is the equivalent to the last 5 Minutes of 3D-Spectrogram): Long Time 3D-Spectrogram

My Graph does not seem to match the 3D-Spectrogram. My frequencies are way off.... What causes this and how can I fix it?

Also I would like the y-Axis in my Graph not in [dB] but in absolute Values. I tried with:

plot(linspace(0, 0.5, len(PSD)), abs(PSD))

but that did not really seem to work. It just drew a hyperbole.

Thank you for your help!

Answer 1

The spectrum package comes with a pburg class than can generate a frequencies array, this is shown below. If you want direct comparison between a spectrogram and AR PSDs, I would take the time definition used to compute the spectrogram to also compute the AR PSD per window.

Also, your spectrogram example image looks focused on very low frequencies, so you may want to increase nfft to increase frequency resolution.

import matplotlib.pyplot as plt
from scipy.signal import spectrogram
import numpy as np
from spectrum import pburg

# Parameter settings
n_seconds = 10
fs = 1000 # sampling rate, in hz
freq = 10

nfft = 4096
nperseg = fs

order = 8

# Simulate 10 hz sine wave with white noise 
x = np.sin(np.arange(0, n_seconds, 1/fs) * freq * 2 * np.pi)
x += np.random.rand(len(x)) / 10

# Compute spectrogram
freqs, times, powers = spectrogram(x, fs=fs, nfft=nfft)

# Get spectrogram time definition
times = (times * fs).astype(int)
window_times = np.array((times-times[0], times+times[0])).T

# Compute Burg's spectrum per window
powers_burg = np.array([pburg(x[t[0]:t[1]], order=order,
                              NFFT=nfft, sampling=fs).psd for t in window_times]).T

freqs_burg = np.array(pburg(x, order=order, NFFT=nfft, sampling=fs).frequencies())

# Plot
inds = np.where(freqs < 20)
inds_burg = np.where(freqs_burg < 20)

fig, axes = plt.subplots(ncols=2, figsize=(10, 5))

axes[0].pcolormesh(times/fs, freqs[inds], powers[inds], shading='gouraud')
axes[1].pcolormesh(times/fs, freqs_burg[inds_burg], powers_burg[inds_burg], shading='gouraud')

axes[0].set_title('Spectrogram')
axes[1].set_title('Burg\'s Spectrogram')