前面使用的信号序列都是原观测序列，没有经过任何预处理的。这一篇来做一下这件工作。应用小波消噪之前，先用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
运行，得：
http://s14/middle/6ad0d3degb987c4c7b6ed&690
图25-1
用emd方法消除趋势项之后的各生理信号序列
下面用小波工具箱图形用户界面GUI对上述信号进行消噪。点击GUI界面左上的Wavelet
1-D调出一维小波分析界面，导入体温序列TW4860_r。先随意选一个小波db4，分解层次选最多层12层。点Analyze按钮开始分析。显示模式Display
mode选Separate Mode（分开形式）。见下图。
http://s3/middle/6ad0d3deg78f40f04b382&690
图25-2
上图中，点More
Display
Options可以选择左右两边的显示内容。我调节的，左边第一列是原始信号TW4860_r，下面12列是各层近似信号；右边第一列是各层细节系数，第二列是重构信号，再以下是各层细节信号。点View
Axes可以将上面各图放大看。各层细节系数放大图如下。
http://s5/middle/6ad0d3degb98896643074&690
图25-3
各层细节系数
颜色越浅的系数越大。大约可见第三层细节系数较邻近层系数要大。
点图25-2中的De-noiss，打开一个消噪界面，见图25-4。
http://s16/middle/6ad0d3degb98896b07c3f&690
图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，打开残差显示界面，见下图。
http://s2/middle/6ad0d3degb98896fe1291&690
图25-5 fixed
form阈值消噪之残差
上图上方是残差（原信号与消噪信号之差，即被“消除”的噪声）。中左是残差的直方图，就是我以前说的“频数函数”图，是数据总体的概率密度函数估计。中右是频数函数的积分函数，即概率分布函数（估计）。下左图是残差自相关，下右图是残差功率谱。最下面是残差各项统计数据。
继续选择阈值方式分别为Rigorous SURE（基于Stein无偏似然估计SURE阈值）、Heuristic
SURE（启发式SURE阈值）、Minimaxi（最小极大方差阈值）。分别得到消噪界面及残差图，见下图25-6～图25-11。
http://s2/middle/6ad0d3degf3f55e8820b1&690
图25-6
rigorous sure阈值消噪
http://s13/middle/6ad0d3degb988978a195c&690
图25-7
rigorous sure阈值消噪之残差
http://s5/middle/6ad0d3degb98897d0adb4&690
图25-8
heuristic sure阈值消噪
http://s7/middle/6ad0d3degb988982371f6&690
图25-9
heuristic sure阈值消噪之残差
http://s1/middle/6ad0d3degb988987e5270&690
图25-10
minimax阈值消噪
http://s5/middle/6ad0d3degb98898a17dc4&690
图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是各信号、残差的功率谱图，蓝线原信号，紫线消噪信号，黄线残差。
http://s1/middle/6ad0d3degb9a3b4515310&690
图25-12
固定形式阈值Fixed form threshold的消噪信号功率谱
http://s7/middle/6ad0d3degb9a3b4740f86&690
图25-13
固定形式阈值Fixed form threshold的消噪残差(噪声)功率谱
从上面两图中看到一个有趣的现象：原信号看不出的一个周期成分，在它的消噪信号与残差中却同时出现了，见上面两图中标尺线左边一点点的谱线。经查，知其基波周期为24个时辰。消噪信号其3倍频、9倍频谱线也很明显，残差功率谱中则只看到3倍频。
这个该如何解释呢？我目前也没有一下子想清楚。
http://s16/middle/6ad0d3degb9a3b4b7b32f&690
图25-14
原信号与固定形式阈值Fixed form threshold消噪信号、残差(噪声)三合一功率谱
http://s3/middle/6ad0d3degb9a3b4e138d2&690
图25-15
图25-14低频端放大
http://s3/middle/6ad0d3degb9a3b5286072&690
图25-16
图25-14高频端放大
http://s3/middle/6ad0d3degb9a3b55040f2&690
图25-17
无偏似然估计Rigorous SURE阈值的消噪信号功率谱
此图中未见24时辰周期成分基、倍频谱线。
http://s8/middle/6ad0d3degb9a3b56a0e27&690
图25-18
无偏似然估计Rigorous SURE阈值的消噪残差功率谱
此图中只见24时辰周期成分基频谱线。
http://s15/middle/6ad0d3degb9a3b59885ae&690
图25-19
原信号与无偏似然估计Rigorous SURE阈值的消噪信号、残差(噪声)三合一功率谱
http://s13/middle/6ad0d3degb9a3b5ee70bc&690
图25-20
最小极大方差阈值Minimaxi的消噪信号功率谱
此图又见周期为24个时辰周期成分之基频、3倍频、9倍频谱线。
http://s14/middle/6ad0d3degb9a3b61f9e7d&690
图25-21
最小极大方差阈值Minimaxi的消噪残差功率谱
此图可见周期为24个时辰周期成分之基频、3倍频谱线。
http://s9/middle/6ad0d3degb9a3b659dee8&690
图25-22
原信号与最小极大方差阈值Minimaxi的消噪信号、残差(噪声)三合一功率谱
看看各信号在时域的波形：
plot(TW4860_r)
hold on
plot(TWden_fix,'g-')
plot(TWden_rig,'r-')
plot(TWden_fix,'k-')
xlim([0,200])
运行，得：
http://s16/middle/6ad0d3degb9a3b7dbad0f&690
图25-23
原信号及3种消噪信号时域波形比较
从上面的功率谱图及时域波形光滑性可以看出，各阈值消噪的效果，与消噪相似性指标反映的消噪效果是一致的。
感觉，消噪信号的应用，应该灵活掌握。有些信号处理方式中应该用消噪信号，另一些信号处理方式中就不一定要使用消噪信号。因为这些处理方式本身可能就含有消噪功能，用已经消噪的信号，也许会丢失一部分信息。
加载中，请稍候......
加载中…
(2012-02-24 19:37:48)