对于做科研工作或工程应用的人来说，拿到一列或数列随机信号观测序列，首先想到的就是看看它们的自相关、互相关函数，自功率谱、互功率谱密度函数估计吧。

    这一篇对上一篇中各数据矩阵TW4615_12、TW2672_12、GY2672_12、DY2672_12、MB2672_12数据作一下相关函数估计，并给出其图形。

MATLAB程序如下(注:为了以示区别,本博客各博文中的程序(包括源代码及注释)全部采用紫色)。

 

%先将各矩阵转换为零均值（时间平均）的一维向量

TW461512_0=TW4615_12';4615天体温

TW461512_0=TW461512_0(:);

TW461512_0=TW461512_0-mean(TW461512_0);

subplot(2,2,1)%绘制子图
plot(TW461512_0)

 

TW267212_0=TW2672_12';%TW4615_12后半部分2672天体温

TW267212_0=TW267212_0(:);

TW267212_0=TW267212_0-mean(TW267212_0);

 

GY267212_0=GY2672_12';%高压

GY267212_0=GY267212_0(:);

GY267212_0=GY267212_0-mean(GY267212_0);

subplot(2,2,2)
plot(GY267212_0)

 

DY267212_0=DY2672_12';%低压

DY267212_0=DY267212_0(:);

DY267212_0=DY267212_0-mean(DY267212_0);

subplot(2,2,3)
plot(DY267212_0)

 

MB267212_0=MB2672_12';%脉搏

MB267212_0=MB267212_0(:);

MB267212_0=MB267212_0-mean(MB267212_0);

subplot(2,2,4)
plot(MB267212_0)

 

    运行，得：

 

                图3-1 各信号原始数据向量图

%计算各列向量自(互)相关函数并画图.由于数据都已经0均值化,故相关函数等价于协方差函数.

R_TW4615=xcorr(TW461512_0,'unbiased')

plot(R_TW4615)

 

    运行，得：

                图3-2 TW461512_0的自相关函数R_TW4615图形

    R_TW4615的长度为55380*2-1=110759（采样点单位为两个小时即一个时辰）。数一数图形包络线上面的“锯齿”，中心线左右大略各为12个多，因此一个“锯齿”代表一年应该没什么问题。

    由于数据比较长，所以看起来黑乎乎的一大块。将数据R_TW4615导入信号处理工具sptool，打开信号观察器窗口Signal Browser，并放大横轴，得图3-2。

                                   图3-3 R_TW4615横轴放大

     可以看出，位于两个波峰处的两根标尺，当中的波形数量为10个，两根标尺之间的距离是dx＝120。因此此图中一个波形的周期是12，即1天，是没有什么问题的。当然，图3-2、图3-3的波形周期都只是凭肉眼观测。严格地说，要确定一列信号中的各频率周期，还要依靠一套算法才行。这个以后再讨论。

     由于周期函数的相关函数也必定是周期函数，因此R_TW4615的周期成分也应该是TW461512_0的周期成分（不过从逻辑上好像不能这么推论。希望内行的朋友赐教：“反之亦真”的必然性在哪？）。

    从TW461512_0的自相关函数图中很容易就看得出的频率成分就这两个。将信号观察器窗口Signal Browser的横轴随意缩放，可以看到R_TW4615的图形常常会出现一些类似“干涉条纹”的现象，如图3-4就是其中的一个。

  

                图3-4 R-TW4615图形中出现的“干涉条纹”现象

     很容易看出来，图3-4中两根标尺之间有3个波形，两根标尺之间的距离大约为dx＝2712，因此一个波形的周期大约是2712/3＝904。这个到底是不是TW461512_0中的一个周期成分呢？我不知道。留待以后慢慢思考。验证。也请内行朋友赐教。这一类的“干涉条纹”图还有很多，无法尽举。我大致数了的，周期成分还有：  516 372.89 354.32 327.6 323.5 319.2 279 262.15 262 258 188.11 186.44 163.71 138.92 131.14 131.08 81.857 77.454 74.182 67.66 65.077 34.604 33.84……

    上面的数据只是我随意缩放横轴后，用鼠标拖动标尺所测量到的各“干涉条纹”图中波形的周期，虽然其小数点后面保留了几位数字，但这并不代表精度。

     下面贴出血压部分的相关函数图。

 

R_GY2672=xcorr(GY267212_0,'unbiased');

plot(R_GY2672)

 

    运行，得：

                图3-5 GY267212_0的自相关函数R_GY2672图形

 

    仍旧将R_GY2672导入信号处理工具sptool，打开信号观察器窗口Signal Browser，并放大横轴与纵轴，得图3-6

                图3-6 R_GY2672横轴纵轴放大

    从图中可以看出，两根标尺之间有20个波形，标尺之间的距离是dx＝240，因此一个波形的周期是240/20＝12，即一天。

 

R_DY2672=xcorr(DY267212_0,'unbiased');

plot(R_DY2672)

 

    运行，得：

 

                图3-7 DY267212_0的自相关函数R_DY2672图形

    上图的放大图就不贴了。 R_GY2672与R_DY2672在信号观察器窗口Signal Browser中，缩放横轴，也都可以看到“干涉条纹”图，只是测起波形周期来没有X_TW4615那么简明，就略去了。 

    下面贴出脉搏的自相关函数图。

 

R_MB2672=xcorr(MB267212_0,'unbiased');

plot(R_MB2672)

 

    运行，得：

 

                图3-8 MB267212_0的自相关函数R_MB2672图形

 

    下面是信号观察器窗口Signal Browser中R_MB2672横轴放大后得到的“干涉条纹”图。

                                 图3-9 R_MB2672的干涉条纹图 

 

    可以看出其一个波形的周期为dx/18=1848/18=102.6667。其余波形周期暂就不数了。

     本来想把自相关、互相关函数图在这一篇中全部贴出的，发现至此内容已经不少了，互相关部分就放到下一篇再贴吧。

 

 

关键词：生命科学 人体节律 生物钟 时间医学  生理信号 信息化 系统医学 系统工程 开放巨系统 混沌理论 复杂性科学 时空还原论 时空整体观 现代信号处理  高阶统计 非平稳信号 循环平稳  非线性 系统辨识 状态估计 模式识别 人工智能 智能科学 控制工程 生物医学工程

 

中医 医理 养生 黄帝内经 难经 经络 气脉 十二经脉 阴阳 五行 五运六气 八字 紫薇 命理 天干 地支 河图  洛书 八卦 六十又四卦  易经 易学 医易  象数 术数  皇极经世  星座 星相  佛教 佛道 道家 丹道  内丹 周易参同契 气功  玄学  神秘学 天人合一  南怀瑾 中国传统文化