前言

小波变换专业处理时变信号！其重要用途包含：突变点检测、时频分析、信号降噪等。本文将详细介绍小波变换的这3种主要用途，借助具体例子来说明并总结相关函数的使用。

间断点检测

现实信号中的间断点是较为常见的，明显的间断点就是信号的"突跳"，反应在数学上就是该点"一阶不可导"！如何突跳很明显，我们可以肉眼识别并剔除；但是如果数据本身幅值变化就很大，那么很多小的时域的间断点凭人眼是很难发现和剔除的！此时我们可以使用小波变换，把原始分解到"小波域"去查看！间断点例子如图1所示：

存在3个间断点的原始信号.png

其生成函数为：

x = zeros(1,500);

% 待检测数据
for t = 1:500
    if (t<200)
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*t);
    elseif (t>=200) & (t<300)
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*200) + 30;
    elseif (t>=300) & (t<400)
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*200) + 10;
    else
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*200) + 10 + 250*sin(2*pi*0.003*t);
    end
end

figure(1)
plot(x); title('原始信号'); grid on;
xlabel('采样点'); ylabel('振幅');

说明：间断点的检查对于原始信号总体上本应是有规律、连续、平滑变化的！但是由于偶然且剧烈的仪器或外界的干扰，导致在某采样点发生了"突跳"！图1中间断点1和2就是这种情况。另外，如果突跳还经常发现在两个频率信号的交界处，图1中的间断点3就是这种情况。

下面我们就用简单的小波分解来检测间断点：一维haar小波1级分解

x = zeros(1,500);

% 待检测数据
for t = 1:500
    if (t<200)
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*t);
    elseif (t>=200) & (t<300)
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*200) + 30;
    elseif (t>=300) & (t<400)
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*200) + 10;
    else
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*200) + 10 + 250*sin(2*pi*0.003*t);
    end
end

% 一级分解: haar
[cA1,cD1] = dwt(x,'haar');
A1 = upcoef('a',cA1,'haar',1);
D1 = upcoef('d',cD1,'haar',1);
figure(1);
subplot(1,3,1); plot(x); title('原始信号'); grid on;
xlabel('采样点'); ylabel('振幅');
subplot(1,3,2); plot(A1); title('时域: 原始信号低频近似部分(点数一样)'); grid on;
xlabel('小波域: x和y轴无量纲');
subplot(1,3,3); plot(D1); title('时域: 原始信号高频细节部分(点数一样)'); grid on;
xlabel('小波域: x和y轴无量纲');

效果：

图2：一维haar小1级分解间断点检测结果

说明：可以看到只做了1级的小波分解，其"高频细节"部分"系数图像"就可以很好的辨识出原始信号的间断点！这里需要注意的是：我们辨识间断点看的就是无量纲的、小波域的系数图像！而不需要做系数重构回到时间域。

总结1：小波域的高频细节的系数图像，专门用来检查原始信号间断点。

下面我们再使用db4小波基做多级(3级)分解看看：

x = zeros(1,500);

% 待检测数据
for t = 1:500
    if (t<200)
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*t);
    elseif (t>=200) & (t<300)
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*200) + 30;
    elseif (t>=300) & (t<400)
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*200) + 10;
    else
        x(t) = 50.*exp(t/300).*sin(2*pi*0.01*200) + 10 + 250*sin(2*pi*0.003*t);
    end
end

% 多级(3级)分解: db4
[C,L] = wavedec(x,3,'db4');
A3=wrcoef('a',C,L,'db4',3); % 低
D3=wrcoef('d',C,L,'db4',3);
D2=wrcoef('d',C,L,'db4',2);
D1=wrcoef('d',C,L,'db4',1); % 3个高
figure(2)
subplot(2,2,1); plot(A3); title('原始信号中的低频信号成分'); grid on;
xlabel('小波域: x和y轴无量纲');
subplot(2,2,2); plot(D3); title('原始信号中的高频信号成分1'); grid on;
xlabel('小波域: x和y轴无量纲');
subplot(2,2,3); plot(D2); title('原始信号中的高频信号成分2'); grid on;
xlabel('小波域: x和y轴无量纲');
subplot(2,2,4); plot(D1); title('原始信号中的高频信号成分3'); grid on;
xlabel('小波域: x和y轴无量纲');

效果：

图3：一维db4小波基3级分解间断点检测结果

说明：看最高分解级里的高频细节图像(右下角)，由于分解的级数较多，其最高级高频细节已经基本属于"只反映间断点"！但是从整体上看很明显可以感觉到分解的次数有些过多！单纯的间断间有了"波动起伏(很像小波基函数)"！所以，针对本例的原始信号，预计做2级小波分解级足够很好检测间断点。

总结2：随着分解级数的增加，最高级高频细节的系数图像几乎只反映间断点位置！

时频分析

我们已经知道，小波变换专业处理时变信号。时变信号就是信号的频率是随着时间变化的，单纯从时域来看只能通过观察图像变紧缩或舒张来判断频率是否变化，这是不精确的！如果能绘制一张横坐标是时间/采样点、纵坐标是频率的图像，那样就可以很好的观察频率随时间的变化！这样的需求就需要做"时频分析"来实现，得到的图像叫作"时频图"。

做时频分析最常用的方法就是"短时傅里叶变换"，其内涵就是用一个"窗函数"不断与原始信号做"滑动卷积"计算。其原理和小波变换的思想十分相近，本文主要介绍用短时傅里叶变换来实现时频分析：

首先观察一个时变信号：

图4：时变信号

其实现的matlab语句为：

clear all; clc;

N = 1024;  % 总采样点
fs = 1000;          % 采样频率
tt = (0:N-1)'/fs;   % 时间刻度
% 构成信号: 注意x和tt都是列向量
f1 = 400 ; f2 = 200 ; f3 = 100 ; f4 = 50 ;
x = sin(2*pi*f1*tt).*(tt<=0.3) + sin(2*pi*f2*tt).*(tt>0.3&tt<=0.6) + ...
    sin(2*pi*f3*tt).*(tt>0.6&tt<=0.8) + sin(2*pi*f4*tt).*(tt>0.8);

figure(1);
plot(tt,x);
axis([0 max(tt) -inf inf]); 
xlabel('时间/s');
ylabel('振幅');
title('原始时域信号');

从图像和实现语句可以看出：该信号有4个频率的变换，反应在时域图像上就是图像变得紧密或舒张！对于这样简单的复合函数信号，我们确实可以直接从时域图像中判断其"频域-时间"的变化关系！但这只是因为原始信号简单。

下面我们用4种不同节点个数的窗函数(汉宁窗)来实现短时傅里叶变换！其用到的函数和实现过程如下matlab程序：

clear all; clc;

N = 1024;  % 总采样点
fs = 1000;          % 采样频率
tt = (0:N-1)'/fs;   % 时间刻度
% 构成信号: 注意x和tt都是列向量
f1 = 400 ; f2 = 200 ; f3 = 100 ; f4 = 50 ;
x = sin(2*pi*f1*tt).*(tt<=0.3) + sin(2*pi*f2*tt).*(tt>0.3&tt<=0.6) + ...
    sin(2*pi*f3*tt).*(tt>0.6&tt<=0.8) + sin(2*pi*f4*tt).*(tt>0.8);

figure(1);
plot(tt,x);
axis([0 max(tt) -inf inf]); 
xlabel('时间/s');
ylabel('振幅');
title('原始时域信号');

% 短时傅里叶变换操作:
figure(2);

% 因为tfr把频率轴限制在-0.5~0.5之间，因此恢复"真实频率"的时候要x2
h = hanning(63);  % 63个节点的汉宁窗
[tfr,t,f] = tfrstft(x,1:N,N,h);  % 参数一个都不能少 
subplot(2,2,1); 
contour( tt, f(1:N/2)*fs*2, abs(tfr(1:N/2,:)) );  % 等值线图
xlabel('时间/s');  ylabel('频率/Hz'); title('63节点汉宁窗时频图');
grid on;

h = hanning(127); % 127个节点的汉宁窗
[tfr,t,f] = tfrstft(x,1:N,N,h);  % 参数一个都不能少 
subplot(2,2,2); 
contour( tt, f(1:N/2)*fs*2, abs(tfr(1:N/2,:)) );  % 等值线图
xlabel('时间/s');  ylabel('频率/Hz'); title('127节点汉宁窗时频图');
grid on;

h = hanning(255); % 255个节点的汉宁窗
[tfr,t,f] = tfrstft(x,1:N,N,h);  % 参数一个都不能少 
subplot(2,2,3); 
contour( tt, f(1:N/2)*fs*2, abs(tfr(1:N/2,:)) );  % 等值线图
xlabel('时间/s');  ylabel('频率/Hz'); title('255节点汉宁窗时频图');
grid on;

h = hanning(511); % 511个节点的汉宁窗
[tfr,t,f] = tfrstft(x,1:N,N,h);  % 参数一个都不能少 
subplot(2,2,4); 
contour( tt, f(1:N/2)*fs*2, abs(tfr(1:N/2,:)) );  % 等值线图
xlabel('时间/s');  ylabel('频率/Hz'); title('511节点汉宁窗时频图');
grid on;

效果：

图5：4种窗函数节点的时频图

说明1：从结果图中可以很明显观察"频率-时间"的变化关系！但是我们注意到窗函数节点不同时，时频图的样式有些差别！当节点/窗口太小时，纵坐标频率分辨率不够(有重叠)！当节点/窗口太大时，横坐标时间分辨率不够(有重叠)！所以不同信号窗函数大小的选取是一个值得考量的问题。但是不管选的好不好，视频图都可以很好反应频率-时间变化关系。

说明2：这个例子的分段函数信号很简单，就是单纯的4个差别较大的频率信号的组合！所以反应在时频图上就是这种"直线阶梯间隔"的样式！对于现实的信号不会用这么明显的间隔区分，时频图各段更多的是连在一起的！就像一个完整连续的xy函数图像一样。时频图的主要目的就是为了观察频率随时间的变换关系。

说明3：等值线图中有反应了3条信息：时间、频率、振幅！

短时傅里叶变换用的函数语法总结，首先总结变换类的函数：tfrstft

N = 1024;

h = hanning(63);  % 63个节点的汉宁窗
[tfr,t,f] = tfrstft(x,1:N,N,h);  % 参数一个都不能少 

% hanning(N)：汉宁窗窗口大小N的设置；
% tfrstft参数：x是原始信号，1:N是采样点数组，N是总点数，h就是上面设置的汉宁窗；
% 左边返回值：tfr是短时傅里叶变换系数，t为系数tfr对应的时刻，f为归一化频率向量；

等值线画图工具：contour

[tfr,t,f] = tfrstft(x,1:N,N,h); 

contour( t, f(1:N/2)*fs*2, abs(tfr(1:N/2,:)) );
% contour参数：
% t是原始信号的时间；
% f(1:N/2)*fs*2归一化频率：根据奈奎斯特定律只取前一半；记得*2才回到原始频率;
% abs(tfr(1:N/2,:))多时傅里叶变换系数：同样的也是只取前一半。

总结3：时频分析的目的是观察原始信号中频率随时间的变换规律，窗口大小是个较为关键的参数。

---

关于时频分析，这里再补充一个实际的案例。案例程序和原始数据这里下载

地震记录到的实际有效数据为：

图6：时频分析实际数据原始信号

用等值线图展示的时频分析：

图7：等值线图展示

mesh二维图像展示：

图8：mesh二维图像展示

mesh三维图像展示：

图9：mesh三维图像展示

小波去噪

噪声具有随机性，很难完全剔除！小波再牛逼，依旧无法除干净！这是前言。最为常用的小波去噪方法为"阈值去噪"，该方法的去噪思路为：噪声一般为高频成分，因此设定一个阈值，让分解得到的各个高频细节部分和这个阈值做对比：超过阈值的高频系数认为是噪声并把值设为0(删除)，不超过则认为是有效信号值不变(保留)。 该方法大体可分为2部分：1. 阈值选取；2. 小波去噪。这2部分matlab都有自带的函数，先列举如下：

• 小波阈值获取函数：ddencmp、thselect、wbmpen等;

• 小波去噪的函数：wden、wdencmp、wthresh等；

后面会具体介绍这些自带函数的使用。"阈值去噪"法又可分为下面的3种不同的"流程"来实现，它们将借助上面提到的函数：

• 强制去噪处理：该流程是将小波分解后的"高频系数"全为0！即滤掉所有的高频成分！然后再做"重构/恢复原始信号"。这种方法操作简单，处理后的信号很平滑，因此很容易缺失部分有用信号。

• 默认阈值去噪处理：该流程先利用ddencmp(阈值获取函数)生成信号的"默认阈值"，然后利用wdencmp(去噪函数)进行去噪处理。用wden函数可以把两步归为一步完成，即wden与wdencmp其实是一样的，只不过wden更全自动一些(一步把"取默认阈值"和"去噪"全在它自己里做了)。

• 人为给定阈值去噪处理：实际问题中阈值可以通过经验公式获得，也就是说阈值是可以人为给定的！人为给定的阈值要比第2种流程中默认的阈值可信度高。若阈值人为给定，后面去噪函数选wthresh。注意：对每一个高频成分(系数)都要手动给阈值，最后再"重构/恢复原始信号"。

注意1：默认法法中，可以用ddencmp + wdencmp搭配的方式完成；也可以用wden一步完成。最常用的是前一种。

注意2：默认法它的每一级分解的阈值都是固定的！但是人为法可以对不同级的高频系数给不同的阈值！完全的自由。

注意3：默认法2与人为法3最大的区别在于阈值是谁给的。那么介于2和3之间的，就是用thselect函数来选择性的给定阈值，该函数内置了多种选项供人去选(都是很常用的，可以就当做经验公式)，但每种算法还是固定的，所以是介于2和3之间的一种"给阈值"的方式。

下面我们将结合具体的例子要使用上面的3种流程和提到的各种自带函数：

（1）上面所有方法综合在一起处理噪声信号的案例：

% 去噪函数2: wdencmp  (和wden是一样的)  
% wdencmp中的参数很多时ddencmp给的！例如: thr、sorh、keepapp

clear; clc;

% 加载信号
load leleccum;
% 采样点
s = leleccum(300:2000);

% 原始信号图像:
figure(1);
plot(s);
xlabel('采样点'); ylabel('振幅');
grid on;
title('原始信号');

% 对原始信号做小波的3级分解, 并做系数提取:
[c,l] = wavedec(s,3,'db4');
a3 = appcoef(c,l,'db4',3);  % 低频近似部分; 3代表层数
d3 = detcoef(c,l,3);
d2 = detcoef(c,l,2);
d1 = detcoef(c,l,1);

figure(2);
% 强制去噪处理: 无需任何新的函数
%  高频"系数"全搞成0, 然后再重构
dd3 = zeros(1,length(d3));
dd2 = zeros(1,length(d2));
dd1 = zeros(1,length(d1));
c1 = [a3 dd3 dd2 dd1];
s1 = waverec(c1,l,'db4');  % 重构
subplot(2,2,1);
plot(s1);
xlabel('采样点'); ylabel('振幅');
grid on;
title('强制去噪处理');

% 默认阈值去噪处理: 用ddencmp给出信号的默认阈值, 然后用wdencmp做去噪
[thr,sorh,keepapp] = ddencmp('den','wv',s);
s2 = wdencmp('gbl', s, 'db4', 3, thr, sorh, keepapp);
subplot(2,2,2);
plot(s2);
xlabel('采样点'); ylabel('振幅');
grid on;
title('默认阈值去噪处理: wdencmp函数');

% 默认阈值去噪处理2: 用wden函数(其实和wdencmp是一样的, 只不过wden更全自动一些)
s3 = wden(s,'minimaxi','s','mln',3,'db4');
subplot(2,2,3);
plot(s3);
xlabel('采样点'); ylabel('振幅');
grid on;
title('默认阈值去噪处理: wden函数');

% 人为给定阈值去噪处理: 阈值纯人为给, 用wthresh函数: 一般是对每一个高频成分(系数)都做阈值处理
softd1 = wthresh(d1, 's', 10);  % 10是纯人给的阈值
softd2 = wthresh(d2, 's', 10);  
softd3 = wthresh(d3, 's', 25);
c2 = [a3 softd3 softd2 softd1];
s4 = waverec(c2,l,'db4');  % 小波重构
subplot(2,2,4);
plot(s4);
xlabel('采样点'); ylabel('振幅');
grid on;
title('人为给定阈值去噪处理');

figure(1)原始信号：

图9：带噪声的原始信号

figure(2)4种处理方法效果图：

图10：4种阈值去噪方式的效果图

说明：可以看出4种阈值去噪方式的效果都很好！都可以把原信号中的"紧密"部分给有效的去掉！这也正好说明原信号的噪声都是集中早高频部分，即高频噪声反应在图像上就是图像突然变的紧密。另外：两种默认阈值去噪的结果几乎是一样的，其实参数给定一样它们的结果就是一样的！再重复一遍上面的注意1：默认法法中，可以用ddencmp + wdencmp搭配的方式完成；也可以用wden一步完成。两种是完全一样的。

下面我们就来总结例子中用到函数的语法：

• 默认阈值去噪处理1：ddencmp + wdencmp

% 信号默认阈值提取：
[thr, sorh, keepapp, crit] = ddencmp(IN1, IN2, X);

% ddencmp参数：
% IN1取'den'或'cmp'：'den'表示进行去噪，'cmp'表示进行压缩；
% IN2取'wv'或'wp'，'wv'表示选择小波，'wp'表示选择小波；
% X为原始信号；

% 左边返回值：
% thr：阈值；
% sorh(soft or hard)：软阈值或硬阈值选择参数；
% keepapp：是否保存低频信号：0或1；
% crit：熵名(只在选择小波包时使用)。

% 本案例中的形式为：
[thr, sorh, keepapp] = ddencmp('den', 'wv', s);

% 小波去噪处理：
s2 = wdencmp('gbl', X, 'wname', N, thr, sorh, keepapp);

% wdencmp参数：
% 'gbl'(gloabl)：每一层都采用同一个阈值进行处理；
% X：原始信号；
% 'waname'：小波基函数；
% N：小波分解的层数；
% thr、sorh、keepapp：ddencmp的返回值；其中keepapp为1表示低频系数不进行阈值量化，为0则要做；

% 左边返回值：s2表示去噪后的结果

% 本案例中的形式为：
s2 = wdencmp('gbl', s, 'db4', 3, thr, sorh, keepapp);

• 默认阈值去噪处理2：wden

s3 = wden(X, TPRT, sorh, scal, N, 'wname');

% wden参数：
% X：原始信号；
% TPTR：阈值选取的规则，它有4种可选：
%      'rigrsure'：是一种基于Stein的无偏似然估计原理的自适应阈值选择；
%      'sqtwolog'：是一种固定的阈值形式；
%      'heursure'：前两种方式的综合，所选择的是最优预测变量阈值；适用于信噪比较小；
%      'minimaxi'：是一种根据最小方差的固定阈值选择形式。
% sorh：是选择软阈值还是硬阈值；一般都选s；
% scal：阈值尺度改变的比例，它有3种选择：
%      'one'：基本模式；可以忽略必须估计的噪声层次；
%      'sln'：非白噪声的基本模式，且在每个分解的层次上都估计噪声的层次，以此来改变阈值的尺度；
%      'mln'：不同分解层次的小波进行独立的噪声估计。
% N：分解的层次；
% 'wname'：小波基函数。

% 左边返回值：s3表示去噪后的结果

% 本案例中的形式：
s3 = wden(s,'minimaxi','s','mln',3,'db4');

补充：关于ddencmp、wdencmp、wden更细节用法参考这篇文章。

• 人为给定阈值去噪处理：wthresh

s4 = wthresh(x, sorh, T);

% wthresh参数：
% x：小波分解后的各个系数(可以是低频和高频都做)
% sorh：是选择软阈值还是硬阈值；一般都选s；
% T：纯人为给定的阈值

% 左边返回值：s4表示去噪后的结果

% 本案例中的形式：
% 先做小波分解和各个系数的提取：
[c,l] = wavedec(s,3,'db4');
a3 = appcoef(c,l,'db4',3);  
d3 = detcoef(c,l,3);
d2 = detcoef(c,l,2);
d1 = detcoef(c,l,1);
% 然后对各个系数(一般只做高频系数)的阈值(人给，可不一样)判断：
softd1 = wthresh(d1, 's', 10);  
softd2 = wthresh(d2, 's', 10);  
softd3 = wthresh(d3, 's', 25);  % 10和25是纯人给的阈值

注意：人为给定法与强制去噪法，都相对于让各级的高频系数与阈值做比较(强制法直接弄成0)，然后把"新的高频系数(保留或去除)"与低频近似部分做重构/恢复原始信号，即完成了去噪。

• 介于默认和人为给定之间的阈值给定方法：thselect

thr = thselect(X, TPTR);

% thselect参数：
% X：原始信号；
% TPTR：阈值选取的规则；

% 左边返回值：thr就是阈值

说明：人为给定阈值中的10和25都是纯人给的，这2个数值可以参考用thselect计算的几种结果。

至此，小波去噪案例和相关函数的使用语法就介绍完毕了。

---

总结

本文介绍了使用小波完成：间断点识别、时频分析(短时傅里叶变换)、信号去噪。下面把文中得到的一些有用结论再列举如下：

• 小波域的高频细节的系数图像，专门用来检查原始信号间断点；
• 随着分解级数的增加，最高级高频细节的系数图像几乎只反映间断点位置！
• 时频分析的目的是观察原始信号中频率随时间的变换规律，窗口大小是个较为关键的参数；现实情况时频图基本上都是连在一起的；
• 小波去噪用到的函数较多，并且很多时配套使用。