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

    这一篇对上一篇中各数据矩阵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。其余波形周期暂就不数了。

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

 

 

（本文首发于：http://blog.sina.com.cn/s/blog_6ad0d3de0100op1f.html

首发时间：2011-01-06 15:04:03）