Cardiorespiratory coupling index

[blueskynot]

I am trying to recreate the code from this paper https://pubmed.ncbi.nlm.nih.gov/16268385/ to calculate cardiopulmonary coupling.
I am trying to calculate the cross-spectral density between ECG and respiratory signals and then the coherence, is there a way to do that in neurokit?
Sharing a code with synthetic data for convenience
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

# Step 1: Generate synthetic NN intervals and respiration data
def generate_synthetic_nn_intervals(num_intervals, mean_rr=0.8, std_rr=0.05):
    nn_intervals = np.random.normal(mean_rr, std_rr, num_intervals)
    return nn_intervals

def generate_synthetic_respiration(num_points, frequency=0.25, sampling_rate=2):
    t = np.arange(num_points) / sampling_rate
    respiration = np.sin(2 * np.pi * frequency * t)
    return respiration

# Generate synthetic NN intervals
num_intervals = 3000  # Higher number to ensure enough data after downsampling
synthetic_nn_intervals = generate_synthetic_nn_intervals(num_intervals)

# Create a time series for the NN intervals
time_nn = np.cumsum(synthetic_nn_intervals)

# Generate synthetic respiration signal
num_points_resp = int(time_nn[-1] * 2)  # Match the time duration of NN intervals after downsampling
synthetic_respiration = generate_synthetic_respiration(num_points_resp)

# Step 2: Resample NN intervals to 2Hz
target_time = np.arange(0, time_nn[-1], 0.5)  # 2Hz sampling -> 0.5 second intervals
interp_func = interp1d(time_nn[:-1], synthetic_nn_intervals, kind='linear', fill_value="extrapolate")
nn_intervals_resampled = interp_func(target_time)

# Ensure respiration data matches the length of resampled NN intervals
respiration_resampled = synthetic_respiration[:len(nn_intervals_resampled)]

# Step 3: Break into windows of 1024 data points
window_size = 1024
step_size = 512
num_windows = (len(nn_intervals_resampled) - window_size) // step_size + 1

# Step 4: Calculate cross-spectral density and coherence for each window
f, Pxy_list, Cxy_list = [], [], []

for i in range(num_windows):
    start = i * step_size
    end = start + window_size
    nn_window = nn_intervals_resampled[start:end]
    resp_window = respiration_resampled[start:end]
    
    csd? coherence?

[DominiqueMakowski]

this method is currently not implemented, but it would make sense to implement it I think