69、长期精密体温数据自相关函数周期分析

        我在2011年本博刚开写的时候，就写了一篇人体生理信号的自相关函数分析的博文《3、自相关函数估计》，其中还贴了截图。那时候体温数据长度还只有13年，截图信息量还很不够。现在9年时间转眼就快过去，数据长度增加到近22年了。现在再来看看体温数据的自相关函数：

x=TW793312;

figure(90)

plot(x)

title('TW793312')

图69-1 体温数据TW793312趋势图

数据起始时间1998/2，结束时间2019/11.采样时间（间隔）2小时。

x=x(:)';%确保x为行向量

x0=x-mean(x);

Rx0=xcorr(x0,'unbiased');

%加个两端收敛的三角窗函数ww:

Lx=length(x);%原数据长度

LRx=2*Lx-1;%自相关函数长度

ww=ones(1,LRx);%长度等于2*Lx-1的单位矩形窗。

n=500;

nn=2*n;

w=triang(nn);%三角窗

ww(1:n)=w(1:n);

ww(LRx-n+1:end)=w(n+1:end);

Rx=Rx0.*ww;

figure(101)

plot(Rx)

title('TW793312自相关函数(加窗)')

图69-2 体温数据TW793312的自相关函数(加窗)

图中一个“锯齿”相当于大约一年、即365*12个数据量（每天12个采样点）。从（上包络线）正中间向（左）右数到（左）右大峰的最高处，大约有16~17个“锯齿”，但这不能说大峰的周期就是16~17年，因为还必须看曲线的下包络线的形状。

下包络线上一个“凹槽”对应于上包络线的一个锯齿。从下包络线正中间凹槽数起，数到（左）右大峰的最低处是12个凹槽，继续数到大峰的最高处是第19个凹槽。但这同样不能说大峰的周期就是12*2=24年或19年。

从图形形状来看，要推知大峰的准确周期，目前的数据长度还是不够的，因为大峰波形还没有完全“发育成熟”。另外还需要一套准确算法，不能就这么凭肉眼数。不管它数据长度够不够，先计算一下。

目前我所知的最好方法还是对数据作emd分解。先看看：

imfRx=emd(Rx);%imfRx是元胞数组

Lv=length(imfRx);分量个数

figure(110)

for j=1:Lv

subplot(4,3,j)

plot(imfRx{j})

end

title('TW793312自相关函数各imf分量')

图69-3 TW793312自相关函数各imf分量

下面看看各imf分量的周期与能量关系。

for j=1:Lv-1%最后一项“残余项”排除

[mn mx z]=extr(imfRx{j});%求极值点与过零点

min{j}=mn;

max{j}=mx;

zp{j}=z;%过零点

zq(j)=(z(end)-z(1))/((length(z)-1)/2);%周期等于最后一个过零点位置减去

%第一个过零点位置，除以过零点个数（减1，或序号直接相减）的一半；

gl(j)=imfRx{j}*imfRx{j}';%第j个分量的能量

end

绘制各分量周期-能量曲线：

figure(120)

stem(log(zq),log(gl))%imf分量的能量与周期的关系

title('TW793312自相关函数各imf分量周期-能量图')

图69-4 自相关函数emd分解各分量周期-能量关系图

左侧为短周期（高频）成分，右侧为长周期（低频）成分。

这张图有两个显著特点：1，能量在第1与第2分量之间发生了很大幅度的跳跃；2，周期在第6、第7分量与第8、第9分量之间都发生了跳跃，尤以后者为大。

emd分解有一个很大的缺点，就是“过度分解”，将一些有物理意义的信号过度分解成了没有物理意义的信号。这个需要凭经验，及事后验证的方法将一些分量重新合并。

经过对各imf分量的分析，我暂时确定第2~6、第8、第10、11分量属于“过度分解。将它们分别与第1、第7、第9分量合并，记为imf16、imf78、imf911

imf16=imfRx{1};

for j=2:6

imf16=imf16+imfRx{j};

end 

figure(16)

plot(imf16)

title('TW793312自相关函数emd第1~6分量合并')

图69-5 imf16图

前面虽然给出了“周期”的计算公式

zq(j)=(z(end)-z(1))/((length(z)-1)/2)

由于实际当中存在“波形发育不完全”，因此通过对波形的观察，不如舍去两端的一些过零点再来计算周期。imf16的周期计算如下：

计算imf16的周期

图69-6 imf16的前部过零点

图69-7 imf16后部过零点

取定前后两个过零点序号与坐标，

图69-8 imf16的周期计算

最后计算结果，imf16周期=11.9999（时辰）

imf78=imfRx{7}+imfRx{8};

figure(78)

plot(imf78)

title('TW793312自相关函数emd第7、8分量合并')

图69-9 imf78

imf78的周期计算

图69-10 imf78前部过零点

图69-11 imf78后部过零点

取定前后两个过零点序号与坐标，

图69-12 imf78周期计算

最后计算结果，imf78周期=365.8417（天）

imf911=imfRx{9}+imfRx{10}+imfRx{11};

figure(911)

plot(imf911)

title('TW793312自相关函数emd第9、10、11分量合并')

图69-13 imf911

图69-14 imf911过零点

取定前后序号2、5两个过零点及其坐标，

图69-15 imf922周期计算

最后计算结果，imf911周期=19.3429（年）

这是一个对人体有着极其重要影响的节律。这样长周期的节律，这样通过实际测量的方式将其呈现出来，相信本人还是全世界第一人。关于长周期节律，以前在书本上所得到的知识，就是我们传统文化中的“天干地支”。“年干”的周期是10年。“年支”的周期是12年。我这里呈现的周期与“天干地支”周期显然密不可分。

Rx的上包络线实质是由每天最高位（大约酉、戌时）的体温所形成，下包络线实质是由每天最低位（大约卯时）的体温，即所谓“基础体温”所形成。会否是Rx的上包络线长期以10*2=20年周期变化，而下包络线长期以12*2=24年周期变化？现在只能猜想，无法断言。

显然，imf78周期=365.8417（天）对应的是“月干”“月支”周期，

imf16周期=11.9999（时辰）对应的是“时干”“时支”周期。其实这很可能仅为imf1的影响所致。imf1的能量太大了，把imf2~imf6的影响全淹没了。而imf2~imf6很可能对应的是“日干”“日支”的影响，只是它的能量太小，无法一眼看出其周期成分、大小。这留待以后慢慢研究。

正如“一天24小时”、“一年365天”这样的节律，你很难说这是人体的节律还是天地大自然的节律，这个“20年”节律，也很难说单纯就是人体节律而与天地大自然无关，而是二者合二为一。这就是“天人合一”的思想与事实。