MP9: Heart Rate Variability

This MP will analyze electrocardiogram (ECG) signals using NeuroKit2.

1. Installing neurokit2 and loading sample data
2. Extracting the HRV signal
3. Spectral analysis using Periodogram
4. Spectral analysis using LPC
5. Spectral analysis using RLS
6. Verifying your work

1. Installing NeuroKit2 and loading sample data

NeuroKit2 seems to work best with python3.13. Type python --version in your command shell. If python --version reports version 3.13, then you should just be able to type jupyter notebook, and then continue with the rest of this notebook. If it reports a version other than python3.13, you can create a virtualenv containing python 3.13 and neurokit2 using this code:

python3 -m venv -p 3.13 nk
source nk/bin/activate  # on Windows, use nk\Scripts\activate
pip install ipykernel
python -m ipykernel install --user --name=nk --display-name="Python (nk)"
jupyter notebook

Then, in the Kernel menu at the top of this page, choose Change Kernel, then choose Python (nk). Once you've done that, you should be able to run the following blocks:

!pip install neurokit2

Requirement already satisfied: neurokit2 in /opt/anaconda3/lib/python3.13/site-packages (0.2.13)
Requirement already satisfied: matplotlib>=3.5.0 in /opt/anaconda3/lib/python3.13/site-packages (from neurokit2) (3.10.6)
Requirement already satisfied: numpy>=2.0.0 in /opt/anaconda3/lib/python3.13/site-packages (from neurokit2) (2.3.5)
Requirement already satisfied: pandas<3.0.0 in /opt/anaconda3/lib/python3.13/site-packages (from neurokit2) (2.3.3)
Requirement already satisfied: pywavelets>=1.4.0 in /opt/anaconda3/lib/python3.13/site-packages (from neurokit2) (1.9.0)
Requirement already satisfied: requests in /opt/anaconda3/lib/python3.13/site-packages (from neurokit2) (2.32.5)
Requirement already satisfied: scikit-learn>=1.0.0 in /opt/anaconda3/lib/python3.13/site-packages (from neurokit2) (1.7.2)
Requirement already satisfied: scipy in /opt/anaconda3/lib/python3.13/site-packages (from neurokit2) (1.16.3)
Requirement already satisfied: setuptools<82.0.0 in /opt/anaconda3/lib/python3.13/site-packages (from neurokit2) (80.9.0)
Requirement already satisfied: python-dateutil>=2.8.2 in /opt/anaconda3/lib/python3.13/site-packages (from pandas<3.0.0->neurokit2) (2.9.0.post0)
Requirement already satisfied: pytz>=2020.1 in /opt/anaconda3/lib/python3.13/site-packages (from pandas<3.0.0->neurokit2) (2025.2)
Requirement already satisfied: tzdata>=2022.7 in /opt/anaconda3/lib/python3.13/site-packages (from pandas<3.0.0->neurokit2) (2025.2)
Requirement already satisfied: contourpy>=1.0.1 in /opt/anaconda3/lib/python3.13/site-packages (from matplotlib>=3.5.0->neurokit2) (1.3.3)
Requirement already satisfied: cycler>=0.10 in /opt/anaconda3/lib/python3.13/site-packages (from matplotlib>=3.5.0->neurokit2) (0.11.0)
Requirement already satisfied: fonttools>=4.22.0 in /opt/anaconda3/lib/python3.13/site-packages (from matplotlib>=3.5.0->neurokit2) (4.60.1)
Requirement already satisfied: kiwisolver>=1.3.1 in /opt/anaconda3/lib/python3.13/site-packages (from matplotlib>=3.5.0->neurokit2) (1.4.9)
Requirement already satisfied: packaging>=20.0 in /opt/anaconda3/lib/python3.13/site-packages (from matplotlib>=3.5.0->neurokit2) (25.0)
Requirement already satisfied: pillow>=8 in /opt/anaconda3/lib/python3.13/site-packages (from matplotlib>=3.5.0->neurokit2) (12.0.0)
Requirement already satisfied: pyparsing>=2.3.1 in /opt/anaconda3/lib/python3.13/site-packages (from matplotlib>=3.5.0->neurokit2) (3.2.5)
Requirement already satisfied: six>=1.5 in /opt/anaconda3/lib/python3.13/site-packages (from python-dateutil>=2.8.2->pandas<3.0.0->neurokit2) (1.17.0)
Requirement already satisfied: joblib>=1.2.0 in /opt/anaconda3/lib/python3.13/site-packages (from scikit-learn>=1.0.0->neurokit2) (1.5.2)
Requirement already satisfied: threadpoolctl>=3.1.0 in /opt/anaconda3/lib/python3.13/site-packages (from scikit-learn>=1.0.0->neurokit2) (3.5.0)
Requirement already satisfied: charset_normalizer<4,>=2 in /opt/anaconda3/lib/python3.13/site-packages (from requests->neurokit2) (3.4.4)
Requirement already satisfied: idna<4,>=2.5 in /opt/anaconda3/lib/python3.13/site-packages (from requests->neurokit2) (3.11)
Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/anaconda3/lib/python3.13/site-packages (from requests->neurokit2) (2.5.0)
Requirement already satisfied: certifi>=2017.4.17 in /opt/anaconda3/lib/python3.13/site-packages (from requests->neurokit2) (2026.6.17)

import neurokit2 as nk
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

data_all = nk.data('bio_resting_8min_200hz.json')

print('Subjects in the dataset include',data_all.keys())
print('Information available for subject S01 includes',data_all['S01'].keys())

Subjects in the dataset include dict_keys(['S01', 'S02', 'S03', 'S04'])
Information available for subject S01 includes Index(['ECG', 'RSP', 'PhotoSensor', 'Participant'], dtype='object')

fig,axs=plt.subplots(2,1,figsize=(14,4),layout='tight')
time_axis = np.arange(len(data_all['S01']['ECG']))/200
axs[0].plot(time_axis, data_all['S01']['ECG'])
axs[0].set_title('S01 ECG')
axs[1].plot(time_axis[:5000], data_all['S01']['ECG'][:5000])
axs[1].set_title('First 5000 samples from S01 ECG')
axs[1].set_xlabel('Time (seconds)')

Text(0.5, 0, 'Time (seconds)')

Neurokit2 already includes implementations of all of the HRV information that we'll be computing in this MP, and more.

First, let's clean the signals, and return the result as a pandas frame:

clean_signals, info = nk.bio_process(ecg=data_all['S01']["ECG"],rsp=data_all['S01']["RSP"],sampling_rate=200)
print('The clean signals are a pandas frame:\n\n')
clean_signals

The clean signals are a pandas frame:

	ECG_Raw	ECG_Clean	ECG_Rate	ECG_Quality	ECG_R_Peaks	ECG_P_Peaks	ECG_P_Onsets	ECG_P_Offsets	ECG_Q_Peaks	ECG_R_Onsets	...	RSP_Rate	RSP_RVT	RSP_Phase	RSP_Phase_Completion	RSP_Symmetry_PeakTrough	RSP_Symmetry_RiseDecay	RSP_Peaks	RSP_Troughs	RSA_P2T	RSA_Gates
0	2.394536e-19	0.002427	73.103941	0.864590	0	0	0	0	0	0	...	10.388129	0.362346	NaN	0.0	0.618852	0.305981	0	0	175.0	8.241524
1	1.281743e-02	0.005413	73.103941	0.864590	0	0	0	0	0	0	...	10.388129	0.362369	NaN	0.0	0.618852	0.305981	0	0	175.0	8.241524
2	1.129138e-02	0.005777	73.103941	0.864590	0	0	0	0	0	0	...	10.388129	0.362395	NaN	0.0	0.618852	0.305981	0	0	175.0	8.241524
3	7.629118e-04	0.002105	73.103941	0.864590	0	0	0	0	0	0	...	10.388129	0.362422	NaN	0.0	0.618852	0.305981	0	0	175.0	8.241524
4	-4.119742e-03	-0.004066	73.103941	0.864590	0	0	0	0	0	0	...	10.388129	0.362452	NaN	0.0	0.618852	0.305981	0	0	175.0	8.241524
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
101102	2.288843e-02	0.040943	69.767442	0.977552	0	0	0	0	0	0	...	9.167303	0.207551	NaN	0.0	0.553307	0.415584	0	0	130.0	8.280915
101103	8.087383e-03	0.032447	69.767442	0.977552	0	0	0	0	0	0	...	9.167303	0.207605	NaN	0.0	0.553307	0.415584	0	0	130.0	8.280915
101104	-7.628217e-04	0.023663	69.767442	0.977552	0	0	0	0	0	0	...	9.167303	0.207662	NaN	0.0	0.553307	0.415584	0	0	130.0	8.280915
101105	-1.144401e-02	0.014762	69.767442	0.977552	0	0	0	0	0	0	...	9.167303	0.207722	NaN	0.0	0.553307	0.415584	0	0	130.0	8.280915
101106	-1.861572e-02	0.006199	69.767442	0.977552	0	0	0	0	0	0	...	9.167303	0.207784	NaN	0.0	0.553307	0.415584	0	0	130.0	8.280915

101107 rows × 32 columns

clean_signals.keys()

Index(['ECG_Raw', 'ECG_Clean', 'ECG_Rate', 'ECG_Quality', 'ECG_R_Peaks',
       'ECG_P_Peaks', 'ECG_P_Onsets', 'ECG_P_Offsets', 'ECG_Q_Peaks',
       'ECG_R_Onsets', 'ECG_R_Offsets', 'ECG_S_Peaks', 'ECG_T_Peaks',
       'ECG_T_Onsets', 'ECG_T_Offsets', 'ECG_Phase_Atrial',
       'ECG_Phase_Completion_Atrial', 'ECG_Phase_Ventricular',
       'ECG_Phase_Completion_Ventricular', 'RSP_Raw', 'RSP_Clean',
       'RSP_Amplitude', 'RSP_Rate', 'RSP_RVT', 'RSP_Phase',
       'RSP_Phase_Completion', 'RSP_Symmetry_PeakTrough',
       'RSP_Symmetry_RiseDecay', 'RSP_Peaks', 'RSP_Troughs', 'RSA_P2T',
       'RSA_Gates'],
      dtype='object')

hrv = nk.hrv(info, sampling_rate=200, show=False)
print('The heart rate variability results are a pandas series, with 95 different statistics based on the HRV of S01:\n\n')
hrv

The heart rate variability results are a pandas series, with 95 different statistics based on the HRV of S01:

	HRV_MeanNN	HRV_SDNN	HRV_SDANN1	HRV_SDNNI1	HRV_SDANN2	HRV_SDNNI2	HRV_SDANN5	HRV_SDNNI5	HRV_RMSSD	HRV_SDSD	...	HRV_CMSEn	HRV_RCMSEn	HRV_CD	HRV_HFD	HRV_KFD	HRV_LZC	HRV_Symbolic_EqualProb4_0V	HRV_Symbolic_EqualProb4_1V	HRV_Symbolic_EqualProb4_2LV	HRV_Symbolic_EqualProb4_2UV
0	820.749186	86.425815	25.983925	83.68819	24.774153	84.635308	NaN	NaN	61.739996	61.790382	...	1.3603	2.179743	1.733391	1.756958	3.102599	0.814582	0.24183	0.504902	0.187908	0.065359

1 rows × 95 columns

print('The hrv function will also create a plot summarizing key results, if you like:\n\n')
nk.hrv(info, sampling_rate=200, show=True)

The hrv function will also create a plot summarizing key results, if you like:

	HRV_MeanNN	HRV_SDNN	HRV_SDANN1	HRV_SDNNI1	HRV_SDANN2	HRV_SDNNI2	HRV_SDANN5	HRV_SDNNI5	HRV_RMSSD	HRV_SDSD	...	HRV_CMSEn	HRV_RCMSEn	HRV_CD	HRV_HFD	HRV_KFD	HRV_LZC	HRV_Symbolic_EqualProb4_0V	HRV_Symbolic_EqualProb4_1V	HRV_Symbolic_EqualProb4_2LV	HRV_Symbolic_EqualProb4_2UV
0	820.749186	86.425815	25.983925	83.68819	24.774153	84.635308	NaN	NaN	61.739996	61.790382	...	1.3603	2.179743	1.733391	1.756958	3.102599	0.814582	0.24183	0.504902	0.187908	0.065359

1 rows × 95 columns

2. Extracting the HRV Signal

ECG contains at least three types of signal, that are usually processed separately:

1. Unpredictable noise (for example, from $t=4.5$ to $t=6$ in the signal above) is caused by movement of the ECG sensors, or sometimes by the activity of the person's muscles other than their heart.
2. The shape of each individual heart pulse contains information about the relatve timing, strength, and duration of the contractions of the various muscles of the heart. In a normal healthy heart, these contractions happen in a stereotypic sequence not under higher neural control, therefore the shape of each pulse tells you about the health of the person's heart.
3. The duration of the period from one R-pulse to the next varies depending on the physical tension of the chest cavity, and depending on the activation of the vagal nerve, therefore the timing of these pulses contains information about the health of the person's nervous system.

The "heart rate variability" signal (HRV) is the duration between successive R-pulses, treated as a signal.

The nk.bio_process function returns many signals, including:

1. A cleaned ECG signal: it uses a bandpass filter to eliminate very-low-freqency noise and high-frequency noise.
2. An R-peak signal: it uses a pattern-matcher to find the R-peak (the electrical pulse showing closure of the heart's ventricles) in each heartbeat, and outputs a signal that equals 1.0 on every beat, and 0.0 everywhere else.
3. The heart rate signal: it is equal to one over the timing between successive R pulses, in beats per minute.

fig, axs = plt.subplots(3,1,figsize=(14,9), layout='tight')
axs[0].plot(time_axis[:12000],clean_signals['ECG_Raw'][:12000],time_axis[:12000],clean_signals['ECG_Clean'][:12000])
axs[0].legend(['Raw','Clean'])
axs[0].set_title('Raw and cleaned ECG signals')
axs[1].plot(time_axis[:12000],clean_signals['ECG_R_Peaks'][:12000])
axs[1].set_title('Detected R pulses')
axs[2].plot(time_axis[:12000],clean_signals['ECG_Rate'][:12000])
axs[2].plot(time_axis[:12000],clean_signals['ECG_Rate'][:12000]*clean_signals['ECG_R_Peaks'][:12000])
axs[2].set_title('Heart Rate Variability (HRV) Signal: R pulses per minute')
axs[2].set_xlabel('Time (seconds)')

Text(0.5, 0, 'Time (seconds)')

As you can see, the HRV signal is the frequency of heartbeats, measured in beats/minute.

• At every R-pulse, the HRV signal equals the time since the most recent R-pulse (in minutes).
• In between R-pulses, it is smoothed using spline interpolation.

You can see, also, that the HRV signal is quasi-periodic, with a period equal to the respiratory rate (the rate at which the person is breathing, which in this case is about $12\frac{\text{breaths}}{\text{minute}}\approx\frac{1}{5}\frac{\text{breath}}{\text{second}}\approx 0.2 \text{Hz}$). This is called respiratory sinus arrhythmia. As you breathe in, your heart rate increases; as you breathe out, your heart rate decreases. Tracking RSA is one of the easiest ways to find out when a person is breathing (direct measurement of lung volume requires them to wear a belt around the chest; RSA can be tracked if they just wear an electrode). RSA also carries information about stress: The difference between inspiration and expiration is large if the person is relaxed, but RSA flattens out if the person is stressed. RSA differs from one person to another, and generally gets smaller as you get older. RSA also differs during sleep depending on whether you are dreaming (REM sleep) or not dreaming (NREM sleep). Thus, roughly:

• The frequency of RSA tells you how rapidly the person is breathing
• The amplitude of RSA tells you how relaxed the person is

These correlates have been reported to be accurate even for people who can't tell you about themselves, e.g., infants and people who are asleep, and therefore it is useful to try to extract this information automatically.

3. Extracting RSA from HRV using FFT

First, let's try extracting the frequency and amplitude of RSA using an FFT.

We have an HRV signal sampled at 200 samples/second because that's how often the ECG was sampled. Most of those samples are just interpolation, though: the HRV signal itself is just one over the inter-R-pulse duration, measured once per R pulse. There are typically one to two R pulses per second, therefore it's traditional to downsample the HRV signal to a sampling rate of either one sample per R pulse, or one sample per second. Since the spline interpolation has already lowpass filtered it, further lowpass filtering is not really necessary: We can just decimate.

hrv = np.array(clean_signals['ECG_Rate'][::200])
hrv = hrv - np.average(hrv)
N = len(hrv)

fig, axs = plt.subplots(1,1,figsize=(14,3),layout='tight')
axs.plot(np.arange(N),hrv)
axs.set_title('HRV signal (R-peak frequency, in beats/minute) for subject S01')
axs.set_xlabel('Time (seconds)')

Text(0.5, 0, 'Time (seconds)')

As you can see, your heart rate varies a lot.

• The huge sinusoidal variation at about 0.2Hz corresponds to breathing. This type of variability is broadly categorized as high-frequency (HF) heart-rate variability. The high-frequency band is typically defined to cover the frequencies from 0.15Hz to 0.4Hz. The high-end cutoff, 0.4Hz, is chosen mostly because it's below the Nyquist rate: If we have one good measurement every 1.0 seconds, then the Nyquist rate is 0.5Hz.
• The low-frequency (LF) band has a standard definition from 0.04-0.15Hz (6.7 to 25 seconds/cycle). Variation in this band is most commonly attributed to baroreflex, the reflex that keeps blood pressure constant while the body moves. The plot above shows little variation in the LF range, possibly because the person is not moving.
• There are smaller variations with a period of 150-200 seconds (about 0.005-0.007Hz). The very-low-frequency (VLF) band has a standard definition from 0.003 to 0.04Hz (25 to 333 seconds/cycle). Variation in this band is commonly attributed to reflexes that control the body temperature and hormonal balance.
• If we had a longer signal, then we might be able to see variations corresponding to sleep/wake cycles in the ultra-low-frequency (ULF) range, under 0.003Hz (longer than 333 seconds/cycle).

To see the energy in each of these frequency ranges, let's compute the FFT.

HRV = np.fft.fft(hrv)
freq_axis = np.arange(N)/N
NF = int(0.4*N)

fig, axs = plt.subplots(1,1,figsize=(14,3),layout='tight')
axs.plot(freq_axis[:NF],np.square(np.abs(HRV[:NF])))
axs.set_title('Power spectrum (squared magnitude FFT) of HRV (R-peak frequency), subject S01')
axs.set_xlabel('Frequency (Hz)')

Text(0.5, 0, 'Frequency (Hz)')

Obviously this spectrum has energy

• in the VLF range (a peak at around 0.01Hz),
• in the LF range (a peak at around 0.8Hz), and
• in the HF range (a peak at around 0.16Hz).

But exactly what are the frequencies and amplitudes of those peaks? It's hard to decide, because the spectrum is not smooth. This is not an artifact -- the spectrum is actually quite complicated -- but it makes it hard for us to summarize the HRV using a few numbers that a doctor can understand. In order to create a human-readable summary, it would be nice to smooth the spectrum.

One method for smoothing the spectrogram is called Welch's method, based on the periodogram. The idea is:

1. Decide how much smoothing you want. Choose a window, $w[n]$, that smooths the spectrum that much, with length $L$.
2. Divide the signal into overlapping frames of length $L$, overlapping by $L/2$.
3. Compute the periodogram, $|X(\omega)|^2$, from each frame.
4. Average the periodograms.

For example, we know that the rectangular window has its first null at $\frac{2\pi}{L}\left[\frac{\text{radians}}{s}\right]=\frac{F_s}{L}\left[\text{Hz}\right]$, and the Hamming window has its first null at $\frac{4\pi}{L}\left[\frac{\text{radians}}{s}\right]=\frac{2F_s}{L}\left[\text{Hz}\right]$. Suppose we look at the plot above and decide that we want to smooth the spectrum with a frequency resolution of about $\frac{2F_s}{L}=0.02$Hz; that means we want a window of length $L=100$ samples.

import importlib, submitted
help(submitted.welch)

Help on function welch in module submitted:

welch(signal, window)
    Divide the signal into frames of length L, overlapping by L/2.
    Compute the periodogram of each, and average them.

    @param:
    signal (T,) - an array-like object containing the input signal
    window (L,) - the window to use on each frame

    @return:
    frames (int(2T/L-1),L) - frames of speech
    pgram (int(2T/L-1),L) - magnitude-squared FFT of each frame
    spectrum (L,) - average periodogram

importlib.reload(submitted)
frames, pgram, spectrum = submitted.welch(hrv, np.hamming(100))
print('frames shape is',frames.shape)
print('pgram shape is',pgram.shape)
print('welch shape is',spectrum.shape)

frames shape is (9, 100)
pgram shape is (9, 100)
welch shape is (100,)

import librosa
fig, axs = plt.subplots(2,1,figsize=(14,6),layout='tight')
freq_axis = np.arange(100)/100
NF = 41

axs[0].plot(freq_axis[:NF],pgram[:,:NF].T)
axs[0].set_title('Periodograms of each frame of HRV')
axs[1].plot(freq_axis[:NF],spectrum[:NF])
axs[1].set_title('Welch-averaged periodogram of HRV')
axs[1].set_xlabel('Frequency (Hz)')

Text(0.5, 0, 'Frequency (Hz)')

Now we can clearly see the three spectral peaks! They are:

print('The VLF peak is at frequency',np.argmax(spectrum[:5])/100,'Hz, and has energy %2.2g bpm^2'%(np.amax(spectrum[:5])))
print('The LF peak is at frequency',(5+np.argmax(spectrum[5:15]))/100,'Hz, and has energy %2.2g bpm^2'%(np.amax(spectrum[5:15])))
print('The HF peak is at frequency',(15+np.argmax(spectrum[15:]))/100,'Hz, and has energy %2.2g bpm^2'%(np.amax(spectrum[15:])))

The VLF peak is at frequency 0.0 Hz, and has energy 9.4e+03 bpm^2
The LF peak is at frequency 0.08 Hz, and has energy 1.2e+04 bpm^2
The HF peak is at frequency 0.16 Hz, and has energy 1.4e+04 bpm^2

4. Spectral Analysis Using LPC

The periodogram is a useful spectral smoothing method if you don't know why the peaks are there. Suppose we hypothesize, however, that the peaks are there because the nerve signals driving heart rate variability have some resonance, and we would like to find the amplitude and frequency of the resonant peaks. In that case, it would make sense to use LPC to model the spectrum, and then compute the frequencies and amplitudes of the resonant peaks as the roots of the LPC polynomial.

Remember that the LPC spectrum is given by $$H(\omega)=\frac{G}{A(\omega)}=\frac{G}{\sum_{m=0}^{\text{order}} a_mz^{-m}}$$ evaluated at $z=e^{j\omega}$. The gain, $G$, is the RMS energy of the LPC residual, the zero'th coefficient is $a_0=1$, and the other coefficients are computed by minimizing the energy of the LPC residual.

Since we only care about the frequency and amplitude of the peaks, not the overall amplitude, we'll just set $G=1$. To find information about $M$ peaks, we need $\text{order}\ge 2M$. Let's set order=10 to see what kind of information we can extract. Feel free to modify the order, in this section; changing the order will change what kind of information LPC extracts.

import librosa
order = 10
coefficients = librosa.lpc(hrv, order=order)
print(coefficients)

[ 1.         -1.18282049  0.73600019 -0.23571973  0.05867662  0.03865279
 -0.26152813  0.23028896 -0.06123724  0.11291819 -0.14226199]

omega_axis = np.pi*np.arange(100)/100
magH = np.array( [ np.abs(1/np.sum([ coefficients[m]*np.exp(-1j*m*omega) for m in range(order+1) ])) for omega in omega_axis ])

fig, axs = plt.subplots(3,1,figsize=(14,8),layout='tight')

axs[0].plot(np.arange(len(hrv)//2)/len(hrv),np.abs(np.fft.fft(hrv))[:len(hrv)//2])
axs[0].set_title('FFT Magnitude Spectrum')
axs[1].plot(omega_axis/(2*np.pi),magH)
axs[1].set_title('LPC Magnitude Spectrum')
axs[2].plot(omega_axis/(2*np.pi),20*np.log10(np.maximum(0.0001,magH)))
axs[2].set_title('LPC Level Spectrum (dB)')
axs[2].set_xlabel('Frequency (Hz)')
axs[2].set_ylabel('dB')

Text(0, 0.5, 'dB')

The FFT versus periodogram versus LPC spectra seem to be detecting these peaks:

• In the VLF range, the FFT finds peaks at about 0.01Hz and 0.02Hz, whereas the periodogram and LPC both lump those two peaks together at frequency $f=0$Hz.
• In the LF range, the FFT has a peak at about 0.08Hz, which is extracted by both LPC and the periodogram.
• In the HF range, the FFT has an obvious peak at 0.16Hz that's captured by both LPC and the periodogram. LPC also seems to detect a second peak at around 0.28Hz; that peak is not obvious in the periodogram or FFT spectra because it has such a low amplitude.

The big benefit of LPC is that, unlike the FFT and the periodogram, we don't need any separate peak-picking algorithm. The LPC coefficients directly tell us the frequency and amplitude of the peaks, if we use np.roots to find the roots of the LPC polynomial:

$$H(\omega)=\frac{1}{\sum_{m=0}^{M}a_mz^{-m}}=\frac{1}{\prod_{k=1}^{M}(1-p_kz^{-1})}=\frac{z^{M}}{\prod_{k=1}^{M}(z-p_k)}$$

where $p_k$ is the $k^{\text{th}}$ root of the denominator polynomial, and we're assuming for now that $G=1$. Notice that, since the polynomial is real-valued, any complex roots must come in complex-conjugate pairs; for convenience, let's call those $p_k=p^*_{M+1-k}=e^{\sigma_1+j\omega_1}$.

The frequency response $H(\omega)$ is computed by setting $z=e^{j\omega}$. It reaches its highest amplitude at the frequency $\omega=\omega_n$ for any one of the poles. At the frequency $\omega_n$,

$$\left|H(\omega_n)\right|=\frac{\left|e^{j\omega_n M}\right|}{\prod_{k=1}^{M}\left|e^{j\omega_n}-e^{\sigma_k+j\omega_k}\right|}=\frac{1}{\prod_{k=1}^{M}\left|e^{j\omega_n}-e^{\sigma_k+j\omega_k}\right|}$$

The $n=k$ term in the denominator is smaller than all the others. All of the other terms in the denominator are constants, roughly in the range between 0.01 and 1.99, so if we make the very rough approximation that $0.01\approx 1.99\approx 1$, then we can reasonably approximate the equation above as

$$\left|H(\omega_n)\right|\approx\frac{1}{\left|e^{j\omega_n}-e^{\sigma_n+j\omega_n}\right|}$$

Using the approximation $e^x\approx 1+x$, we can further simplify the equation to

$$\left|H(\omega_n)\right|\approx\frac{1}{\left|1+j\omega_n-1-\sigma_n-j\omega_n\right|}=\frac{1}{\sigma_n}$$

The real part of the pole is thus shown to be one over the amplitude of the pole! This constant is usually called the 3dB half-bandwidth of the pole, because the spectrum drops to $20\log_{10}|H(\omega)|=20\log_{10}|H(\omega_n)|-3\text{dB}$ (i.e., $|H(\omega)|=\frac{1}{\sqrt{2}}|H(\omega_n)|$) at frequencies of $\omega=\omega_n\pm\sigma_n$:

$$\left|H(\omega_n\pm\sigma_n)\right|\approx\frac{1}{\left|e^{j(\omega_n\pm\sigma_n)}-e^{\sigma_n+j\omega_n}\right|} \approx\frac{1}{\left|1+j(\omega_n\pm\sigma_n)-1-\sigma_n-j\omega_n\right|}=\frac{1}{\left|\sigma_n\pm j\sigma_n\right|}=\frac{1}{\sigma_n\sqrt{2}}$$

(We call $\sigma_n$ the half-bandwidth because the bandwidth is $(\omega_n+\sigma_n)-(\omega_n-\sigma_n)=2\sigma_n$.)

We can thus find the frequencies and the bandwidths of the poles by computing the imaginary part and real part of the logarithm of the pole:

• Frequency: $f_k=\frac{F_s\omega_k}{2\pi}=\frac{F_s}{2\pi}\times \Im\left\{\ln p_k\right\}$
• Bandwidth: $b_k=\frac{F_s\sigma_k}{\pi}=-\frac{F_s}{\pi}\times \Re\left\{\ln p_k\right\}$

import importlib, submitted
importlib.reload(submitted)
help(submitted.lpc_roots)

Help on function lpc_roots in module submitted:

lpc_roots(coefficients, sample_rate)
    Calculate the frequencies and bandwidths of the LPC poles, in Hertz.

    @param:
    coefficients (order+1,) - the coefficients of the denominator polynomial of an autoregressive spectrum
    sample_rate (1,) - sampling rate of the signal, in samples/second

    @return:
    frequency (order,) - pole frequencies, in Hertz
    bandwidth (order,) - bandwidths, in Hertz

importlib.reload(submitted)
frequency, bandwidth = submitted.lpc_roots(coefficients, 1)
print('The pole frequencies in Hz are:',frequency)
print('The corresponding bandwidths, in Hz, are:',bandwidth)

The pole frequencies in Hz are: [ 0.5         0.3621031  -0.3621031   0.28273807 -0.28273807  0.16655826
 -0.16655826  0.08037029 -0.08037029  0.        ]
The corresponding bandwidths, in Hz, are: [0.06554805 0.10406696 0.10406696 0.07082351 0.07082351 0.02213422
 0.02213422 0.05433299 0.05433299 0.0524679 ]

Since it's a $10^{\text{th}}$ order polynomial, there are 10 poles. The poles at 0Hz and 0.5Hz don't need to be matched with complex conjugates, because $e^{0}$ and $e^{\pi}$ are both real numbers. The other 8 poles each come in complex-conjugate pairs.

• The VLF pole is at 0Hz has bandwidth 0.05Hz.
• The LF poles at $\pm 0.08$Hz have bandwidth 0.05Hz.
• There are HF poles at $\pm 0.17$Hz with high amplitude (BW 0.02Hz).
• There are HF poles at $\pm 0.28$ and $\pm 0.36$Hz with low amplitude (BW 0.07 and 0.1, respectively).

5. Spectral Analysis using RLS

LPC is great for getting very precise frequencies and bandwidths. The disadvantage of LPC, however, is that it needs a lot of samples to compute a good spectral estimate (number of samples is at least twice the LPC order, and it's better if the number of samples is much larger, as it was in our analysis above).

Recursive Least Squares (RLS) updates an all-pole spectral model once for every sample of the original signal. One sample is not enough to calculate a spectrum; instead, RLS recursively updates the model, once per sample, with most of the weight applied to the most recent sample.

Like LPC, RLS estimates $s[n]=\sum_{m=1}^p w_m[n] s[n-m]$, where $w_m[n]=-a_m[n]$ is the $m^{\text{th}}$ predictor coefficient for $1\le m\le\text{order}$ at time $n$. Like LPC, the coefficient vector $\mathbf{w}[n]=[w_1[n],\ldots,w_M[n]]^T$ is computed as $\mathbf{w}[n]=\mathbf{R}^{-1}[n]\mathbf{v}[n]$, where

• $\mathbf{R}[n]=\mathbb{E}\left[\mathbf{x}[n]\mathbf{x}^T[n]\right]$ is the autocorrelation matrix of vector $\mathbf{x}[n]=[s[n-1],\ldots,s[n-M]]^T$ and
• $\mathbf{v}[n]=\mathbb{E}\left[s[n]\mathbf{x}[n]\right]$ is its cross-correlation vector.

Unlike LPC, RLS computes $\mathbf{R}[n]$ and $\mathbf{v}[n]$ recursively as

$$\mathbf{R}[n]=\mathbf{x}[n]\mathbf{x}^T[n]+\lambda\mathbf{R}[n-1]$$ $$\mathbf{v}[n]=s[n]\mathbf{x}[n]+\lambda\mathbf{v}[n-1]$$

where $\lambda<1$ is a memory weight. Larger values of $\lambda$ mean that $\mathbf{R}[n]$ and $\mathbf{v}[n]$ remember more of the past, which reduces their ability to rapidly respond to changes in the signal, but also makes them less responsive to noise.

Updating $\mathbf{R}[n]$ every time, and then inverting it, tends to cause numerical instability. Instead, with some algebra, we can rewrite the update in terms of a gain vector $\mathbf{g}[n]$ and a precision matrix $\mathbf{P}[n]$, so that $\mathbf{w}[n]$ can be computed without any matrix inversions:

$$\mathbf{g}[n]=\frac{\mathbf{P}[n-1]\mathbf{x}[n]/\lambda}{1+\mathbf{x}^T[n]\mathbf{P}[n-1]\mathbf{x}[n]/\lambda},~~~\mathbf{P}[n]=\frac{1}{\lambda}\left(\mathbf{P}[n-1]-\mathbf{g}[n]\mathbf{x}^T[n]\mathbf{P}[n-1]\right)$$

$$\mathbf{w}[n] = \mathbf{w}[n-1]+\mathbf{g}[n]\left(s[n]-\mathbf{x}^T[n]\mathbf{w}[n-1]\right)$$

Even with this formulation, the recursion can sometimes be unstable. A few extra tricks we can use to force it to be stable:

1. Start the recursion at $n=\text{order}$, so that $\mathbf{x}[n]$ is well-defined.
2. Start the recursion with $\mathbf{P}[\text{order}-1]=\mathbf{I}$, the identity matrix.

import importlib, submitted
importlib.reload(submitted)
help(submitted.rls)

Help on function rls in module submitted:

rls(signal, P_init, memory_weight)
    Calculate the RLS predictor coefficients for each sample of the signal,
    starting from sample n=order.

    @param:
    signal (T,) = signal to be predicted from its own past samples
    P_init (order,order) = initial precision matrix before recursion begins
    memory_weight (1,) = lambda: past autocorrelation is downweighted by this factor

    @return:
    predictor (T-order,order) = RLS predictor weights
      - signal[n+order] is predicted by predictor[n,:]@signal[n:n+order]

importlib.reload(submitted)
predictor = submitted.rls(hrv, np.eye(10), 0.95)
print('Signal length is',len(hrv))
print('Predictor coefficient matrix has shape',predictor.shape)

Signal length is 506
Predictor coefficient matrix has shape (496, 10)

rlsf = np.zeros((496,10))
rlsb = np.zeros((496,10))
sgram = np.zeros((496,100))
for n in range(len(predictor)):
    a = np.concatenate(([1], -predictor[n,:]))
    sgram[n,:] = np.array([20*np.log10(np.abs(1/np.sum([a[m]*np.exp(-1j*m*omega) for m in range(order+1)]))) for omega in omega_axis ])
    rlsf[n,:], rlsb[n,:] = submitted.lpc_roots(a, 1.0)

fig, axs = plt.subplots(3,1,figsize=(14,9),layout='tight')
axs[0].plot(hrv)
axs[0].set_title('HRV signal, in beats/minute relative to average')
axs[1].imshow(sgram.T, aspect='auto', origin='lower')
axs[1].set_title('Spectrogram display computed from RLS predictor coefficients')
axs[1].set_ylabel('FFT bin')
axs[2].plot(np.arange(order,len(hrv)),np.abs(rlsf),'.')
axs[2].set_title('Pole frequencies, in Hertz, computed using RLS from each sample of the signal')
axs[2].set_xlabel('Time (seconds)')

Text(0.5, 0, 'Time (seconds)')

What do you see? Here are a few things that I notice:

• The person being recorded was pretty quiet and relaxed, with a nice average heart rate of about 70bpm and a nice big RSA peak at 0.17Hz (about six seconds/breath), except for some changes about 150, 220, and 330 seconds after the start of the recording. From the data shown above it's impossible to say for sure whether the change is caused by movement of the subject's body, or by some type of sensor noise. If you wanted to investigate, you could plot the raw ECG at those times, zoomed in close enough to see individual heartbeats, to see if the ECG is noisy or clean.
• RLS (and maybe LPC too) has a suspicious tendency to choose peaks that are uniformly spaced around the unit circle (10 peaks = one every 0.1Hz). From the spectrum we know that the peaks at 0.08Hz (LF peak) and 0.17Hz (HF peak) are valid, but we might wonder if the higher-frequency peaks really correspond to any real-world heart-rate-variability phenomenon.

6. Verify your work

You can verify your work by running python grade.py before submitting to the autograder:

!python grade.py

.....
----------------------------------------------------------------------
Ran 5 tests in 0.020s

OK