DCCA互相关系数 实例 python

 数据来源：

PhysioBank ATM

下面资源还在审核的话，就用上面的链接下载

MIT- BIH 正常窦性心律数据 nsrdb-16265

---

针对数据ECG1和ECG2进行DCCA分析。

绿色为ECG1和ECG2基于DCCA的线性拟合

【需要去掉注释才能显示】

红色为ECG1基于DFA的线性拟合

蓝色为ECG2基于DFA的线性拟合

---

需要注意如下：

1. 数据以对数曲线表示。

2. 这里时间长time=3 间隔s=0.03 数据长度只取到3000 为了避免错误，请确保【数据长度>=10^time】

3. 基于FDA算ECG1和ECG2的DCCA，整体和DFA代码结构一致。

4. 拟合是线性的，可根据自己数据情况考虑更换二次（比如金融数据）及以上。具体方法：DFA的第三个参数：线性=1，二次=2，以此类推

5. 由于在sqrt开方负数会存在虚部，因此防止np.sqrt时出现nan无效数据，这里采用cmath.sqrt保留虚数，与正常DFA里的sqrt有区别。

---

DFA实例在这 DFA去趋势波动分析 实例 python_Blossom Flight的博客-CSDN博客

以下DCCA代码：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import cmath

# def DFA(data,ni,fittime):
#     n = len(data)//ni
#     nf = int(n*ni)
# 
#     n_mean =np.mean(data[:nf])
#     y = []
#     y_hat = []
#     for i in range(nf):
#         y.append(np.sum(data[:i+1]-n_mean))
#     for i in range(int(n)):
#         x = np.arange(1,ni+1,1)
#         y_temp = y[int(i*ni+1)-1:int((i+1)*ni)]
#         coef = np.polyfit(x,y_temp,deg=fittime)
#         y_hat.append(np.polyval(coef,x))
#     fn = np.sqrt(sum((np.asarray(y)-np.asarray(y_hat).reshape(-1))**2)/nf)
#     return fn

def DCCA(data1,data2,ni,fittime):
    n = len(data1) // ni
    nf = int(n * ni)

    n1_mean = np.mean(data1[:nf])
    n2_mean = np.mean(data2[:nf])
    y1 = []
    y2 = []
    y1_hat = []
    y2_hat = []
    for i in range(nf):
        y1.append(np.sum(data1[:i+1]-n1_mean))
        y2.append(np.sum(data2[:i+1]-n2_mean))
    for i in range(int(n)):
        x = np.arange(1,ni+1,1)
        y1_temp = y1[int(i*ni+1)-1:int((i+1)*ni)]
        coef1 = np.polyfit(x,y1_temp,deg=fittime)
        y1_hat.append(np.polyval(coef1,x))
        y2_temp = y1[int(i * ni + 1) - 1:int((i + 1) * ni)]
        coef2 = np.polyfit(x, y2_temp, deg=fittime)
        y2_hat.append(np.polyval(coef2, x))
    fn = cmath.sqrt(sum((np.asarray(y1) - np.asarray(y1_hat).reshape(-1))*(np.asarray(y2) - np.asarray(y2_hat).reshape(-1))) / nf)
    return fn


df = pd.read_csv('nsrdb-16265.csv')
ECG = df['ECG1'].values[:3000]

time = 3
s = 0.03
t = np.arange(1,time,s)
n = []
for i in range(len(t)):
    n.append(int(10**t[i]))
# f = []
# for i in range(len(n)):
#     f.append(DFA(ECG,n[i],1))
# coef = np.polyfit(np.log10(n),np.log10(f),1)
# y_hat = np.polyval(coef,np.log10(n))
# for i in range(len(n)):
#     y_hat[i] = 10**y_hat[i]
# plt.plot(n,y_hat,c='r')
# 
# 
ECG2 = df['ECG2'].values[:3000]
# 
# f2 = []
# for i in range(len(n)):
#     f2.append(DFA(ECG2,n[i],1))
# coef2 = np.polyfit(np.log10(n),np.log10(f2),1)
# y2_hat = np.polyval(coef2,np.log10(n))
# for i in range(len(n)):
#     y2_hat[i] = 10**y2_hat[i]
# plt.plot(n,f2,'o')
# plt.plot(n,y2_hat,'b')

f_dcca = []
for i in range(len(n)):
    f_dcca.append(DCCA(ECG,ECG2,n[i],1))
coef_dcca = np.polyfit(np.log10(n),np.log10(f_dcca),1)
ydcca_hat = np.polyval(coef_dcca,np.log10(n))
for i in range(len(n)):
    ydcca_hat[i] = 10**ydcca_hat[i]
plt.plot(n,ydcca_hat,'g')
plt.xlabel('n')
plt.ylabel('F(n)')
plt.xscale('log')
plt.yscale('log')
plt.show()