科学网 › 标签 › 生理信号

标签: 生理信号

相关帖子	版块	作者	回复/查看	最后发表

没有相关内容

相关日志

25、GUI界面一维小波消噪: baishp 2012-11-30 00:01; 前面使用的信号序列都是原观测序列，没有经过任何预处理的。这一篇来做一下这件工作。应用小波消噪之前，先用EMD方法消除一下趋势项。用EMD方法消除趋势项可以说是迄今为止最好的消除趋势项的方法。 TW4860_r=TW486012-imf_TW486012(end,:); MB2917_r=MB291712-imf_MB291712(end,:); GY2917_r=GY291712-imf_GY291712(end,:); DY2917_r=DY291712-imf_DY291712(end,:); JY2917_r=JY291712-imf_JY291712(end,:); CY2917_r=CY291712-imf_CY291712(end,:); TW4860_r的下标“_r”表示原观测序列TW486012去除了emd分解中的残余项。余类推。不过，所得新序列的均值可能并不严格等于0，而是有一点点偏差，这是由emd分解的特点决定的。暂不管这个。 M_r=cell(6,3); M_r{1,1}=TW4860_r; M_r{1,2}='TW4860_ r'; M_r{1,3}= ; M_r{2,1}=MB2917_r; M_r{2,2}='MB2917_ r'; M_r{2,3}= ; M_r{3,1}=GY2917_r; M_r{3,2}='GY2917_ r'; M_r{3,3}= ; M_r{4,1}=DY2917_r; M_r{4,2}='DY2917_ r'; M_r{4,3}= ; M_r{5,1}=JY2917_r; M_r{5,2}='JY2917_ r'; M_r{5,3}= ; M_r{6,1}=CY2917_r; M_r{6,2}='CY2917_ r'; M_r{6,3}= ; for i=1:6 subplot(3,2,i); plot(M_r{i,1});title(M_r{i,2}); axis(M_r{i,3}); end 运行，得：图25-1 用emd方法消除趋势项之后的各生理信号序列下面用小波工具箱图形用户界面GUI对上述信号进行消噪。点击GUI界面左上的Wavelet 1-D调出一维小波分析界面，导入体温序列TW4860_r。先随意选一个小波db4，分解层次选最多层12层。点Analyze按钮开始分析。显示模式Display mode选Separate Mode（分开形式）。见下图。图25-2 上图中，点More Display Options可以选择左右两边的显示内容。我调节的，左边第一列是原始信号TW4860_r，下面12列是各层近似信号；右边第一列是各层细节系数，第二列是重构信号，再以下是各层细节信号。点View Axes可以将上面各图放大看。各层细节系数放大图如下。图25-3 各层细节系数颜色越浅的系数越大。大约可见第三层细节系数较邻近层系数要大。点图25-2中的De-noiss，打开一个消噪界面，见图25-4。图25-4 fixed form阈值消噪界面上有三类参数需要选择：软阈值soft硬阈值hard选择，既然大多数人都认为软阈值作用方式比较好，那我就先选定软阈值；噪声结构选择Select noise structure，选Scaled white noise(这个到底是指“强度非归一化白噪声”，还是指“有色高斯噪声”？我一直没找到确定的答案。)吧；阈值选取方式Select thresholding method有7种。下面分别选前4种来进行消噪，看看消噪的效果如何。先选固定形式阈值Fixed form threshold。选定后，点消噪De-noiss，界面即进行消噪处理。界面左侧给出各层细节系数与阈值线，（中）右上面是原始信号（红色）与消噪后的信号（紫色），中间是原始信号分解后各层细节系数，下面是经阈值作用后各层细节系数。点De-noiss右边的残差Residuals，打开残差显示界面，见下图。图25-5 fixed form阈值消噪之残差上图上方是残差（原信号与消噪信号之差，即被“消除”的噪声）。中左是残差的直方图，就是我以前说的“频数函数”图，是数据总体的概率密度函数估计。中右是频数函数的积分函数，即概率分布函数（估计）。下左图是残差自相关，下右图是残差功率谱。最下面是残差各项统计数据。继续选择阈值方式分别为Rigorous SURE（基于Stein无偏似然估计SURE阈值）、Heuristic SURE（启发式SURE阈值）、Minimaxi（最小极大方差阈值）。分别得到消噪界面及残差图，见下图25-6～图25-11。图25-6 rigorous sure阈值消噪图25-7 rigorous sure阈值消噪之残差图25-8 heuristic sure阈值消噪图25-9 heuristic sure阈值消噪之残差图25-10 minimax阈值消噪图25-11 minimax阈值消噪之残差可以看出图25-8、图25-9与图25-6、图25-7完全是一样的。因为启发式SURE阈值Heuristic SURE，就是在固定形式阈值Fixed form threshold与无偏似然估计SURE阈值Rigorous SURE之间作一个选择。本例中选择的是无偏似然估计SURE阈值，所以两个SURE阈值消噪结果是一样的。接下来Heuristic SURE阈值方式就不再单独分析了。图25-5与图25-11中残差的直方图，看起来符合白噪声的概率密度函数即正态分布。而图25-7（图25-9）的直方图不符合正态分布。从下方的自相关函数与功率谱图中知道它们都含有至少一个周期成分。这个后面再谈。衡量消噪效果有两个原则：消噪信号与原信号的相似性、消噪信号的光滑性。但这两个原则是一对矛盾的统一体：太相似就不会太光滑，太光滑就不会太相似。最终的取舍只能看具体情况，折中选择了。衡量消噪信号与原信号的相似性，有几个衡量指标，如残差的标准差（standard.dev）、残差的绝对偏差（mean abs. dev.）、残差模数（L2 norm）或消噪信号与原信号的模数比。从图25-5、图25-7、图25-9、图25-11下方的统计数据中可以知道，前3项指标分别为：阈值： Fixed SURE Minimaxi 标准差： 15.23、 2.673、 11.97 平均绝对差：12.31、 2.267、9.662 模数： 3679、 645.5、 2890 指数越小表示相似程度越高，否则越低。因此Rigorous SURE（Heuristic SURE）阈值消噪的相似性最好，Minimaxi阈值其次，Fixed form阈值最差。再看看消噪信号与原信号的模数比，其含义类似于信号能量比。将上面第1、第2、第4种阈值消噪所得的消噪信号导出，分别记为TWden_fix、TWden_rig、TWden_mnmx。残差也导出，分别记为TWres_fix、TWres_rig、TWres_mnmx。 %模数比： normr1=norm(TWden_fix)/norm(TW4860_r) %MATLAB打印 normr1 = 0.7205 normr2=norm(TWden_rig)/norm(TW4860_r) %MATLAB打印 normr2 = 0.9846 normr4=norm(TWden_mnmx)/norm(TW4860_r) %MATLAB打印 normr4 = 0.7891 总的感觉，第1、4种方法，消噪“太厉害”，消噪信号与原信号相似性差，但光滑度肯定很好；第2（3）种方法，消噪“太弱”，消噪信号与原信号相似性肯定很好，但光滑度会差。虽然说“太相似就不会太光滑，太光滑就不会太相似”，但不能说“太不相似就会太光滑，太不光滑就会太相似”，因此信号光滑度应该另有一衡量指标才好。目前我还没有看到光滑度计算公式，不过从功率谱的本义，光滑信号低频端功率谱值高，高频端功率谱值低；不光滑信号则与此相反。现在把消噪信号TWden_fix、TWden_rig、TWden_mnmx、残差TWres_fix、TWres_rig、TWres_mnmx都导入信号观测器sptool，并加载于其信号功率谱处理器Spectrum Viewer，从直观上来看看各信号的光滑度。下面图25-12～图25-22是各信号、残差的功率谱图，蓝线原信号，紫线消噪信号，黄线残差。图25-12 固定形式阈值Fixed form threshold的消噪信号功率谱图25-13 固定形式阈值Fixed form threshold的消噪残差(噪声)功率谱从上面两图中看到一个有趣的现象：原信号看不出的一个周期成分，在它的消噪信号与残差中却同时出现了，见上面两图中标尺线左边一点点的谱线。经查，知其基波周期为24个时辰。消噪信号其3倍频、9倍频谱线也很明显，残差功率谱中则只看到3倍频。这个该如何解释呢？我目前也没有一下子想清楚。图25-14 原信号与固定形式阈值Fixed form threshold消噪信号、残差(噪声)三合一功率谱图25-15 图25-14低频端放大图25-16 图25-14高频端放大图25-17 无偏似然估计Rigorous SURE阈值的消噪信号功率谱此图中未见24时辰周期成分基、倍频谱线。图25-18 无偏似然估计Rigorous SURE阈值的消噪残差功率谱此图中只见24时辰周期成分基频谱线。图25-19 原信号与无偏似然估计Rigorous SURE阈值的消噪信号、残差(噪声)三合一功率谱图25-20 最小极大方差阈值Minimaxi的消噪信号功率谱此图又见周期为24个时辰周期成分之基频、3倍频、9倍频谱线。图25-21 最小极大方差阈值Minimaxi的消噪残差功率谱此图可见周期为24个时辰周期成分之基频、3倍频谱线。图25-22 原信号与最小极大方差阈值Minimaxi的消噪信号、残差(噪声)三合一功率谱看看各信号在时域的波形： plot(TW4860_r) hold on plot(TWden_fix,'g-') plot(TWden_rig,'r-') plot(TWden_fix,'k-') xlim( ) 运行，得：图25-23 原信号及3种消噪信号时域波形比较从上面的功率谱图及时域波形光滑性可以看出，各阈值消噪的效果，与消噪相似性指标反映的消噪效果是一致的。感觉，消噪信号的应用，应该灵活掌握。有些信号处理方式中应该用消噪信号，另一些信号处理方式中就不一定要使用消噪信号。因为这些处理方式本身可能就含有消噪功能，用已经消噪的信号，也许会丢失一部分信息。（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de010138ox.html 首发时间：2012-02-24 19:37:48）; 个人分类: 斤斤计较|4407 次阅读|0 个评论

24、基于各生理信号之间相干函数的频率选择: baishp 2012-11-29 23:52; 这一篇来看看各生理信号之间的相干函数。相干函数值大的频率，属人体固有的节律频率的可能性大，反之则小。但是，即使相干函数值很小，也不能完全排除该频率仍然是人体固有的节律频率的可能性。 %下面将脉搏、收缩压、舒张压、均压、差压信号序列0均值化，并在序列 %前面补“0”至与TWXM4860120同样长度 MB4860120= ; GY4860120= ; DY4860120= ; JY4860120= ; CY4860120= ; TWXM4860120= ;%组成矩阵 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( ) subplot(2,3,2) plot(C{1,3})% title('体温收缩压相干函数Ctg') xlim( ) subplot(2,3,3) plot(C{1,4})% title('体温舒张压相干函数Ctd') xlim( ) subplot(2,3,4) plot(C{2,3})% title('脉搏收缩压相干函数Cmg') xlim( ) subplot(2,3,5) plot(C{2,4})% title('脉搏舒张压相干函数Cmd') xlim( ) subplot(2,3,6) plot(C{3,4})% title('收缩压舒张压相干函数Cgd') xlim( ) 运行，得：图24-1 figure(2) subplot(2,3,1) plot(C{1,2})% title('体温脉搏相干函数Ctm') xlim( ) subplot(2,3,2) plot(C{1,5})% title('体温均压相干函数Ctj') xlim( ) subplot(2,3,3) plot(C{1,6})% title('体温差压相干函数Ctc') xlim( ) subplot(2,3,4) plot(C{2,5})% title('脉搏均压相干函数Cmj') xlim( ) subplot(2,3,5) plot(C{2,6})% title('脉搏差压相干函数Cmc') xlim( ) subplot(2,3,6) plot(C{5,6})% title('均压差压相干函数Cjc') xlim( ) 运行，得：图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= ; mv截图，得：图24-3 上图第1～6行、7～12行分别是均值、标准差。如果将各相干函数相加，再将其某阈值以上的函数值对应的自变量频率点选为人体节律频率，那么原各相干函数对入选的频率的影响是不一样的。这样显然不合理。应该让各相干函数以同样程度影响力来影响节律频率的选取。将各相干函数乘以某个系数后再相加也不行，不是太大，就是太小，很难确定加权系数大小的。可以先求相干函数的频数，从它的频数函数自变量中确定一个阈值，使得在阈值以上的频率数量相等。将这些频率选为人体节律频率，那么各相干函数对入选的频率的影响程度就是一样的了。 N=cell(6,6); for i=1:5 for j=i+1:6 =hist(C{i,j}, );%各相干函数的频数函数 end end figure(4) subplot(2,3,1) barh(xout,N{1,2})% title('体温脉搏相干函数之频数ntm') axis( ) subplot(2,3,2) barh(xout,N{1,3})% title('体温收缩压相干函数之频数ntg') axis( ) subplot(2,3,3) barh(xout,N{1,4})% title('体温舒张压相干函数之频数ntd') axis( ) subplot(2,3,4) barh(xout,N{2,3})% title('脉搏收缩压相干函数之频数nmg') axis( ) subplot(2,3,5) barh(xout,N{2,4})% title('脉搏舒张压相干函数之频数nmd') axis( ) subplot(2,3,6) barh(xout,N{3,4})% title('收缩压舒张压相干函数之频数ngd') axis( ) 运行，得：图24-4 各相干函数的频数函数 figure(5) subplot(2,3,1) barh(xout,N{1,2})% title('体温脉搏相干函数之频数ntm') axis( ) subplot(2,3,2) barh(xout,N{1,5})% title('体温均压相干函数之频数ntj') axis( ) subplot(2,3,3) barh(xout,N{1,6})% title('体温舒张压相干函数之频数ntc') axis( ) subplot(2,3,4) barh(xout,N{2,5})% title('脉搏均压相干函数之频数nmj') axis( ) subplot(2,3,5) barh(xout,N{2,6})% title('脉搏差压相干函数之频数nmc') axis( ) subplot(2,3,6) barh(xout,N{5,6})% title('均压差压相干函数之频数njc') axis( ) 运行，得：图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}); =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( ) subplot(2,3,2) stem(F{1,3})% title('体温收缩压相干函数Ctg入选频率') axis( ) subplot(2,3,3) stem(F{1,4})% title('体温舒张压相干函数Ctd入选频率') axis( ) subplot(2,3,4) stem(F{2,3})% title('脉搏收缩压相干函数Cmg入选频率') axis( ) subplot(2,3,5) stem(F{2,4})% title('脉搏舒张压相干函数Cmd入选频率') axis( ) subplot(2,3,6) stem(F{3,4})% title('收缩压舒张压相干函数Cgd入选频率') axis( ) 运行，得：图24-6 figure(7) subplot(2,3,1) stem(F{1,2})% title('体温脉搏相干函数Ctm入选频率') axis( ) subplot(2,3,2) stem(F{1,5})% title('体温均压相干函数Ctj入选频率') axis( ) subplot(2,3,3) stem(F{1,6})% title('体温差压相干函数Ctc入选频率') axis( ) subplot(2,3,4) stem(F{2,5})% title('脉搏均压相干函数Cmj入选频率') axis( ) subplot(2,3,5) stem(F{2,6})% title('脉搏差压相干函数Cmc入选频率') axis( ) subplot(2,3,6) stem(F{5,6})% title('均压差压相干函数Cjc入选频率') axis( ) 运行，得：图24-7 下面将上述各频率向量串合在一起。此串合向量Vf取倒数就得到周期数据，以后可以将其并入《21、人体节律周期初步统计排查》的周期数据向量中一起作统计。 Vf= ; figure(8) stem(Vf) title('各相干函数入选频率串合') 运行，得：图24-8 为了保持信息来源的均衡性，收缩压、舒张压对均压、差压的相干函数分析我都没有再作，而加入统计的体温对脉搏的频率向量F{1,2}数加倍。频率数据串合向量Vf现在就可以作一下频数统计： =sort(Vf); figure(9) = hist(Bf, ); subplot(1,2,1) barh(xoutf,nf) title('所有入选频率的频数统计') ylim( ) subplot(1,2,2) bar(Bf) title('所有入选频率的排序图') axis( ) 运行，得：图24-9 这个方法可以最大限度消除传感器噪声谱，但是不能消除谐波频率。（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de01012lby.html 首发时间：2012-02-04 11:45:16）; 个人分类: 斤斤计较|2848 次阅读|0 个评论

21、人体节律周期初步统计排查: baishp 2012-11-29 23:08; 上一篇博文求出了各生理信号序列imf分量的平均周期，但是靠单一的周期数据是难以确定它就是人体节律周期的。对已经求出的那些周期数据，单凭眼睛也很难判断2个及以上的数据是否为同一周期数据。这就需要将那些数据统一放在一起并按大小次序排列起来进行观察。 VT= ;%将上篇博文中各分量周期时辰数连接为一个向量。由于人体节律周期在各通道信号中只表现为有无、强弱，而其数值是不变的，也不随信号处理方式的不同而改变，所以可以将更多不同处理方式中所获得的周期信息放在一起进行比较。这一篇博文中对信号序列再作两种处理，并将所得到的周期数据与上一篇博文中的周期数据合并观察。（一）先看看信号序列自相关函数的imf分量平均周期（自相关函数图形参见第三篇博文）。 m函数文件： function = 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= ; %脉搏 TydhV_xMB = TydhV_imfxcv( MB486012 ) VT= ; %收缩压 TydhV_xGY = TydhV_imfxcv( GY486012 ) VT= ; %舒张压 TydhV_xDY = TydhV_imfxcv( DY486012 ) VT= ; %均压 TydhV_xJY = TydhV_imfxcv( JY486012 ) VT= ; %差压 TydhV_xCY = TydhV_imfxcv( CY486012 ) VT= ; （二）再看看信号序列平方函数的imf分量平均周期 function = 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= ; %脉搏 TydhV_MB2 = TydhV_imfx2( MB486012 ) VT= ; %收缩压 TydhV_GY2 = TydhV_imfx2( GY486012 ) VT= ; %舒张压 TydhV_DY2 = TydhV_imfx2( DY486012 ) VT= ; %均压 TydhV_JY2 = TydhV_imfx2( JY486012 ) VT= ; %差压 TydhV_CY2 = TydhV_imfx2( CY486012 ) VT= ; stem(VT) 运行，得：图21-1.各信号序列及其自相关函数、平方函数的imf分量的平均周期串联图 %将周期数据向量VT按从小到大的次序排列： =sort(VT);%index是各数据在原向量VT中的序号 logb=log10(b);%取对数 stem(logb) 运行，得：图21-2.各信号序列及其自相关函数、平方函数的imf分量的平均周期排序图。上图纵坐标是周期数据的10的幂指数。从这张图中就可以看出人体生理信号很明显的节律性了。图中的台阶，表示从不同的方法处理不同通道信号所得周期数据的聚集性。在低周期端尤其明显。这是因为，对于一定长度的实验数据来说，周期越短，参见统计的波形数量就越多，统计就越精确。随着实验数据的进一步增加与信号处理方式的增加，我相信长周期端的曲线台阶也会变得越来越明显。如果人体生理没有节律性，得到的周期数据就只是些随机数，没有聚集性。其周期排序图将是一条很光滑的曲线（直线）。 stem(index) 运行，得：图21-3.各周期数据原序号通过周期数据原序号索引可以回查排序后每一个周期数据的来源（何种处理方式处理何种信号所得数据）。这个不多说了。 %周期数据频数统计 = hist(logb, ); 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= ; 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）; 个人分类: 斤斤计较|2789 次阅读|0 个评论

20、(二)生理信号imf分量平均周期的估计与分析: baishp 2012-11-29 23:04; 这一篇来看看其它生理信号imf分量的平均周期。（一）先看看脉搏信号序列MB291712的imf分量imf_MB291712的平均周期： M函数文件： function = TydhV_imf( imf_MB291712 ) P=cell(size(imf_MB291712,1)-1,2); for i=1:size(imf_MB291712,1)-1 = findpeaks(imf_MB291712(i,:)); if length(P{i,2})1 T_peak(i)=(P{i,2}(end)-P{i,2}(1))/(length(P{i,2}-1));%最后一个峰值点 %序号减去第一个峰值点序号，除以两峰值点之间峰峰段的段数。 end end B=cell(size(imf_MB291712,1)-1,2); for i=1:size(imf_MB291712,1)-1 = findpeaks(-imf_MB291712(i,:)); B{i,1}=-B{i,1}; if length(B{i,2})1 T_bottom(i)=(B{i,2}(end)-B{i,2}(1))/(length(B{i,2}-1));%最后一个谷值点 %序号减去第一个谷值点序号，除以两谷值点之间谷谷段的段数。 end end for i=1:min(length(T_peak),length(T_bottom)) T_m_MB(i)=(T_peak(i)+T_bottom(i))/2; end if length(T_peak)length(T_bottom) T_m_MB(length(T_peak))=T_peak(end); else T_m_MB(length(T_bottom))=T_bottom(end); end for i=1:length(T_m_MB) Var_MB(i)=var(imf_MB291712(i,:)); end subplot(2,1,1) bar(T_m_MB) title('T _ m _ MB') grid subplot(2,1,2) bar(Var_MB) title('Var _ MB') grid %将周期时辰数转换成年、日、小时表达方式。 T2_m_MB=T_m_MB*2; y=fix(T2_m_MB/8765.8128); d0=rem(T2_m_MB,8765.8128); d=fix(d0/24); h=rem(d0,24); TydhV_MB= ; TydhV_MB=TydhV_MB'; end 以文件名“TydhV_imf”保存(TydhV－周期、年、日、时、方差）。 TydhV_MB = TydhV_imf( imf_MB291712 ); 运行，得：图20-1 imf_MB291712各分量的平均周期与方差只画了各分量的周期与方差Bar图，上篇中其余图省去。变量TydhV_MB截图如下：图20-2 imf_MB291712各分量平均周期的年、日、小时数上图第一列是各分量周期时辰数，2、3、4列是年、日、小时数，最后一列是方差数。这些数据现在保存在这里，不作详细分析。只有第10个分量的周期136天14.55小时我觉得很有意思。因为在开通此博客之前，它在我的研究中曾经出现过。那是在我对体温信号序列作循环平稳性研究时出现的。本来，以前画出图形来，自己看看也就罢了，没有特意保存，程序也是随意处置（因为出于方便的原因，程序是处于不断改动当中，所以觉得保存也没什么意义）。但是这个136.5天的分量周期，在跟一个朋友的QQ聊天中留下了记录。见下面截图：图20-3 用不同的方法处理不同通道的信号，得到了同一结果，那么这个结果就可以确定为研究对象本身的属性，而不是由于数学处理方式的不当，凭空产生的一个虚假信息。从这个意义上来说，这里保存的处理结果，将来都有可能被证明是人体生命本身的一种属性。因此，不能轻易认为某个处理结果“没有意义”。（二）下面看看收缩压信号序列GY291712的imf分量imf_GY291712的平均周期： TydhV_GY = T TydhV_imf ( imf_GY291712 ); 运行，得：图20-4 imf_GY291712各分量的平均周期与方差变量TydhV_GY截图如下：图20-5 imf_GY291712各分量平均周期的年、日、小时数（三）下面看看舒张压信号序列DY291712的imf分量imf_DY291712的平均周期： TydhV_DY = TydhV_imf ( imf_DY291712 ); 运行，得：图20-6 imf_DY291712各分量的平均周期与方差变量TydhV_DY截图如下：图20-7 imf_DY291712各分量平均周期的年、日、小时数（四）下面看看均压信号序列JY291712的imf分量imf_JY291712的平均周期： TydhV_JY = TydhV_imf ( imf_JY291712 ); 运行，得：图20-8 imf_JY291712各分量的平均周期与方差变量TydhV_JY截图如下：图20-9 imf_JY291712各分量平均周期的年、日、小时数（五）下面看看差压信号序列JCY291712的imf分量imf_CY291712的平均周期： TydhV_CY = TydhV_imf ( imf_CY291712 ); 运行，得：图20-10 imf_CY291712各分量的平均周期与方差变量TydhV_CY截图如下：图20-11 imf_CY291712各分量平均周期的年、日、小时数（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de010126ij.html 首发时间：2012-01-07 13:55:40）; 个人分类: 斤斤计较|3514 次阅读|0 个评论

13、EMD分解_生理信号的HHT变换（一）: baishp 2012-11-28 23:17; 在我十多年前刚开始测量生理信号并思考如何样来研究的时候，那时候就想到多个复变函数合成的模型。因为看传统文化的书籍，中医或易经易学，牵涉到定量部分的内容，思考一下就知道，千言万语，无非讲的就是一个高度非线性、非平稳复变函数合成的问题（当然同时含有模式识别的问题）。只是古代文化环境、用词，与现代科技术语天差地别，所以初看不容易搞清楚事情的本质。后来有了matlab软件，知道有一个Hilbert变换函数，就对信号数据作了一下Hilbert变换，想看看变换得到的解析信号的相位与自己身体种种感觉、症状的关系。结果与实际相差太大，只能放弃。近年来现代信号处理学科高速发展。其中基于EMD（经验模式分解）的 Hilbert-Huang变换（HHT），由于能够分解出具有物理意义的分量，较之以往只能在数学框架内实现分解的信号处理方式，是一种本质性的飞跃，因此现在在很多领域都得到了很好的应用。它也可以说是最符合我事先设想的一项信号处理技术。现在就来研究一下生理信号的Hilbert-Huang变换。将我2010年10月10日～2011年6月11日共245天的体温、血压及脉搏观测数据，合并以前的观测数据整理，得到1998年2月20日～2011年6月11日共4860天的体温信号序列TW486012、2003年6月17日～2011年6月11日共2917天的收缩压信号序列GY291712、舒张压信号序列DY291712、脉搏信号序列MB291712（都没有0均值化）。绘制一下各信号序列图： subplot(2,2,1)%绘制子图 plot(TW486012) subplot(2,2,2) plot(GY291712) subplot(2,2,3) plot(DY291712) subplot(2,2,4) plot(MB291712) 运行，得：图13-1 至2011年6月11日为止的各信号序列图下面先对各信号序列作EMD分解。这里说明一下，我用的HHT-EMD工具箱，是网上已经得到广泛认可并得到广泛流传的、由G.Rilling编写的HHT-EMD函数工具箱。（1）对体温序列TW486012作EMD分解 imf_TW486012=emd(TW486012); 运行后得到一个15*58320的矩阵imf_TW486012。这表示TW486012分解成了14个基本模式分量（或称“本征模函数”，intrisic mode function,简称imf）与1个残余项。58320是信号序列（时间）长度。看看各IMF分量图： %TW486012的IMF分量图 for k=1:15 subplot(5,3,k) plot(imf_TW486012(k,:)) end 运行，得：图13-2 体温序列TW486012的IMF分量（2）对收缩压序列GY291712作EMD分解 imf_GY291712=emd(GY291712); 运行后得到一个12*35004的矩阵imf_GY291712。这表示GY291712分解成了11个基本模式分量（IMF）与1个残余项。35004是序列（时间）长度。 %GY291712的IMF分量图 for k=1:12 subplot(4,3,k) plot(imf_GY291712(k,:)) end 运行，得：图13-3 收缩压信号序列GY291712的IMF分量（3）对舒张压序列DY291712作EMD分解 imf_DY291712=emd(DY291712); 运行后得到一个13*35004的矩阵imf_DY291712。这表示DY291712分解成了12个基本模式分量（IMF）与1个残余项。35004是序列（时间）长度。 %DY291712的IMF分量图 for k=1:13 subplot(5,3,k) plot(imf_DY291712(k,:)) end 运行，得：图13-4 舒张压序列DY291712的IMF分量（4）对均压序列JY291712作EMD分解 JY291712=(GY291712+DY291712)/2; CY291712=(GY291712-DY291712)/2; 看看均压JY291712、差压CY291712的基本模式分量。之所以总是要把这两个变量数据拿来研究，目的还是要找出影响血压变化的最直接的因素。 imf_JY291712=emd(JY291712); 运行后得到一个14*35004的矩阵imf_JY291712。这表示JY291712分解成了13个基本模式分量（IMF）与1个残余项。35004是序列（时间）长度。 %JY291712的IMF分量图 for k=1:14 subplot(5,3,k) plot(imf_JY291712(k,:)) end 运行，得：图13-5 均压序列JY291712的IMF分量（5）对差压序列CY291712作EMD分解 imf_CY291712=emd(CY291712); 运行后得到一个14*35004的矩阵imf_CY291712。这表示CY291712分解成了13个基本模式分量（IMF）与1个残余项。35004是序列（时间）长度。 %CY291712的IMF分量图 for k=1:14 subplot(5,3,k) plot(imf_CY291712(k,:)) end 运行，得：图13-6 差压序列CY291712的IMF分量（6）对脉搏序列MB291712作EMD分解 imf_MB291712=emd(MB291712); 运行后得到一个14*35004的矩阵imf_MB291712。这表示MB291712分解成了13个基本模式分量（IMF）与1个残余项。35004是序列（时间）长度。 %MB291712的IMF分量图 for k=1:14 subplot(5,3,k) plot(imf_MB291712(k,:)) end 运行，得：图13-7 脉搏序列MB291712的IMF分量（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de0100zhfb.html 首发时间：2011-09-13 12:21:54）; 个人分类: 斤斤计较|5356 次阅读|0 个评论

5、互相关函数估计: baishp 2012-11-27 18:33; 这一篇来看看各信号序列体温TW267212_0，高压GY267212_0，低压DY267212_0，(均压JY267212_0，差压CY267212_0)，脉搏MB267212_0的互相关函数图。 %体温与高压的互相关函数图 RR_TWGY2672=xcorr(TW267212_0,GY267212_0,'unbiased'); plot(RR_TWGY2672) title('TW267212_0,GY267212_0的互相关函数RR_TWGY2672图') 运行,得: 图5-1 TW267212_0,GY267212_0的互相关函数RR_TWGY2672 将RR_TWGY2672导入信号处理工具sptool，打开信号观察器窗口Signal Browser，并放大横轴128倍，纵轴3.56倍,得图5-2。图5-2 RR_TWGY2672横轴放大128倍,纵轴放大3.56倍可以看到其波峰位置大致可以分为两类，一类位于图形上包络线上,一类位于上下包络线之间且靠近下包络线。RR_TWGY2672的长度是32064×2-1＝64127，但信号观察器窗口Signal Browser信号第一个数据横坐标为“0”，因此图中左标尺的位置32063，实际上就是时延为0时的互相关函数值。互相关函数为什么在时延为0处也有一个特别大的峰值呢？很奇怪。图中波峰位于上包络线上的波形周期为dx/35=418/35=~11.9429。可见此波形周期并不必然等于12即1天。这是为什么呢？波峰位置位于上下包络线之间，我大致看了一下，这种现象主要位于曲线中部。在曲线两边，我看到的波峰位置都在上包络线上。见图5-3。图5-3 RR_TWGY2672标尺换位图中波形周期为dx/32=385/32=12.03125。再看体温与低压的互相关函数图： RR_TWDY2672=xcorr(TW267212_0,DY267212_0,'unbiased'); plot(RR_TWDY2672) 运行,得: 图5-4 TW267212_0,DY267212_0的互相关函数RR_TWDY2672 图5-5 RR_TWDY2672在中心线附近的波形放大图在信号观察器Signal Browser中放大RR_TWDY2672，可以看出一个波形的周期为dx/40=480/40=12。此图中左边的标尺位于一个明显较大的波峰上，其横坐标32064，比图5-2右移一个采样单位。图5-6 RR_TWDY2672在中心线左侧某处的波形放大图可以看出每一个波形有两个并列的波峰。这样的波形周期为dx/22=276/22=12.5455。在整个横轴上，波形由等距单个波峰变为近距并列双峰，这种变化不止一次。由于数据太长，手工查起来比较费事，就不仔细清查了。反正已经知道，在体温-高压互相关函数RR_TWGY2672与体温-低压互相关函数RR_TWDY2672中，周期在T=12左右摆动的波形，是由（至少）两个周期有微小差异（或相位摆动）的波形叠加而形成的。下面看体温与均压、差压的互相关函数图： RR_TWJY2672=xcorr(TW267212_0,JY267212_0,'unbiased'); plot(RR_TWJY2672) 运行，得：图5-7 TW267212_0与JY267212_0的互相关函数RR_TWJY2672 图5-8 RR_TWJY2672放大图其左标尺位置32063是时延等于0处，此处波峰明显较附近其它处大。此处波形周期为dx/20=241/20=12.05。 RR_TWCY2672=xcorr(TW267212_0,CY267212_0,'unbiased'); plot(RR_TWCY2672) 运行，得：图5-9 TW267212_0与CY267212_0的互相关函数RR_TWCY2672 图5-10 RR_TWCY2672放大图其左标尺位置32065是时延等于0处右移两个采样单位，此处波谷（绝对值）明显较附近其它处小。此处波形周期为dx/20=239/20=11.95。体温与血压的互相关函数图就贴到这里。下面看看高压与低压、均压与差压的互相关函数图。它们的互相关系数在上一篇中已经给出过了。 RR_GYDY2672=xcorr(GY267212_0,DY267212_0,'unbiased'); plot(RR_GYDY2672) 运行，得：图5-11 GY267212_0与DY267212_0的互相关函数RR_GYDY2672 图5-12 RR_GYDY2672放大图此图中，左标尺位置32063为时延等于0处，且此处的波峰值既较附近波峰值大很多。波峰位于上包络线上的波形周期为dx/30=360/30=12。 RR_JYCY2672=xcorr(JY267212_0,CY267212_0,'unbiased'); plot(RR_JYCY2672) 运行，得：图5-13 JY267212_0与CY267212_0的互相关函数RR_JYCY2672 图5-14 RR_JYCY2672放大图此图中，左标尺位置32063为时延等于0处，且此处的波峰值较附近波峰值大很多。此图看不出明显的周期性波形了。下面贴出均压、差压与脉搏的互相关函数图： RR_JYMB2672=xcorr(JY267212_0,MB267212_0,'unbiased'); plot(RR_JYMB2672) 运行，得：图5-15 JY267212_0,MB267212_0的互相关函数RR_JYMB2672 图5-16 RR_JYMB2672放大图 RR_CYMB2672=xcorr(CY267212_0,MB267212_0,'unbiased'); plot(RR_CYMB2672) 运行，得：图5-17 CY267212_0,MB267212_0的互相关函数RR_CYMB2672 图5-18 RR_CYMB2672放大图下面贴出体温与脉搏的互相关函数图： RR_TWMB2672=xcorr(TW267212_0,MB267212_0,'unbiased'); plot(RR_TWMB2672) 运行，得：图5-19 TW267212_0,MB267212_0的互相关函数RR_TWMB2672 图5-20 RR_TWMB2672放大图图5-21RR_TWMB2672干涉条纹图最后再看看体温、均压、差压与脉搏的互相关系数： c=corrcoef( ) 运行，得： c = 1.0000 0.0181 -0.0721 0.6069 0.0181 1.0000 0.5091 0.1765 -0.0721 0.5091 1.0000 -0.0227 0.6069 0.1765 -0.0227 1.0000 说明： TW267212_0，JY267212_0的互相关系数c_TJ=0.0181; TW267212_0，CY267212_0的互相关系数c_TC=-0.0721; TW267212_0，MB267212_0的互相关系数c_TM=0.6069; JY267212_0，CY267212_0的互相关系数c_JC=0.5091; JY267212_0，MB267212_0的互相关系数c_JM=0.1765; C Y267212_0，MB267212_0 的互相关系数c_CM=-0.0227; 基本上可以认为：体温与均压、差压都没有相关性；差压与脉搏也没有相关性。我为什么不厌其烦地在这里贴这些图、记这些数？因为某些现象我现在还不能很清楚地解释其原因，需要将来进一步的研究，如功率谱分析、时频分析……之后，再回过头来看。我把它记下来，既可提醒自己以后继续思考，也可以作为资料提供给将来其他的同道爱好者研究。我希望人们将来研究生理学、生命科学的热情，就象迄今为止人们研究物理学一样，有一事不知即为耻。作为“开放巨系统”中的一个“子系统”，这些资料是不可能再生的。时间不会倒流，空间不可重叠，其中每一张图片、每一个数据都是不可复制的。俗话又说，生死无常，生命脆弱。如果哪一天呼啦啦我的房子倒了，那么这些资料就等于是从废墟中抢救出来的“珍稀贵重物品”了。嘻！（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de0100os3d.html 首发时间：2011-01-10 15:53:09）; 个人分类: 斤斤计较|4690 次阅读|0 个评论

4、谈谈我对人体血压的看法: baishp 2012-11-27 16:26; 在观察各信号序列互相关函数之前，我想谈谈我对人体血压的看法。我认为影响血压高低，主要有两个因素：一是血液总量相对于人体血管（及其它容血脏器）容积的多少。如果血液量大，那么血压就上升，反之则下降；二是人体血管（及其它容血脏器）的柔韧性（或称老化、硬化程度）。如果血管柔韧性好，那么收缩压会下降，舒张压会上升，反之收缩压会上升，舒张压会下降。根据我这个观点，收缩压与舒张压都不能单独反映上面两个因素中的任何一个。如果将收缩压与舒张压作一下变换，就能够很好单独反映上面的两个因素。这个变换就是：取收缩压与舒张压的平均值（我简称其为“均压”（JY))、收缩压与舒张压的差值（我简称其为“差压”（CY))。均压的高低可单独反映血液量的多少，差压的大小也可单独反映血管硬化程度。读者可能会注意到，在我上一篇博文中，我的高压、低压的自相关函数图图3-5、图3-7极为相似。这说明它们的变化步调很相似。它们一定有很高的互相关性。现在来看看我的均压自相关函数图： JY267212_0=(GY267212_0+DY267212_0)/2;%均压 R_JY2672=xcorr(JY267212_0,'unbiased'); plot(R_JY2672) 运行,得: 图4-1 均压JY267212_0的自相关函数R_JY2672 图4-2 R_JY2672放大图可以看出与高压自相关函数图3-5、低压自相关函数图3-7都极为相似。从生理学角度，这是很容易理解的。再来看看我的差压自相关函数图： CY267212_0=(GY267212_0-DY267212_0)/2;%差压.为了保持此变换亦具可逆性,故亦除以2. R_CY2672=xcorr(CY267212_0,'unbiased'); plot(R_CY2672) 运行,得: 图4-3 差压CY267212_0的自相关函数R_CY2672 图4-4 R_CY2672放大图这个与高压自相关函数图3-5、低压自相关函数图3-7就不那么相似了。下面看看高压、低压与均压、差压的互相关系数： c=corrcoef( ) 运行，得： c = 1.0000 0.8556 0.9736 0.6922 0.8556 1.0000 0.9512 0.2186 0.9736 0.9512 1.0000 0.5091 0.6922 0.2186 0.5091 1.0000 说明： GY267212_0，DY267212_0的互相关系数c_GD=0.8556; GY267212_0，JY267212_0的互相关系数c_GJ=0.9736; GY267212_0，CY267212_0的互相关系数c_GC=0.6922; DY267212_0，JY267212_0的互相关系数c_DJ=0.9512; DY267212_0，CY267212_0的互相关系数c_DC=0.2186; J Y267212_0，CY267212_0 的互相关系数c_JC=0.5091; 高压、低压的互相关系数达0.8556，说明二者变化同步性强，老是单独分析高压、低压变量，有些浪费精力。同样的信息含量，均压、差压的互相关系数只有0.5091，因此在分析均压、差压时能够更容易地看到更多的东西(当然，如果将各列信号作正交变换，变换后的信号相关系数都等于0，但其生理学意义就没有了）。而且由于均压、差压的生理学意义比高压、低压的生理学意义更明确，所以本人以后在对血压数据进行分析时，可能更多地使用均压、差压的概念。虽然高压、低压都同等地为均压、差压贡献了信息量，但高压、低压与均压的互相关系数都很高，而与差压的互相关系数都很低。这说明本人过去这些年来的血压变化，主要是基于均压变化而变化，而不是基于差压的变化而变化。如果我的关于均压、差压的生理学意义的论述是对的，那就说明我的血压问题在于血液量的多少而不在于血管的硬化。血液量的多少，这个在西医名称中应该叫做“贫血”程度吧？“贫血”对应于中医名词，我觉得应该是“气血两虚”，而不是单纯的“血虚”。那些造电子血压计的厂家，应该在他们的产品上设立一个均压、差压读取按钮，让使用者常常看一下：我今天贫血了吗？我今天血管硬化了吗？厂家如果因我这一条建议而赚了大钱，记得给我一点利润分成哦！（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de0100oqzn.html 首发时间：2011-01-09 11:58:55）; 个人分类: 斤斤计较|3030 次阅读|0 个评论

3、自相关函数估计: baishp 2012-11-27 12:35; 对于做科研工作或工程应用的人来说，拿到一列或数列随机信号观测序列，首先想到的就是看看它们的自相关、互相关函数，自功率谱、互功率谱密度函数估计吧。这一篇对上一篇中各数据矩阵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中的一个周期成分呢？我不知道。留待以后慢慢思考。验证。也请内行朋友赐教。这一类的“干涉条纹”图还有很多，无法尽举。我大致数了的，周期成分还有： 516372.89354.32327.6323.5319.2279262.15262258188.11186.44163.71138.92131.14131.0881.85777.45474.18267.6665.07734.60433.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）; 个人分类: 斤斤计较|5862 次阅读|0 个评论

2、观测数据_用于研究的“原材料”: baishp 2012-11-27 02:59; （1）体温数据的原始记录。不失一般性，略去其十位数“3”及小数点，如36.85℃即记为685。（2）血压、脉搏数据的原始记录。血压值开始一段时间使用千帕kpa作单位（高压部分略去百位数“1”），后换算成并一直改用为毫米汞柱mmHg。脉搏值单位为次/分钟。（3）将体温数据经3次样条插值处理并整理后，得到一个4615×12的体温数据矩阵TW4615_12。其中4615来源于1998、2、20～2010、10、9之间的天数，12为每天的插值点数（自临晨1点至晚23点，平均每2个小时取一个插值点）。（4）将血高压（为简便计，我常将收缩压称为“高压”）数据经3次样条插值处理并整理后，得到一个2672×12的血高压数据矩阵GY2672_12。其中2672来源于2003、6、17～2010、10、9之间的天数，12为每天的插值点数（自临晨1点至晚23点，平均每2个小时取一个插值点）。（5）将血低压（为简便计，我常将舒张压称为“低压”）数据经3次样条插值处理并整理后，得到一个2672×12的血低压数据矩阵DY2672_12。其中2672来源于2003、6、17～2010、10、9之间的天数，12为每天的插值点数（自临晨1点至晚23点，平均每2个小时取一个插值点）。（6）将脉搏数据经3次样条插值处理并整理后，得到一个2672×12的脉搏数据矩阵MB2672_12。其中2672来源于2003、6、17～2010、10、9之间的天数，12为每天的插值点数（自临晨1点至晚23点，平均每2个小时取一个插值点）。（7）截取矩阵TW4615_12后面自2003、6、17～2010、10、9日的部分，得到一个2672×12的体温数据矩阵TW2672_12。此矩阵数据用来与GY2672_12、DY2672_12、MB2672_12作横向比较。（本文首发网址： http://blog.sina.com.cn/s/blog_6ad0d3de0100oofr.html 首发时间：2011-01-04 20:03:01）; 个人分类: 斤斤计较|4582 次阅读|0 个评论

1、研究缘起_从人体生理信号开始研究，直趋宇宙天地人生大奥秘: baishp 2012-11-27 00:40; 寒来暑往，斗转星移，宇宙天地的运动就是律动。人生于天地之间，无时无刻不在与天地进行着节律互动。其实说“互动”只是一种方便的说法，严格地说，人体与宇宙天地都受着同样的一种力量支配，呈现着相关性很强的节律。人们最了解、最不容置疑的一个事实，就是女人们每月一次的例事。如果我们眼睛看到了，身体感觉到了，我们就承认，否则就不承认，那真是愚蠢之极。因为我们眼力有限，感觉也迟钝。我们平时看不到、感觉不到的现象，到了某个时候却会变成撕心裂肺、天崩地裂的大事。所以我们应该尽量全面地了解人体节律。节律是我们生老病死的主宰，也是我们研究生老病死的枢纽。根据我们传统文化里的“干支”理论，人体有一部分的节律周期，与我们肉眼看得到的宇宙天体运行周期相同，如地支中的“月支”与“时支”；有一部分的节律周期，却无法从宇宙天体运行周期中得出，如天干、“年支”与“日支”。人体到底有哪些节律周期呢？这仍是一个谜。即便现有的干支理论，其实也是尚未得到现代科学意义上证实的事情。只是老祖宗留给我们这么个东西，我们就用了。中医医理。命理八字，学得好的人，用起来十中八九也不是极罕见的事。但另一方面，我们却从来没有看到过一个百发百中的活神仙。这说明我们身体的节律，始终有一部分还未被人了解、认识。在现代，涉及到中医、人体生命科学、人体节律研究的，可以大致分为两大阵营。一是在现代医学、系统科学知识环境下成长起来的人，可称之为“科班”出身的人，包括一些很有名的大科学家（如钱学森先生）。这一部分人又始终停留在两个极端：要么就是宇宙开放巨系统、复杂性科学、混沌理论……这些乍听起来玄乎得很的理论。也不知他们到底要拿什么来做研究的“原材料”。你总得要有观测数据吧？没有数据的研究，只能永远停留在哲学思辩上。这些“科班”人士要么就是停留在另一个极端，眼睛看着具体一时的个体，几秒钟、几分钟的心电图、脑电图……而至今看不到这两个极端互相靠拢的一丝丝迹象。这些“科班”人士还有一个共同的特点，就是对传统文化总是无法做到完全的开放，完全的兼容并蓄。说先入为主也好，思想桎梏也罢，都可以。记得有位学者，对传统文化与中医都很推崇，但他在一次演讲中却说“八字是迷信”，令人唏嘘不已（因为最近没有搜到这段话了，故略去其名）。岂不知现在的子平八字，它最初的理论根据正是源于中医始祖典籍《黄帝内经》。另一阵营的人，就是民间研究命理八字、易经易学、医理医易的这些人，也可称为草根人士。这些人思想开放、活跃，可惜大多没受过现代科学知识的熏陶，无法用现代科学的手段来深入研究他们极其感兴趣的这门学问。不过也正是这些人，延续、保存了我们传统文化最宝贵的这一部分内容，功为至巨。这两个阵营，如果能够互相接头、融合，对中医与生命科学的研究，必将产生本质性的、飞跃性的推动作用。南怀瑾先生在讲解《黄帝内经》时指出：《黄帝内经》的宗旨，最重要者是“举痛论篇”中所说的三要义：“（一）善言天者，必有验於人。（二）善言古者，必有合於今。（三）善言人者，必有厌於己。如此则道不惑而要数极，所谓明也。”意思凡是谈到天地之道，都必定会在人类身上反映出来；凡是古时的天地规律，亦必适用于现代；凡是在别人身上存在的规律，也必定会存在于在自己身上。所谓“道不远人”吧。因此要研究宇宙天地之道，并不需要用高倍天文望远镜观测十万亿光年之外的宇宙深处。只要把自己研究透了，也等于把宇宙研究透了。老子“道可道，非常道”，“惚兮恍兮，其中有象；恍兮惚兮，其中有物”……他从哪里去弄高倍天文望远镜啊？也无法想象他能够用复杂巨系统理论来研究宇宙。博主由于个人命运的关系，在十多年前对宇宙天地、生命、生死之道产生了很大兴趣，阅读了一些传统文化典籍，也领悟到探求这个“道”，重点在对自身的研究；另一方面，由于博主工科出身，对传统文化（中医医理、八字命理）中既有的研究方式，也就是只有定性的论述，没有定量的研究手段，极不适应。因此打算另辟蹊径，将现代科技手段与传统文化结合起来研究。要对自身作定量的研究，除了不断测量自身的生理信号，研究这些生理信号外，还有什么途径？没有。人体生理信号有无数，但要长期、大量测量，必须是很容易获得的信号数据，必须是非侵入式测量。结合博主所能获得的测量工具，即定在体温、血压（包括收缩压与舒张压）与脉搏这三类四种信号。关于“道”，真是一个说不完的话题。这个“开篇词”就写到这里吧，以后慢慢道来。也希望以后与博友们有很好的互动，彼此帮助，共同进步。（本文首发网址： http://blog.sina.com.cn/s/blog_6ad0d3de0100oofe.html 首发时间：2011-01-04 19:49:21）; 个人分类: 夸夸其谈|1112 次阅读|0 个评论

更多...

帐号		自动登录	找回密码
密码			注册

关闭 安全验证

标签: 生理信号

相关帖子

相关日志

关闭安全验证