上一篇博文求出了各生理信号序列imf分量的平均周期，但是靠单一的周期数据是难以确定它就是人体节律周期的。对已经求出的那些周期数据，单凭眼睛也很难判断2个及以上的数据是否为同一周期数据。这就需要将那些数据统一放在一起并按大小次序排列起来进行观察。

VT=[TydhV_TW(:,1);TydhV_MB(:,1);TydhV_GY(:,1);TydhV_DY(:,1);TydhV_JY(:,1);TydhV_CY(:,1)];%将上篇博文中各分量周期时辰数连接为一个向量。

    由于人体节律周期在各通道信号中只表现为有无、强弱，而其数值是不变的，也不随信号处理方式的不同而改变，所以可以将更多不同处理方式中所获得的周期信息放在一起进行比较。

    这一篇博文中对信号序列再作两种处理，并将所得到的周期数据与上一篇博文中的周期数据合并观察。

    （一）先看看信号序列自相关函数的imf分量平均周期（自相关函数图形参见第三篇博文）。

m函数文件：

function [ TydhV_xTW ] = TydhV_imfxcv( TW486012 )
%向量求自相关、作emd分解，再求各分量平均周期

x_TWTW=xcov(TW486012,'unbiased');
x_TW=x_TWTW(length(TW486012):end);%只取时延不小于0部分

a=ones(1,length(TW486012));
a(end-299:end)=linspace(1,0.2,300);%构造一窗函数以消除自相关函数

%末端数据发散严重问题
x_TW=x_TW.*a;
imf_xTW=emd(x_TW);
%for k=1:size(imf_xTW,1)
%  subplot(ceil(size(imf_xTW,1)/3),3,k)
%  plot(imf_xTW(k,:))
%end%图就不作了
TydhV_xTW=TydhV_imf(imf_xTW);%上篇博文中自写程序
end

    以文件名“TydhV_imfxcv”保存。

%体温

TydhV_xTW = TydhV_imfxcv( TW486012 )

VT=[VT;TydhV_xTW(:,1)];

%脉搏

TydhV_xMB = TydhV_imfxcv( MB486012 )

VT=[VT;TydhV_xMB(:,1)];

%收缩压

TydhV_xGY = TydhV_imfxcv( GY486012 )

VT=[VT;TydhV_xGY(:,1)];

%舒张压

TydhV_xDY = TydhV_imfxcv( DY486012 )

VT=[VT;TydhV_xDY(:,1)];

%均压

TydhV_xJY = TydhV_imfxcv( JY486012 )

VT=[VT;TydhV_xJY(:,1)];

%差压

TydhV_xCY = TydhV_imfxcv( CY486012 )

VT=[VT;TydhV_xCY(:,1)];

    （二）再看看信号序列平方函数的imf分量平均周期

function [ TydhV_TW2 ] = TydhV_imfx2( TW486012 )

TW2=TW486012.^2;
TW2=TW2-mean(TW2);%0均值一下以减少emd运算量
imf_TW2=emd(TW2);
TydhV_TW2=TydhV_imf(imf_TW2);

end

 以文件名“TydhV_imfx2”保存。

%体温

TydhV_TW2 = TydhV_imfx2( TW486012 )

VT=[VT;TydhV_TW2(:,1)];

%脉搏

TydhV_MB2 = TydhV_imfx2( MB486012 )

VT=[VT;TydhV_MB2(:,1)];

%收缩压

TydhV_GY2 = TydhV_imfx2( GY486012 )

VT=[VT;TydhV_GY2(:,1)];

%舒张压

TydhV_DY2 = TydhV_imfx2( DY486012 )

VT=[VT;TydhV_DY2(:,1)];

%均压

TydhV_JY2 = TydhV_imfx2( JY486012 )

VT=[VT;TydhV_JY2(:,1)];

%差压

TydhV_CY2 = TydhV_imfx2( CY486012 )

VT=[VT;TydhV_CY2(:,1)];

stem(VT)

    运行，得：

    图21-1.各信号序列及其自相关函数、平方函数的imf分量的平均周期串联图
 

%将周期数据向量VT按从小到大的次序排列：

[b index]=sort(VT);%index是各数据在原向量VT中的序号

logb=log10(b);%取对数
stem(logb)

    运行，得：

    图21-2.各信号序列及其自相关函数、平方函数的imf分量的平均周期排序图。

    上图纵坐标是周期数据的10的幂指数。从这张图中就可以看出人体生理信号很明显的节律性了。图中的台阶，表示从不同的方法处理不同通道信号所得周期数据的聚集性。在低周期端尤其明显。这是因为，对于一定长度的实验数据来说，周期越短，参见统计的波形数量就越多，统计就越精确。随着实验数据的进一步增加与信号处理方式的增加，我相信长周期端的曲线台阶也会变得越来越明显。如果人体生理没有节律性，得到的周期数据就只是些随机数，没有聚集性。其周期排序图将是一条很光滑的曲线（直线）。

stem(index)

    运行，得：

    图21-3.各周期数据原序号

    通过周期数据原序号索引可以回查排序后每一个周期数据的来源（何种处理方式处理何种信号所得数据）。这个不多说了。

  
%周期数据频数统计

[n,xout] = hist(logb,[0.5:0.04:4.5]);

subplot(1,2,1) 

barh(xout,n)

subplot(1,2,2)

bar(logb)

    运行，得：

    图21-4 各周期数据排序及频数统计图

x=zeros(34,6);
x(1:34,1)=n(1:34);

x(1:34,2)=xout(1:34);
x(1:34,3)=n(35:68);
x(1:34,4)=xout(35:68);
x(1:33,5)=n(69:101);
x(1:33,6)=xout(69:101);

    x截图：

    图21-5 频数及周期幂指数变量截图

    统计频数在5个及以上的周期有12个，其周期幂指数分别为0.7000、0.9400、1.1400、1.1800、1.3800、1.4200、1.6600、1.7400、2.0200、3.3800、3.6200、3.7800。将上面12个幂指数换算成周期：

z=[0.7000 0.9400 1.1400 1.1800 1.3800 1.4200 1.6600 1.7400 2.0200 3.3800 3.6200 3.7800];

y=10.^z;

    运行，得y值：

5.01187233627272
8.70963589956081
13.8038426460288
15.1356124843621
23.9883291901949
26.3026799189538
45.7088189614875
54.9540873857625
104.712854805090
2398.83291901949
4168.69383470336
6025.59586074358

    以上为时辰数。2398.83291901949/12 = 199.9027天。

    这样统计出来的结果，其可靠性就比较高了。当然，这次参加频数统计的周期数据样本还不够多。我以后会将各种信号处理方式所得到的周期数据，添加到其向量VT中，继续这项统计。

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

首发时间：2012-01-10 11:49:45）