|||
这一篇来看看各生理信号之间的相干函数。相干函数值大的频率,属人体固有的节律频率的可能性大,反之则小。但是,即使相干函数值很小,也不能完全排除该频率仍然是人体固有的节律频率的可能性。
%下面将脉搏、收缩压、舒张压、均压、差压信号序列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)
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 11:17
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社