这一篇来看看各生理信号之间的相干函数。相干函数值大的频率，属人体固有的节律频率的可能性大，反之则小。但是，即使相干函数值很小，也不能完全排除该频率仍然是人体固有的节律频率的可能性。

 

 

%下面将脉搏、收缩压、舒张压、均压、差压信号序列0均值化，并在序列

%前面补“0”至与TWXM4860120同样长度

MB4860120=[zeros(1,length(TWXM4860120)-length(MB291712)),MB291712-mean(MB291712)];

 

GY4860120=[zeros(1,length(TWXM4860120)-length(GY291712)),GY291712-mean(GY291712)];

 

DY4860120=[zeros(1,length(TWXM4860120)-length(DY291712)),DY291712-mean(DY291712)];

 

JY4860120=[zeros(1,length(TWXM4860120)-length(JY291712)),GJY291712-mean(GJY291712)];

 

CY4860120=[zeros(1,length(TWXM4860120)-length(CY291712)),CY291712-mean(CY291712)];

 

TWXM4860120=[TW4860120;MB4860120;GY4860120;DY4860120;JY4860120;CY4860120];%组成矩阵

 

Nfft=2^16;%快速傅立叶长度
Fs=Nfft;%采样频率

window=hanning(2^9*48);%

noverlap=2^9*36;

 

    相干函数估计中的功率谱估计采用的是Welch法，其加窗原则见第6篇博文《功率谱密度函数估计》中的Welch法。

 

C=cell(6,6);%1～6行及列分别代表体温、脉搏、收缩压、舒张压、均压、差压。

for i=1:5
  for j=i+1:6
    C{i,j}=cohere(TWXM4860120(i,:),TWXM4860120(j,:),Nfft,Fs,window,noverlap);%相干函数估计
  end
end

 

figure(1)

subplot(2,3,1)
plot(C{1,2})%
title('体温脉搏相干函数Ctm')
xlim([0 3.5e4])

 

subplot(2,3,2)
plot(C{1,3})%
title('体温收缩压相干函数Ctg')
xlim([0 3.5e4])

 

subplot(2,3,3)
plot(C{1,4})%
title('体温舒张压相干函数Ctd')
xlim([0 3.5e4])

 

subplot(2,3,4)
plot(C{2,3})%
title('脉搏收缩压相干函数Cmg')
xlim([0 3.5e4])

 

subplot(2,3,5)
plot(C{2,4})%
title('脉搏舒张压相干函数Cmd')
xlim([0 3.5e4])

 

subplot(2,3,6)
plot(C{3,4})%
title('收缩压舒张压相干函数Cgd')
xlim([0 3.5e4])

 

    运行，得： 

   图24-1

 

figure(2)
subplot(2,3,1)
plot(C{1,2})%
title('体温脉搏相干函数Ctm')
xlim([0 3.5e4])

 

subplot(2,3,2)
plot(C{1,5})%
title('体温均压相干函数Ctj')
xlim([0 3.5e4])

 

subplot(2,3,3)
plot(C{1,6})%
title('体温差压相干函数Ctc')
xlim([0 3.5e4])

 

subplot(2,3,4)
plot(C{2,5})%
title('脉搏均压相干函数Cmj')
xlim([0 3.5e4])

 

subplot(2,3,5)
plot(C{2,6})%
title('脉搏差压相干函数Cmc')
xlim([0 3.5e4])

 

subplot(2,3,6)
plot(C{5,6})%
title('均压差压相干函数Cjc')
xlim([0 3.5e4]) 

 

    运行，得：

   图24-2

 

 

    可以看出各相干函数大小是不一样的。 看看它们的均值与方差：

 

m=zeros(6,6);

v2=zeros(6,6);

v=zeros(6,6);

for i=1:5
  for j=i+1:6

    m(i,j)=mean(C{i,j});%相干函数均值
    v2(i,j)=var(C{i,j});%相干函数方差

    v(i,j)=v2(i,j)^0.5;%相干函数标准差

  end
end

mv=[m;v];   

 

    mv截图，得： 

   图24-3
 

    上图第1～6行、7～12行分别是均值、标准差。

    如果将各相干函数相加，再将其某阈值以上的函数值对应的自变量频率点选为人体节律频率，那么原各相干函数对入选的频率的影响是不一样的。这样显然不合理。应该让各相干函数以同样程度影响力来影响节律频率的选取。将各相干函数乘以某个系数后再相加也不行，不是太大，就是太小，很难确定加权系数大小的。

    可以先求相干函数的频数，从它的频数函数自变量中确定一个阈值，使得在阈值以上的频率数量相等。将这些频率选为人体节律频率，那么各相干函数对入选的频率的影响程度就是一样的了。

 

N=cell(6,6);
for i=1:5
  for j=i+1:6
    [N{i,j},xout]=hist(C{i,j},[0:0.001:1]);%各相干函数的频数函数
  end
end

 

figure(4)

subplot(2,3,1)
barh(xout,N{1,2})%
title('体温脉搏相干函数之频数ntm')
axis([0 100 -0.1 1.1])

subplot(2,3,2)
barh(xout,N{1,3})%
title('体温收缩压相干函数之频数ntg')
axis([0 100 -0.1 1.1])

subplot(2,3,3)
barh(xout,N{1,4})%
title('体温舒张压相干函数之频数ntd')
axis([0 100 -0.1 1.1])

subplot(2,3,4)
barh(xout,N{2,3})%
title('脉搏收缩压相干函数之频数nmg')
axis([0 100 -0.1 1.1])

subplot(2,3,5)
barh(xout,N{2,4})%
title('脉搏舒张压相干函数之频数nmd')
axis([0 100 -0.1 1.1])

subplot(2,3,6)
barh(xout,N{3,4})%
title('收缩压舒张压相干函数之频数ngd')
axis([0 100 -0.1 1.1])

 

    运行，得：

   图24-4 各相干函数的频数函数

 

figure(5)

subplot(2,3,1)
barh(xout,N{1,2})%
title('体温脉搏相干函数之频数ntm')
axis([0 100 -0.1 1.1])

subplot(2,3,2)
barh(xout,N{1,5})%
title('体温均压相干函数之频数ntj')
axis([0 100 -0.1 1.1])

subplot(2,3,3)
barh(xout,N{1,6})%
title('体温舒张压相干函数之频数ntc')
axis([0 100 -0.1 1.1])

subplot(2,3,4)
barh(xout,N{2,5})%
title('脉搏均压相干函数之频数nmj')
axis([0 100 -0.1 1.1])

subplot(2,3,5)
barh(xout,N{2,6})%
title('脉搏差压相干函数之频数nmc')
axis([0 100 -0.1 1.1])

subplot(2,3,6)
barh(xout,N{5,6})%
title('均压差压相干函数之频数njc')
axis([0 100 -0.1 1.1])

 

    运行，得：

 

   图24-5 各相干函数的频数函数

 

    图24-4与图24-5右下角的子图分别是收缩压与舒张压相干函数的频数统计图、均压与差压相干函数的频数统计图。可见二图有很大的差别。它说明了收缩压与舒张压的相干函数，平均比均压与差压的相干函数大很多。相干函数就是信号序列在频域的相关性。在时域的相关性强，则在频域相关性也强，反之亦然。这与我在第4篇博文《谈谈我对人体血压的看法》中对它们的互相关系数的计算结果是一致的。

    下面对各相干函数，自最大值“1”向下，选取数量相等的频率点：

 

CS=cell(6,6);%各频数函数的累加函数
h=zeros(6,6);%各累加函数自变量的某阈值
F=cell(6,6);%原各相干函数在阈值h以上部分的函数点位置索引
K=Nfft/2^9;%频数函数在阈值h以上部分的频率数量（累加值）。此数视情况可以随时任意改变
for i=1:5
  for j=i+1:6
    CS{i,j}=cumsum(N{i,j});
    [mmin,h(i,j)]=min(abs(CS{i,j}-(Nfft/2-fix(K))));
    F{i,j}=find(C{i,j}>h(i,j)*0.001);
  end
end

 

figure(6)
subplot(2,3,1)
stem(F{1,2})%
title('体温脉搏相干函数Ctm入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,2)
stem(F{1,3})%
title('体温收缩压相干函数Ctg入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,3)
stem(F{1,4})%
title('体温舒张压相干函数Ctd入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,4)
stem(F{2,3})%
title('脉搏收缩压相干函数Cmg入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,5)
stem(F{2,4})%
title('脉搏舒张压相干函数Cmd入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,6)
stem(F{3,4})%
title('收缩压舒张压相干函数Cgd入选频率')
axis([0 K+8 0 Nfft/2+8])

 

    运行，得：

   图24-6

 

figure(7)

subplot(2,3,1)
stem(F{1,2})%
title('体温脉搏相干函数Ctm入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,2)
stem(F{1,5})%
title('体温均压相干函数Ctj入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,3)
stem(F{1,6})%
title('体温差压相干函数Ctc入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,4)
stem(F{2,5})%
title('脉搏均压相干函数Cmj入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,5)
stem(F{2,6})%
title('脉搏差压相干函数Cmc入选频率')
axis([0 K+8 0 Nfft/2+8])

subplot(2,3,6)
stem(F{5,6})%
title('均压差压相干函数Cjc入选频率')
axis([0 K+8 0 Nfft/2+8])

 

    运行，得：

   图24-7

 

    下面将上述各频率向量串合在一起。此串合向量Vf取倒数就得到周期数据，以后可以将其并入《21、人体节律周期初步统计排查》的周期数据向量中一起作统计。

Vf=[F{1,2};F{1,3};F{1,4};F{2,3};F{2,4};F{3,4};F{1,2};F{1,5};F{1,6};F{2,5};F{2,6};F{5,6}];
figure(8)
stem(Vf)
title('各相干函数入选频率串合')

 

    运行，得：

   图24-8

   

    为了保持信息来源的均衡性，收缩压、舒张压对均压、差压的相干函数分析我都没有再作，而加入统计的体温对脉搏的频率向量F{1,2}数加倍。

    频率数据串合向量Vf现在就可以作一下频数统计：

 

[Bf index]=sort(Vf);
figure(9)

[nf,xoutf] = hist(Bf,[0:32:Nfft/2]);

 

subplot(1,2,1)
barh(xoutf,nf)
title('所有入选频率的频数统计')
ylim([0 Nfft/2+8])

 

subplot(1,2,2)
bar(Bf)
title('所有入选频率的排序图')
axis([0 length(Vf)+8 0 Nfft/2+8])

 

    运行，得：

   图24-9

 

    这个方法可以最大限度消除传感器噪声谱，但是不能消除谐波频率。    

 

 

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

首发时间：2012-02-04 11:45:16）