|||
各位未来的博友,龙年快乐!
从第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)
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 11:18
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社