各位未来的博友，龙年快乐！

 

    从第6篇博文《功率谱密度函数估计》中可知，生理信号（体温）的功率谱，除了周期为12时辰的基波及其谐波谱线，其余是大量的低幅波。这些低幅波真的没有任何分析价值吗？我认为不是的。本来，从一般的控制系统来说，这些低幅波都属于噪声谱（由系统噪声谱及观测噪声谱组成）。噪声谱中除了与白噪声相应的波谱外，还会有若干的周期波谱。（人体）系统谱中周期波谱就是我的研究对象，所以也可以称其为信号谱。信号谱当中包括若干基波及频率等于基波频率整数倍的谐波。我现在可以确认，除了那几条线谱，信号谱中还有很多的低幅波。这篇博文的目的，就是尝试从总谱中剔除噪声谱及信号谱中的谐波成分，只提取信号谱中的基波频率。按说这是信号处理技术中的一个基本问题，由于本人才疏学浅，怎么也找不到别人在这方面的论文，只好自己来尝试一下。

 

    先来看看从体温信号序列TW486012的功率谱中可不可以提取基波频率（暂不管信号谱、噪声谱）。

 

Nfft=2^16;%快速傅立叶算法点数
Order=10000;%参数估计法阶数，选取原则参见第6篇博文

Fs=Nfft;%采样频率取Nfft，为使P_TW序号等于圆频率(再加“1”）。

TW4860120=TW486012-mean(TW486012);%取0均值
P_TW=pyulear(TW486012,Order,Nfft,Fs);%功率谱估计
[PKS_TW,LOCS_TW] = findpeaks(P_TW);%求功率谱中所有谱峰点位置
Sf=zeros(1,length(P_TW)-1);
for i=1:length(P_TW)-1%在整个功率谱频率范围内从第二个频率

%（“1HZ”）开始搜索基波频率

    b=0;
for j=i:i:length(P_TW)-1;%搜索基波及其所有的谐波频率
    b=b+min(abs(LOCS_TW-1-j));%基波及其所有的谐波频率点距其最近

%的谱峰点的距离累加
end  
    Sf(i)=b/fix((length(P_TW)-1)/i);%上述累加对基、谐波数取平均值
end
plot(Sf)

 
    图23-1

 

    由于这条曲线Sf是用来衡量、选取基波的，我就把它叫做基波选取线或者基波选取函数吧。

    图23-2

 

 

    从这个放大的图形来看，怎样从曲线Sf中选取基波频率点，的确是件很难的事。首先，基波频率点肯定是Sf的谷值点，而且Sf的谷值点越小，属基波频率的可能性就越大。但谷值点不一定就是基波频率。怎样构造一个阈值曲线，使得在此曲线以下的Sf点被选取呢？

    图23-1，Sf的右半部分是高频带，频率全部高于(length(P_TW)-1)/2（周期小于4个采样点即8小时）。根据Sf的计算方法，功率谱P_TW右半部分的每一个峰值点处，Sf值都等于0，都符合基波频率的入选条件。但仅此一个“基波”，没有谐波。这样的高频，就我的研究范围来说，已经没有什么意义了，而且大多属于噪声，故舍去。

    从Sf的左半部看，其粗轮廓跟指数曲线极接近。因此我打算构造一个指数曲线，来拟合Sf的左下部轮廓。指数曲线参数的选择，应当满足两个条件：1，在此曲线以下，Sf的谷点数量符合预设范围；2，此曲线以下，Sf的谷点在横轴方向上的分布大致均匀。

 

SfL=Sf(1:length(Sf)/2);%取Sf线的左半部分
plot(SfL)

 

a=0.9997;
k=2.0;
d=-0.07;
u=[0:length(SfL)-1];
v1=a.^u;
v2=k*v1;
Ex=v2-v2(1)+SfL(1)+d;

hold on
plot(Ex,'.-k');

 

 

    图23-3

 

    经过几次人机交互绘图，当选取a=0.9997;k=2.0;d=-0.07时，所得指数曲线Ex（上图黑线）比较吻合SfL线轮廓。a、k、d单独变化时，见红、绿、紫线。

 

axis([1,300,2,3])%局部放大

    图23-4 前图局部放大

 

    下面看看a=0.9997;k=2.0;d=-0.07时SfL位于指数曲线Ex以下的谷点数量及其在左右两半部分的分布：

SfL_Ex=SfL-Ex;
[PKS_TWSfL_Ex,LOCS_TWSfL_Ex] = findpeaks(-SfL_Ex);%谷值变峰值
BTM_TWSfL_Ex=-PKS_TWSfL_Ex;%峰值变谷值

Sfd=find(BTM_TWSfL_Ex<0);%位于指数曲线Ex以下的SfL谷值点

n=length(Sfd)%位于指数曲线Ex以下的SfL谷值点数量

 

    运行，得:

n =

    179

    
    下面nl是SfL的左半部分位于指数曲线Ex以下的谷点数量

SfL_ExL=SfL_Ex(1:length(SfL_Ex)/2);
[PKS_TWSfL_ExL,LOCS_TWSfL_ExL] = findpeaks(-SfL_ExL);
BTM_TWSfL_ExL=-PKS_TWSfL_ExL;%峰值变谷值

Ssdl=find(BTM_TWSfL_ExL<0);
nl=length(Ssdl);

r=nl/n%SfL位于指数曲线Ex以下的谷点数量,左半部分占全局的比例

 

    运行，得:

r =

    0.469273743016760

 

    通过对指数曲线Ex的参数a、k、d的搜索计算，可以使n、r最大限度接近预设值。在一定范围内，a对Ex的影响主要体现在曲线的左半部分，随着a的增加，左半部分的基波频率点较右半部分增加的更快；k对Ex的影响主要体现在曲线的右半部分，随着k的增加，右半部分的基波频率点较左半部分减少的更快；d对Ex的影响则是全局性的，d增加，整个SfL曲线上面基波频率点增加。

    固定取d=0，只搜索a与k值。现取a=[0.999555:0.000005:0.999800],k=[1.96:0.01:2.40]，循环计算相应的n、r值。

 

d=0;
u=[0:length(SfL)-1];
Mn=zeros(45,50);
Mr=zeros(45,50);
for i=1:45
    k=1.95+0.01*i;
    for j=1:50
    a=0.999550+0.000005*j;  
        v1=a.^u;
        v2=k*v1;

        Ex=v2-v2(1)+SfL(1)+d;
       
        SfL_Ex=SfL-Ex;
        [PKS_TWSfL_Ex,LOCS_TWSfL_Ex] = findpeaks(-SfL_Ex);
        BTM_TWSfL_Ex=-PKS_TWSfL_Ex;%峰值变谷值
        Mn(i,j)=length(find(BTM_TWSfL_Ex<0));%位于指数曲线以下的SfL谷值点数

        SfL_ExL=SfL_Ex(1:length(SfL_Ex)/2);
        [PKS_TWSfL_ExL,LOCS_TWSfL_ExL] = findpeaks(-SfL_ExL);
        BTM_TWSfL_ExL=-PKS_TWSfL_ExL;%峰值变谷值
        Mr(i,j)=length(find(BTM_TWSfL_ExL<0));
        Mr(i,j)=Mr(i,j)/Mn(i,j);
    end
end 

 

    Mn、Mr就是谷点总数n、谷点左右比例r与a、k的关系函数。以下绘图。先作两张曲面图，看看Mn、Mr总体情况。

 

subplot(1,2,1) 

surf(Mn)

subplot(1,2,2)

surf(Mr)

    运行，得： 

    图23-5

 

    下面作Mn、Mr的三维等高线图：

 

figure

subplot(1,2,1)

csn=contour3(Mn,[150:10:250]);

clabel(csn,[150:10:250])

 

subplot(1,2,2)

csr=contour3(Mr,[0.4:0.02:0.6])

clabel(csr,[0.4:0.02:0.6])

 

    运行，得：

    图23-6

    这个三维等高线图看起来很不爽眼。没办法，程序就是这样。

    幸亏有平面等高线图，就是上面的三维等高线在X-Y平面上的投影。

 

figure 

subplot(1,2,1)

contourf(Mn,[150:10:250])

clabel(csn,[150:10:250])

 

subplot(1,2,2)

contourf(Mr,[0.4:0.02:0.6])

clabel(csr,[0.40:0.02:0.60]) 

 

    运行，得：

    图23-7

 

    但是这两张图不能画在一个窗口里。 要画在一个窗口里，不能使用填充色。

 

figure 

contour(Mn,[150:10:250],'b')

clabel(csn,[150:10:250])

hold on

contour(Mr,[0.40:0.02:0.60],'r')

clabel(csr,[0.40:0.02:0.60])

 

    运行，得：

    图23-8

 

    Mn、Mr的等高线中，各自满足基波选取两个预设条件之一的两根等高线，如果有交点，则此交点坐标就是Mn、Mr同时满足基波选取两个预设条件的a、k值。 

 

    由于手头一时怎么也找不到等高线函数句柄矩阵csn、csr的结构说明，所以只好人工在csn、csr中查找n=180与r=0.50两条等高线的数据子矩阵，记为csn180、csr050。作图：

 

figure
plot(csn180(1,:),csn180(2,:))
hold on
plot(csr050(1,1:101),csr050(2,1:101),'r')
plot(csr050(1,102:end),csr050(2,102:end),'r')

 

    运行，得：

    图23-9

    这图只是为了说明子矩阵csn180、csr050没有找错。由于上面两条曲线是由离散数据组成，它们的交点并不一定在数据集中，本项研究中亦无必要求此严格数学意义上的交点。下面求csr050上离csn180最近的点（可称准交点）的坐标。

 

M=zeros(size(csn180,2),size(csr050,2));

for i=1:size(csn180,2);
  for j=1:size(csr050,2);
    M(i,j)=(csr050(1,j)-csn180(1,i)).^2+(csr050(2,j)-csn180(2,i)).^2;
  end
end

[b,index]=min(M);
[bb,iindex]=min(b);

pj=iindex;%就是csr050上距csn180最近的点(准交点)的列序号

pi=index(iindex);%就是csn180上距csr050最近的点(准交点)的列序号

x_r050=csr050(1,pj)%csr050准交点x坐标
y_r050=csr050(2,pj);)%csr050准交点y坐标

x_n180=csn180(1,pi)%csn180准交点x坐标
y_n180=csn180(2,pi)%csn180准交点y坐标

 

    运行，得： 

x_r050 =

   23.7561983471074

y_r050 =

    17

x_n180 =

   23.5000000000000

y_n180 =

    17 

 

    两条曲线上的“准交点”纵坐标是一样的，横坐标略有差异。csr050上的“准交点”的a、k值分别为：

    a_r050=0.999550+0.000005*x_r050=0.999668780991736；

    k_r050=1.95+0.01*y_r050=2.12000000000000。

     

下面就新取定的参数重运行前面的程序：

a=0.999668780991736;

k=2.12000000000000;

d=0;

u=[0:length(SfL)-1]; 
v1=a.^u;
v2=k*v1;

Ex=v2-v2(1)+SfL(1)+d;

SfL_Ex=SfL-Ex;
[PKS_TWSfL_Ex,LOCS_TWSfL_Ex] = findpeaks(-SfL_Ex);%谷值变峰值
BTM_TWSfL_Ex=-PKS_TWSfL_Ex;%峰值变谷值

Sfd=find(BTM_TWSfL_Ex<0)%位于指数曲线Ex以下的SfL谷值点位置

n=length(Sfd);%位于指数曲线以下的SfL谷值点数

 

SfL_ExL=SfL_Ex(1:length(SfL_Ex)/2);
[PKS_TWSfL_ExL,LOCS_TWSfL_ExL] = findpeaks(-SfL_ExL);
BTM_TWSfL_ExL=-PKS_TWSfL_ExL;%峰值变谷值
Sfdl=find(BTM_TWSfL_ExL<0);
nl=length(Sfdl);
r=nl/n;

 

    运行，得n=182，r=0.5000。

 

     

Sfdd=LOCS_TWSfL_Ex(Sfd);%这就是此次选取的基波频率：

Sfdd=[Sfdd,[0 0]];%已知Sfdd的长度为n=182，补“0”为方便截图

TSfdd=Nfft./(Sfdd*12);%换算成周期，天数
Sfdd_T=[Sfdd;TSfdd];

Sfdd_T=Sfdd_T';

Sfdd_T=Sfdd_T(:);

Sfdd_T=reshape(Sfdd_T,46,8);变成46*8的矩阵，便于截图。 

 

    对Sfdd_T截图，得：

    图23-10

 

    上表中左边4列是频率，右边4列是周期（天数）。这里面有多少是真正的基波频率，最终还需要统计确认。这篇论文，我认为基波频率选取函数Sf的构造原理是对的，关键是阈值函数Ex的构造。这有待于更深的理论研究。

 

 

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

首发时间：2012-01-23 00:00:01）