博文

25、GUI界面一维小波消噪

已有 4386 次阅读 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}=[0,6e4,-200,275];

M_r{2,1}=MB2917_r;

M_r{2,2}='MB2917_ r';

M_r{2,3}=[0,3.6e4,-50,50];

M_r{3,1}=GY2917_r;

M_r{3,2}='GY2917_ r';

M_r{3,3}=[0,3.6e4,-50,50];

M_r{4,1}=DY2917_r;

M_r{4,2}='DY2917_ r';

M_r{4,3}=[0,3.6e4,-50,50];

M_r{5,1}=JY2917_r;

M_r{5,2}='JY2917_ r';

M_r{5,3}=[0,3.6e4,-50,50];

M_r{6,1}=CY2917_r;

M_r{6,2}='CY2917_ r';

M_r{6,3}=[0,3.6e4,-50,50];

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是各信号、残差的功率谱图，蓝线原信号，紫线消噪信号，黄线残差。