小波变换检测信号突变点的MATLAB实现

以前在不经意间也有接触过求突变点的问题。在我看来，与其说是求突变点，不如说是咱们经常玩的"找不一样"。给你两幅图像，让你找出两个图像中不一样的地方，我认为这其实也是找突变点在生活中的应用之一吧。回到找突变点位置上，之前本身有过一个傻傻的方法：就是直接求先后两个采样的的差分值，最后设置一个阈值，若是差分值大于这个阈值则该点是突变点。但这个方法问题很大，实际中突变点幅值有大有小，你怎么能肯定阈值究竟是多少呢？还有可能信号原本的差分值就比你那突变点的差分值还要大。因此这种方法在信号或噪声稍微复杂一点就行不通了。

这几天看到了一种船新的信号突变点检测的方法-基于小波变换的信号突变点检测。因而乎去学习了一下小波变换的相关知识和应用，学习的不是很深刻，但也模模糊糊感受到了小波变换确实是检测突变点的一大利器，下面分为两个大部分总结一下我所学习到的小波变换求突变点的实现过程和相关知识理论。

小波变换求信号突变点实现

我喜欢直接从应用入手，或者应用结合理论。一步一步分析代码，看数据和图像的变化比一步一步推公式有趣的多(虽然多是错误的呀)。因而在这里我就先直接上代码和图像了，这样先让咱们对整个过程有个感性的认识。

原始信号的生成

首先生成原始信号，这里随便什么信号均可以，那我就生成一个正弦信号吧，具体信号参数见代码注释。

clear all; close all; clc;    
Fs = 1000;                 % 采样频率1000Hz
Ts = 1 / Fs;               % 采样时间间隔1ms
L = 1000;                  % 采样点数1000
t = (0 : L - 1) * Ts;      % 采样时间。1000个点，每一个点1ms，至关于采集了1s
x = sin(2 * pi * 10 * t);  % 原始正弦信号，频率为10Hz，振幅为1

添加突变点

第二步咱们要人为添加突变点了，为了看起来直观就暂时不添加噪声了。此处咱们添加两个突变点，将第233个点的幅度在本来基础上增长0.5，将第666个点的幅度在本来基础上增长0.1，代码和添加后信号图像以下：

x(233) = x(233) + 0.5;
x(666) = x(666) + 0.1;

能够看到一个突变点很明显，而另外一个却不是那么的明显，可能肉眼看的话都会忽略掉这个突变点。

对信号作傅里叶变换

可能有人要问，既然咱们作的是小波变换，为何又要对信号作傅里叶变换呢？其实咱们确实能够不用作傅里叶变换的，可是为了与小波变换作对比，分析各自的优点和劣势，咱们仍是看一下该突变信号的傅里叶变换。

Y = fft(x,1024);
f = Fs * (0 : (L / 2)) / L;
P2 = abs(Y / L);
P1 = P2(1 : L / 2 + 1);
plot(f,P1) 
title('突变信号的单边幅度频谱')
xlabel('f(Hz)')
ylabel('|P1(f)|')
axis([0,100,0,0.5])

补充

下面咱们再给一个原始信号的fft幅度谱来作对比。从肉眼来看，咱们能够发现原始信号和添加突变信号的频域差异不大，只是突变信号的频谱在高频部分多了些抖动。因此要从傅里叶变换的频域来检测突变信号是不合适的。具体缘由在第二部分会有总结，主要是两个变换选取“基”的不一样而致使的。

对信号作小波变换

重头戏小波变换来了，这里咱们用两种小波变换的方法检测突变点，一是连续小波变换；二是离散小波变换，这里只会简略说明一下图像，能够结合第二部分原理一块儿查看。

连续小波变换

咱们对突变信号进行连续小波变换，实现代码和图像以下：

cw1 = cwt(x,1:32,'sym2','plot'); % 对信号作连续小波变换
title('连续小波变换');

cwt(Continuous wavelet transform)函数表示进行连续小波变换，主要是处理一维的数据，好比咱们这个数据。参数x是输入的信号；1：32表示尺度参数Scales的取值范围为(1:32)；'sym2'表示咱们用的小波是sym2小波；'plot'是画出连续小波变换系数的意思。运行图像以下：

不一样于傅里叶变换只有w一个自变量，小波变换有两个自变量，分别是a(尺度参数)和b(位移参数)。从上图咱们能够看出在小波位移到第233个点和第666个点，且a较小时，能够看到一条较亮的白条，能够暂且理解成小波在这个位移和尺度上与信号相关性较大。在某个位置出现小波与信号相关性激增的缘由就是信号在这个位置出现了突变，因而咱们就愉快的找到了两个突变点的位置。

离散小波变换

因为连续小波变换的位移参数(b)和尺度参数(a)都是连续变化的，特别是尺度参数的连续变化，会带来巨大的计算量，因而利用离散小波变换来处理信号，这里仍是主要说代码和图像，具体实现原理在第二部分有粗浅介绍。

[d,a]=wavedec(x,3,'db4');           %对原始信号进行3层离散小波分解
a3=wrcoef('a',d,a,'db4',3);         %重构第3层近似系数
d3=wrcoef('d',d,a,'db4',3);         %重构第3层细节系数  
d2=wrcoef('d',d,a,'db4',2);         %重构第2层细节系数  
d1=wrcoef('d',d,a,'db4',1);         %重构第1层细节系数  

subplot(411);plot(a3);ylabel('近似信号a3');   %画出各层小波系数
title('小波分解示意图');
subplot(412);plot(d3);ylabel('细节信号d3');
subplot(413);plot(d2);ylabel('细节信号d2');
subplot(414);plot(d1);ylabel('细节信号d1');
xlabel('时间');

wavedec(wavelet decomposition)函数表示进行离散多辨小波分解，x是待处理的输入信号；3表示进行3层分解，'db4'也是一个小波的名字。处理完毕后获得一、二、3层的细节系数(details)和第3层的近似系数(Approximations)。画出这些系数的图像以下：

由上图可明显看出，除去开头和结尾的一些比较大的点外，在第一、二、3层的细节信号中，最大值点偏偏是第233点和第666点，因而也能够愉快的能够肯定这两个点便是突变信号的位置了。

这里还能够稍微注意一下近似信号a3，它相似于原始信号滤去了高频成分的样子，它是怎么得来的你看了第二部分就知道了！

总结

在这一部分中咱们直抓要害，知道了怎么经过小波变换来检测信号的突变点，MATLAB的函数用起来就是爽有木有。可是能应用是一回事，咱们仍是尽可能多了解一下小波变换的原理为好。小波是怎么构造的，它的性质有什么？连续小波变换的图像是怎么计算出来的呢?离散小波变换的每一层又是怎么算出来的呢？只有学习了它们背后的支撑运算的数学公式，咱们才能算真正理解了它。

---

小波变换的基础知识